# source:src/libcfa/rational.c@f621a148

Last change on this file since f621a148 was f621a148, checked in by Peter A. Buhr <pabuhr@…>, 7 years ago

change to implementation type for rational and add to test suite

• Property mode set to `100644`
File size: 5.4 KB
Line
1//
2// Cforall Version 1.0.0 Copyright (C) 2016 University of Waterloo
3//
4// The contents of this file are covered under the licence agreement in the
5// file "LICENCE" distributed with Cforall.
6//
7// rational.c --
8//
9// Author           : Peter A. Buhr
10// Created On       : Wed Apr  6 17:54:28 2016
13// Update Count     : 51
14//
15
16#include "rational"
17#include "fstream"
18#include "stdlib"
19#include "math"                                                                                 // floor
20
21
22// constants
23
24struct Rational 0 = {0, 1};
25struct Rational 1 = {1, 1};
26
27
28// helper routines
29
30// Calculate greatest common denominator of two numbers, the first of which may be negative. Used to reduce rationals.
31// alternative: https://en.wikipedia.org/wiki/Binary_GCD_algorithm
32static RationalImpl gcd( RationalImpl a, RationalImpl b ) {
33        for ( ;; ) {                                                                            // Euclid's algorithm
34                RationalImpl r = a % b;
35          if ( r == 0 ) break;
36                a = b;
37                b = r;
38        } // for
39        return b;
40} // gcd
41
42static RationalImpl simplify( RationalImpl *n, RationalImpl *d ) {
43        if ( *d == 0 ) {
44                serr | "Invalid rational number construction: denominator cannot be equal to 0." | endl;
45                exit( EXIT_FAILURE );
46        } // exit
47        if ( *d < 0 ) { *d = -*d; *n = -*n; }                           // move sign to numerator
48        return gcd( abs( *n ), *d );                                            // simplify
49} // Rationalnumber::simplify
50
51
52// constructors
53
54void ?{}( Rational * r ) {
55        r{ 0, 1 };
56} // rational
57
58void ?{}( Rational * r, RationalImpl n ) {
59        r{ n, 1 };
60} // rational
61
62void ?{}( Rational * r, RationalImpl n, RationalImpl d ) {
63        RationalImpl t = simplify( &n, &d );                            // simplify
64        r->numerator = n / t;
65        r->denominator = d / t;
66} // rational
67
68
69// getter for numerator/denominator
70
71RationalImpl numerator( Rational r ) {
72        return r.numerator;
73} // numerator
74
75RationalImpl denominator( Rational r ) {
76        return r.denominator;
77} // denominator
78
79[ RationalImpl, RationalImpl ] ?=?( * [ RationalImpl, RationalImpl ] dest, Rational src ) {
80        return *dest = src.[ numerator, denominator ];
81}
82
83// setter for numerator/denominator
84
85RationalImpl numerator( Rational r, RationalImpl n ) {
86        RationalImpl prev = r.numerator;
87        RationalImpl t = gcd( abs( n ), r.denominator );                // simplify
88        r.numerator = n / t;
89        r.denominator = r.denominator / t;
90        return prev;
91} // numerator
92
93RationalImpl denominator( Rational r, RationalImpl d ) {
94        RationalImpl prev = r.denominator;
95        RationalImpl t = simplify( &r.numerator, &d );                  // simplify
96        r.numerator = r.numerator / t;
97        r.denominator = d / t;
98        return prev;
99} // denominator
100
101
102// comparison
103
104int ?==?( Rational l, Rational r ) {
105        return l.numerator * r.denominator == l.denominator * r.numerator;
106} // ?==?
107
108int ?!=?( Rational l, Rational r ) {
109        return ! ( l == r );
110} // ?!=?
111
112int ?<?( Rational l, Rational r ) {
113        return l.numerator * r.denominator < l.denominator * r.numerator;
114} // ?<?
115
116int ?<=?( Rational l, Rational r ) {
117        return l < r || l == r;
118} // ?<=?
119
120int ?>?( Rational l, Rational r ) {
121        return ! ( l <= r );
122} // ?>?
123
124int ?>=?( Rational l, Rational r ) {
125        return ! ( l < r );
126} // ?>=?
127
128
129// arithmetic
130
131Rational -?( Rational r ) {
132        Rational t = { -r.numerator, r.denominator };
133        return t;
134} // -?
135
136Rational ?+?( Rational l, Rational r ) {
137        if ( l.denominator == r.denominator ) {                         // special case
138                Rational t = { l.numerator + r.numerator, l.denominator };
139                return t;
140        } else {
141                Rational t = { l.numerator * r.denominator + l.denominator * r.numerator, l.denominator * r.denominator };
142                return t;
143        } // if
144} // ?+?
145
146Rational ?-?( Rational l, Rational r ) {
147        if ( l.denominator == r.denominator ) {                         // special case
148                Rational t = { l.numerator - r.numerator, l.denominator };
149                return t;
150        } else {
151                Rational t = { l.numerator * r.denominator - l.denominator * r.numerator, l.denominator * r.denominator };
152                return t;
153        } // if
154} // ?-?
155
156Rational ?*?( Rational l, Rational r ) {
157        Rational t = { l.numerator * r.numerator, l.denominator * r.denominator };
158        return t;
159} // ?*?
160
161Rational ?/?( Rational l, Rational r ) {
162        if ( r.numerator < 0 ) {
163                r.numerator = -r.numerator;
164                r.denominator = -r.denominator;
165        } // if
166        Rational t = { l.numerator * r.denominator, l.denominator * r.numerator };
167        return t;
168} // ?/?
169
170
171// conversion
172
173double widen( Rational r ) {
174        return (double)r.numerator / (double)r.denominator;
175} // widen
176
177// http://www.ics.uci.edu/~eppstein/numth/frap.c
178Rational narrow( double f, RationalImpl md ) {
179        if ( md <= 1 ) {                                                                        // maximum fractional digits too small?
180                return (Rational){ f, 1};                                               // truncate fraction
181        } // if
182
183        // continued fraction coefficients
184        RationalImpl m00 = 1, m11 = 1, m01 = 0, m10 = 0;
185        RationalImpl ai, t;
186
187        // find terms until denom gets too big
188        for ( ;; ) {
189                ai = (RationalImpl)f;
190          if ( ! (m10 * ai + m11 <= md) ) break;
191                t = m00 * ai + m01;
192                m01 = m00;
193                m00 = t;
194                t = m10 * ai + m11;
195                m11 = m10;
196                m10 = t;
197                t = (double)ai;
198          if ( f == t ) break;                                                          // prevent division by zero
199                f = 1 / (f - t);
200          if ( f > (double)0x7FFFFFFF ) break;                          // representation failure
201        }
202        return (Rational){ m00, m10 };
203} // narrow
204
205
206// I/O
207
208forall( dtype istype | istream( istype ) )
209istype * ?|?( istype *is, Rational *r ) {
210        RationalImpl t;
211        is | &(r->numerator) | &(r->denominator);
212        t = simplify( &(r->numerator), &(r->denominator) );
213        r->numerator /= t;
214        r->denominator /= t;
215        return is;
216} // ?|?
217
218forall( dtype ostype | ostream( ostype ) )
219ostype * ?|?( ostype *os, Rational r ) {
220        return os | r.numerator | '/' | r.denominator;
221} // ?|?
222
223// Local Variables: //
224// tab-width: 4 //
225// End: //
Note: See TracBrowser for help on using the repository browser.