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