Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed
  • src/libcfa/rational.c

    r561f730 rf621a148  
    1010// Created On       : Wed Apr  6 17:54:28 2016
    1111// Last Modified By : Peter A. Buhr
    12 // Last Modified On : Sun May 14 17:25:19 2017
    13 // Update Count     : 131
     12// Last Modified On : Thu Apr 27 17:05:06 2017
     13// Update Count     : 51
    1414//
    1515
     
    1717#include "fstream"
    1818#include "stdlib"
     19#include "math"                                                                                 // floor
     20
     21
     22// constants
     23
     24struct Rational 0 = {0, 1};
     25struct Rational 1 = {1, 1};
     26
    1927
    2028// helper routines
     
    2230// Calculate greatest common denominator of two numbers, the first of which may be negative. Used to reduce rationals.
    2331// alternative: https://en.wikipedia.org/wiki/Binary_GCD_algorithm
    24 forall ( otype RationalImpl | arithmetic( RationalImpl ) )
    2532static RationalImpl gcd( RationalImpl a, RationalImpl b ) {
    2633        for ( ;; ) {                                                                            // Euclid's algorithm
    2734                RationalImpl r = a % b;
    28           if ( r == (RationalImpl){0} ) break;
     35          if ( r == 0 ) break;
    2936                a = b;
    3037                b = r;
     
    3340} // gcd
    3441
    35 forall ( otype RationalImpl | arithmetic( RationalImpl ) )
    36 static RationalImpl simplify( RationalImpl * n, RationalImpl * d ) {
    37         if ( *d == (RationalImpl){0} ) {
     42static RationalImpl simplify( RationalImpl *n, RationalImpl *d ) {
     43        if ( *d == 0 ) {
    3844                serr | "Invalid rational number construction: denominator cannot be equal to 0." | endl;
    3945                exit( EXIT_FAILURE );
    4046        } // exit
    41         if ( *d < (RationalImpl){0} ) { *d = -*d; *n = -*n; } // move sign to numerator
     47        if ( *d < 0 ) { *d = -*d; *n = -*n; }                           // move sign to numerator
    4248        return gcd( abs( *n ), *d );                                            // simplify
    4349} // Rationalnumber::simplify
     
    4652// constructors
    4753
    48 forall ( otype RationalImpl | arithmetic( RationalImpl ) )
    49 void ?{}( Rational(RationalImpl) * r ) {
    50         r{ (RationalImpl){0}, (RationalImpl){1} };
     54void ?{}( Rational * r ) {
     55        r{ 0, 1 };
    5156} // rational
    5257
    53 forall ( otype RationalImpl | arithmetic( RationalImpl ) )
    54 void ?{}( Rational(RationalImpl) * r, RationalImpl n ) {
    55         r{ n, (RationalImpl){1} };
     58void ?{}( Rational * r, RationalImpl n ) {
     59        r{ n, 1 };
    5660} // rational
    5761
    58 forall ( otype RationalImpl | arithmetic( RationalImpl ) )
    59 void ?{}( Rational(RationalImpl) * r, RationalImpl n, RationalImpl d ) {
     62void ?{}( Rational * r, RationalImpl n, RationalImpl d ) {
    6063        RationalImpl t = simplify( &n, &d );                            // simplify
    6164        r->numerator = n / t;
     
    6669// getter for numerator/denominator
    6770
    68 forall ( otype RationalImpl | arithmetic( RationalImpl ) )
    69 RationalImpl numerator( Rational(RationalImpl) r ) {
     71RationalImpl numerator( Rational r ) {
    7072        return r.numerator;
    7173} // numerator
    7274
    73 forall ( otype RationalImpl | arithmetic( RationalImpl ) )
    74 RationalImpl denominator( Rational(RationalImpl) r ) {
     75RationalImpl denominator( Rational r ) {
    7576        return r.denominator;
    7677} // denominator
    7778
    78 forall ( otype RationalImpl | arithmetic( RationalImpl ) )
    79 [ RationalImpl, RationalImpl ] ?=?( * [ RationalImpl, RationalImpl ] dest, Rational(RationalImpl) src ) {
     79[ RationalImpl, RationalImpl ] ?=?( * [ RationalImpl, RationalImpl ] dest, Rational src ) {
    8080        return *dest = src.[ numerator, denominator ];
    8181}
     
    8383// setter for numerator/denominator
    8484
    85 forall ( otype RationalImpl | arithmetic( RationalImpl ) )
    86 RationalImpl numerator( Rational(RationalImpl) r, RationalImpl n ) {
     85RationalImpl numerator( Rational r, RationalImpl n ) {
    8786        RationalImpl prev = r.numerator;
    8887        RationalImpl t = gcd( abs( n ), r.denominator );                // simplify
     
    9291} // numerator
    9392
    94 forall ( otype RationalImpl | arithmetic( RationalImpl ) )
    95 RationalImpl denominator( Rational(RationalImpl) r, RationalImpl d ) {
     93RationalImpl denominator( Rational r, RationalImpl d ) {
    9694        RationalImpl prev = r.denominator;
    9795        RationalImpl t = simplify( &r.numerator, &d );                  // simplify
     
    104102// comparison
    105103
    106 forall ( otype RationalImpl | arithmetic( RationalImpl ) )
    107 int ?==?( Rational(RationalImpl) l, Rational(RationalImpl) r ) {
     104int ?==?( Rational l, Rational r ) {
    108105        return l.numerator * r.denominator == l.denominator * r.numerator;
    109106} // ?==?
    110107
    111 forall ( otype RationalImpl | arithmetic( RationalImpl ) )
    112 int ?!=?( Rational(RationalImpl) l, Rational(RationalImpl) r ) {
     108int ?!=?( Rational l, Rational r ) {
    113109        return ! ( l == r );
    114110} // ?!=?
    115111
    116 forall ( otype RationalImpl | arithmetic( RationalImpl ) )
    117 int ?<?( Rational(RationalImpl) l, Rational(RationalImpl) r ) {
     112int ?<?( Rational l, Rational r ) {
    118113        return l.numerator * r.denominator < l.denominator * r.numerator;
    119114} // ?<?
    120115
    121 forall ( otype RationalImpl | arithmetic( RationalImpl ) )
    122 int ?<=?( Rational(RationalImpl) l, Rational(RationalImpl) r ) {
    123         return l.numerator * r.denominator <= l.denominator * r.numerator;
     116int ?<=?( Rational l, Rational r ) {
     117        return l < r || l == r;
    124118} // ?<=?
    125119
    126 forall ( otype RationalImpl | arithmetic( RationalImpl ) )
    127 int ?>?( Rational(RationalImpl) l, Rational(RationalImpl) r ) {
     120int ?>?( Rational l, Rational r ) {
    128121        return ! ( l <= r );
    129122} // ?>?
    130123
    131 forall ( otype RationalImpl | arithmetic( RationalImpl ) )
    132 int ?>=?( Rational(RationalImpl) l, Rational(RationalImpl) r ) {
     124int ?>=?( Rational l, Rational r ) {
    133125        return ! ( l < r );
    134126} // ?>=?
     
    137129// arithmetic
    138130
    139 forall ( otype RationalImpl | arithmetic( RationalImpl ) )
    140 Rational(RationalImpl) +?( Rational(RationalImpl) r ) {
    141         Rational(RationalImpl) t = { r.numerator, r.denominator };
    142         return t;
    143 } // +?
    144 
    145 forall ( otype RationalImpl | arithmetic( RationalImpl ) )
    146 Rational(RationalImpl) -?( Rational(RationalImpl) r ) {
    147         Rational(RationalImpl) t = { -r.numerator, r.denominator };
     131Rational -?( Rational r ) {
     132        Rational t = { -r.numerator, r.denominator };
    148133        return t;
    149134} // -?
    150135
    151 forall ( otype RationalImpl | arithmetic( RationalImpl ) )
    152 Rational(RationalImpl) ?+?( Rational(RationalImpl) l, Rational(RationalImpl) r ) {
     136Rational ?+?( Rational l, Rational r ) {
    153137        if ( l.denominator == r.denominator ) {                         // special case
    154                 Rational(RationalImpl) t = { l.numerator + r.numerator, l.denominator };
     138                Rational t = { l.numerator + r.numerator, l.denominator };
    155139                return t;
    156140        } else {
    157                 Rational(RationalImpl) t = { l.numerator * r.denominator + l.denominator * r.numerator, l.denominator * r.denominator };
     141                Rational t = { l.numerator * r.denominator + l.denominator * r.numerator, l.denominator * r.denominator };
    158142                return t;
    159143        } // if
    160144} // ?+?
    161145
    162 forall ( otype RationalImpl | arithmetic( RationalImpl ) )
    163 Rational(RationalImpl) ?-?( Rational(RationalImpl) l, Rational(RationalImpl) r ) {
     146Rational ?-?( Rational l, Rational r ) {
    164147        if ( l.denominator == r.denominator ) {                         // special case
    165                 Rational(RationalImpl) t = { l.numerator - r.numerator, l.denominator };
     148                Rational t = { l.numerator - r.numerator, l.denominator };
    166149                return t;
    167150        } else {
    168                 Rational(RationalImpl) t = { l.numerator * r.denominator - l.denominator * r.numerator, l.denominator * r.denominator };
     151                Rational t = { l.numerator * r.denominator - l.denominator * r.numerator, l.denominator * r.denominator };
    169152                return t;
    170153        } // if
    171154} // ?-?
    172155
    173 forall ( otype RationalImpl | arithmetic( RationalImpl ) )
    174 Rational(RationalImpl) ?*?( Rational(RationalImpl) l, Rational(RationalImpl) r ) {
    175         Rational(RationalImpl) t = { l.numerator * r.numerator, l.denominator * r.denominator };
     156Rational ?*?( Rational l, Rational r ) {
     157        Rational t = { l.numerator * r.numerator, l.denominator * r.denominator };
    176158        return t;
    177159} // ?*?
    178160
    179 forall ( otype RationalImpl | arithmetic( RationalImpl ) )
    180 Rational(RationalImpl) ?/?( Rational(RationalImpl) l, Rational(RationalImpl) r ) {
    181         if ( r.numerator < (RationalImpl){0} ) {
     161Rational ?/?( Rational l, Rational r ) {
     162        if ( r.numerator < 0 ) {
    182163                r.numerator = -r.numerator;
    183164                r.denominator = -r.denominator;
    184165        } // if
    185         Rational(RationalImpl) t = { l.numerator * r.denominator, l.denominator * r.numerator };
     166        Rational t = { l.numerator * r.denominator, l.denominator * r.numerator };
    186167        return t;
    187168} // ?/?
     
    190171// conversion
    191172
    192 // forall ( otype RationalImpl | arithmetic( RationalImpl ) )
    193 // double widen( Rational(RationalImpl) r ) {
    194 //      return (double)r.numerator / (double)r.denominator;
    195 // } // widen
    196 
    197 // // http://www.ics.uci.edu/~eppstein/numth/frap.c
    198 // forall ( otype RationalImpl | arithmetic( RationalImpl ) )
    199 // Rational(RationalImpl) narrow( double f, RationalImpl md ) {
    200 //      if ( md <= 1 ) {                                                                        // maximum fractional digits too small?
    201 //              return (Rational(RationalImpl)){ f, 1};                 // truncate fraction
    202 //      } // if
    203 
    204 //      // continued fraction coefficients
    205 //      RationalImpl m00 = 1, m11 = 1, m01 = 0, m10 = 0;
    206 //      RationalImpl ai, t;
    207 
    208 //      // find terms until denom gets too big
    209 //      for ( ;; ) {
    210 //              ai = (RationalImpl)f;
    211 //        if ( ! (m10 * ai + m11 <= md) ) break;
    212 //              t = m00 * ai + m01;
    213 //              m01 = m00;
    214 //              m00 = t;
    215 //              t = m10 * ai + m11;
    216 //              m11 = m10;
    217 //              m10 = t;
    218 //              t = (double)ai;
    219 //        if ( f == t ) break;                                                          // prevent division by zero
    220 //        f = 1 / (f - (double)t);
    221 //        if ( f > (double)0x7FFFFFFF ) break;                          // representation failure
    222 //      }
    223 //      return (Rational(RationalImpl)){ m00, m10 };
    224 // } // narrow
     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
    225204
    226205
    227206// I/O
    228207
    229 forall ( otype RationalImpl | arithmetic( RationalImpl ) )
    230 forall( dtype istype | istream( istype ) | { istype * ?|?( istype *, RationalImpl * ); } )
    231 istype * ?|?( istype * is, Rational(RationalImpl) * r ) {
     208forall( dtype istype | istream( istype ) )
     209istype * ?|?( istype *is, Rational *r ) {
    232210        RationalImpl t;
    233211        is | &(r->numerator) | &(r->denominator);
     
    238216} // ?|?
    239217
    240 forall ( otype RationalImpl | arithmetic( RationalImpl ) )
    241 forall( dtype ostype | ostream( ostype ) | { ostype * ?|?( ostype *, RationalImpl ); } )
    242 ostype * ?|?( ostype * os, Rational(RationalImpl ) r ) {
     218forall( dtype ostype | ostream( ostype ) )
     219ostype * ?|?( ostype *os, Rational r ) {
    243220        return os | r.numerator | '/' | r.denominator;
    244221} // ?|?
Note: See TracChangeset for help on using the changeset viewer.