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