// // Cforall Version 1.0.0 Copyright (C) 2016 University of Waterloo // // The contents of this file are covered under the licence agreement in the // file "LICENCE" distributed with Cforall. // // rational.c -- // // Author : Peter A. Buhr // Created On : Wed Apr 6 17:54:28 2016 // Last Modified By : Peter A. Buhr // Last Modified On : Sat Feb 8 17:56:36 2020 // Update Count : 187 // #include "rational.hfa" #include "fstream.hfa" #include "stdlib.hfa" forall( otype RationalImpl | arithmetic( RationalImpl ) ) { // helper routines // Calculate greatest common denominator of two numbers, the first of which may be negative. Used to reduce // rationals. alternative: https://en.wikipedia.org/wiki/Binary_GCD_algorithm static RationalImpl gcd( RationalImpl a, RationalImpl b ) { for ( ;; ) { // Euclid's algorithm RationalImpl r = a % b; if ( r == (RationalImpl){0} ) break; a = b; b = r; } // for return b; } // gcd static RationalImpl simplify( RationalImpl & n, RationalImpl & d ) { if ( d == (RationalImpl){0} ) { abort | "Invalid rational number construction: denominator cannot be equal to 0."; } // exit if ( d < (RationalImpl){0} ) { d = -d; n = -n; } // move sign to numerator return gcd( abs( n ), d ); // simplify } // Rationalnumber::simplify // constructors void ?{}( Rational(RationalImpl) & r ) { r{ (RationalImpl){0}, (RationalImpl){1} }; } // rational void ?{}( Rational(RationalImpl) & r, RationalImpl n ) { r{ n, (RationalImpl){1} }; } // rational void ?{}( Rational(RationalImpl) & r, RationalImpl n, RationalImpl d ) { RationalImpl t = simplify( n, d ); // simplify r.[numerator, denominator] = [n / t, d / t]; } // rational void ?{}( Rational(RationalImpl) & r, zero_t ) { r{ (RationalImpl){0}, (RationalImpl){1} }; } // rational void ?{}( Rational(RationalImpl) & r, one_t ) { r{ (RationalImpl){1}, (RationalImpl){1} }; } // rational // getter for numerator/denominator RationalImpl numerator( Rational(RationalImpl) r ) { return r.numerator; } // numerator RationalImpl denominator( Rational(RationalImpl) r ) { return r.denominator; } // denominator [ RationalImpl, RationalImpl ] ?=?( & [ RationalImpl, RationalImpl ] dest, Rational(RationalImpl) src ) { return dest = src.[ numerator, denominator ]; } // ?=? // setter for numerator/denominator RationalImpl numerator( Rational(RationalImpl) r, RationalImpl n ) { RationalImpl prev = r.numerator; RationalImpl t = gcd( abs( n ), r.denominator ); // simplify r.[numerator, denominator] = [n / t, r.denominator / t]; return prev; } // numerator RationalImpl denominator( Rational(RationalImpl) r, RationalImpl d ) { RationalImpl prev = r.denominator; RationalImpl t = simplify( r.numerator, d ); // simplify r.[numerator, denominator] = [r.numerator / t, d / t]; return prev; } // denominator // comparison int ?==?( Rational(RationalImpl) l, Rational(RationalImpl) r ) { return l.numerator * r.denominator == l.denominator * r.numerator; } // ?==? int ?!=?( Rational(RationalImpl) l, Rational(RationalImpl) r ) { return ! ( l == r ); } // ?!=? int ??( Rational(RationalImpl) l, Rational(RationalImpl) r ) { return ! ( l <= r ); } // ?>? int ?>=?( Rational(RationalImpl) l, Rational(RationalImpl) r ) { return ! ( l < r ); } // ?>=? // arithmetic Rational(RationalImpl) +?( Rational(RationalImpl) r ) { return (Rational(RationalImpl)){ r.numerator, r.denominator }; } // +? Rational(RationalImpl) -?( Rational(RationalImpl) r ) { return (Rational(RationalImpl)){ -r.numerator, r.denominator }; } // -? Rational(RationalImpl) ?+?( Rational(RationalImpl) l, Rational(RationalImpl) r ) { if ( l.denominator == r.denominator ) { // special case return (Rational(RationalImpl)){ l.numerator + r.numerator, l.denominator }; } else { return (Rational(RationalImpl)){ l.numerator * r.denominator + l.denominator * r.numerator, l.denominator * r.denominator }; } // if } // ?+? Rational(RationalImpl) ?-?( Rational(RationalImpl) l, Rational(RationalImpl) r ) { if ( l.denominator == r.denominator ) { // special case return (Rational(RationalImpl)){ l.numerator - r.numerator, l.denominator }; } else { return (Rational(RationalImpl)){ l.numerator * r.denominator - l.denominator * r.numerator, l.denominator * r.denominator }; } // if } // ?-? Rational(RationalImpl) ?*?( Rational(RationalImpl) l, Rational(RationalImpl) r ) { return (Rational(RationalImpl)){ l.numerator * r.numerator, l.denominator * r.denominator }; } // ?*? Rational(RationalImpl) ?/?( Rational(RationalImpl) l, Rational(RationalImpl) r ) { if ( r.numerator < (RationalImpl){0} ) { r.[numerator, denominator] = [-r.numerator, -r.denominator]; } // if return (Rational(RationalImpl)){ l.numerator * r.denominator, l.denominator * r.numerator }; } // ?/? // I/O forall( dtype istype | istream( istype ) | { istype & ?|?( istype &, RationalImpl & ); } ) istype & ?|?( istype & is, Rational(RationalImpl) & r ) { is | r.numerator | r.denominator; RationalImpl t = simplify( r.numerator, r.denominator ); r.numerator /= t; r.denominator /= t; return is; } // ?|? forall( dtype ostype | ostream( ostype ) | { ostype & ?|?( ostype &, RationalImpl ); } ) { ostype & ?|?( ostype & os, Rational(RationalImpl) r ) { return os | r.numerator | '/' | r.denominator; } // ?|? void ?|?( ostype & os, Rational(RationalImpl) r ) { (ostype &)(os | r); ends( os ); } // ?|? } // distribution } // distribution forall( otype RationalImpl | arithmetic( RationalImpl ) | { RationalImpl ?\?( RationalImpl, unsigned long ); } ) Rational(RationalImpl) ?\?( Rational(RationalImpl) x, long int y ) { if ( y < 0 ) { return (Rational(RationalImpl)){ x.denominator \ -y, x.numerator \ -y }; } else { return (Rational(RationalImpl)){ x.numerator \ y, x.denominator \ y }; } // if } // conversion forall( otype RationalImpl | arithmetic( RationalImpl ) | { double convert( RationalImpl ); } ) double widen( Rational(RationalImpl) r ) { return convert( r.numerator ) / convert( r.denominator ); } // widen forall( otype RationalImpl | arithmetic( RationalImpl ) | { double convert( RationalImpl ); RationalImpl convert( double ); } ) Rational(RationalImpl) narrow( double f, RationalImpl md ) { // http://www.ics.uci.edu/~eppstein/numth/frap.c if ( md <= (RationalImpl){1} ) { // maximum fractional digits too small? return (Rational(RationalImpl)){ convert( f ), (RationalImpl){1}}; // truncate fraction } // if // continued fraction coefficients RationalImpl m00 = {1}, m11 = { 1 }, m01 = { 0 }, m10 = { 0 }; RationalImpl ai, t; // find terms until denom gets too big for ( ;; ) { ai = convert( f ); if ( ! (m10 * ai + m11 <= md) ) break; t = m00 * ai + m01; m01 = m00; m00 = t; t = m10 * ai + m11; m11 = m10; m10 = t; double temp = convert( ai ); if ( f == temp ) break; // prevent division by zero f = 1 / (f - temp); if ( f > (double)0x7FFFFFFF ) break; // representation failure } // for return (Rational(RationalImpl)){ m00, m10 }; } // narrow // Local Variables: // // tab-width: 4 // // End: //