Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/libcfa/rational.c

    r3ce0d440 r09687aa  
    1010// Created On       : Wed Apr  6 17:54:28 2016
    1111// Last Modified By : Peter A. Buhr
    12 // Last Modified On : Sat Jun  2 09:24:33 2018
    13 // Update Count     : 162
     12// Last Modified On : Wed Dec  6 23:13:58 2017
     13// Update Count     : 156
    1414//
    1515
    … …  
    1818#include "stdlib"
    1919
    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 };
    123                 return t;
    124         } // +?
    125 
    126         Rational(RationalImpl) -?( Rational(RationalImpl) r ) {
    127                 Rational(RationalImpl) t = { -r.numerator, r.denominator };
    128                 return t;
    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 };
    153                 return t;
    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
     20// helper routines
     21
     22// Calculate greatest common denominator of two numbers, the first of which may be negative. Used to reduce rationals.
     23// alternative: https://en.wikipedia.org/wiki/Binary_GCD_algorithm
     24forall( otype RationalImpl | arithmetic( RationalImpl ) )
     25static 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
     35forall( otype RationalImpl | arithmetic( RationalImpl ) )
     36static RationalImpl simplify( RationalImpl & n, RationalImpl & d ) {
     37        if ( d == (RationalImpl){0} ) {
     38                serr | "Invalid rational number construction: denominator cannot be equal to 0." | endl;
     39                exit( EXIT_FAILURE );
     40        } // exit
     41        if ( d < (RationalImpl){0} ) { d = -d; n = -n; }        // move sign to numerator
     42        return gcd( abs( n ), d );                                                      // simplify
     43} // Rationalnumber::simplify
     44
     45
     46// constructors
     47
     48forall( otype RationalImpl | arithmetic( RationalImpl ) )
     49void ?{}( Rational(RationalImpl) & r ) {
     50        r{ (RationalImpl){0}, (RationalImpl){1} };
     51} // rational
     52
     53forall( otype RationalImpl | arithmetic( RationalImpl ) )
     54void ?{}( Rational(RationalImpl) & r, RationalImpl n ) {
     55        r{ n, (RationalImpl){1} };
     56} // rational
     57
     58forall( otype RationalImpl | arithmetic( RationalImpl ) )
     59void ?{}( Rational(RationalImpl) & r, RationalImpl n, RationalImpl d ) {
     60        RationalImpl t = simplify( n, d );                                      // simplify
     61        r.numerator = n / t;
     62        r.denominator = d / t;
     63} // rational
     64
     65
     66// getter for numerator/denominator
     67
     68forall( otype RationalImpl | arithmetic( RationalImpl ) )
     69RationalImpl numerator( Rational(RationalImpl) r ) {
     70        return r.numerator;
     71} // numerator
     72
     73forall( otype RationalImpl | arithmetic( RationalImpl ) )
     74RationalImpl denominator( Rational(RationalImpl) r ) {
     75        return r.denominator;
     76} // denominator
     77
     78forall( otype RationalImpl | arithmetic( RationalImpl ) )
     79[ RationalImpl, RationalImpl ] ?=?( & [ RationalImpl, RationalImpl ] dest, Rational(RationalImpl) src ) {
     80        return dest = src.[ numerator, denominator ];
     81}
     82
     83// setter for numerator/denominator
     84
     85forall( otype RationalImpl | arithmetic( RationalImpl ) )
     86RationalImpl numerator( Rational(RationalImpl) r, RationalImpl n ) {
     87        RationalImpl prev = r.numerator;
     88        RationalImpl t = gcd( abs( n ), r.denominator );                // simplify
     89        r.numerator = n / t;
     90        r.denominator = r.denominator / t;
     91        return prev;
     92} // numerator
     93
     94forall( otype RationalImpl | arithmetic( RationalImpl ) )
     95RationalImpl denominator( Rational(RationalImpl) r, RationalImpl d ) {
     96        RationalImpl prev = r.denominator;
     97        RationalImpl t = simplify( r.numerator, d );                    // simplify
     98        r.numerator = r.numerator / t;
     99        r.denominator = d / t;
     100        return prev;
     101} // denominator
     102
     103
     104// comparison
     105
     106forall( otype RationalImpl | arithmetic( RationalImpl ) )
     107int ?==?( Rational(RationalImpl) l, Rational(RationalImpl) r ) {
     108        return l.numerator * r.denominator == l.denominator * r.numerator;
     109} // ?==?
     110
     111forall( otype RationalImpl | arithmetic( RationalImpl ) )
     112int ?!=?( Rational(RationalImpl) l, Rational(RationalImpl) r ) {
     113        return ! ( l == r );
     114} // ?!=?
     115
     116forall( otype RationalImpl | arithmetic( RationalImpl ) )
     117int ?<?( Rational(RationalImpl) l, Rational(RationalImpl) r ) {
     118        return l.numerator * r.denominator < l.denominator * r.numerator;
     119} // ?<?
     120
     121forall( otype RationalImpl | arithmetic( RationalImpl ) )
     122int ?<=?( Rational(RationalImpl) l, Rational(RationalImpl) r ) {
     123        return l.numerator * r.denominator <= l.denominator * r.numerator;
     124} // ?<=?
     125
     126forall( otype RationalImpl | arithmetic( RationalImpl ) )
     127int ?>?( Rational(RationalImpl) l, Rational(RationalImpl) r ) {
     128        return ! ( l <= r );
     129} // ?>?
     130
     131forall( otype RationalImpl | arithmetic( RationalImpl ) )
     132int ?>=?( Rational(RationalImpl) l, Rational(RationalImpl) r ) {
     133        return ! ( l < r );
     134} // ?>=?
     135
     136
     137// arithmetic
     138
     139forall( otype RationalImpl | arithmetic( RationalImpl ) )
     140Rational(RationalImpl) +?( Rational(RationalImpl) r ) {
     141        Rational(RationalImpl) t = { r.numerator, r.denominator };
     142        return t;
     143} // +?
     144
     145forall( otype RationalImpl | arithmetic( RationalImpl ) )
     146Rational(RationalImpl) -?( Rational(RationalImpl) r ) {
     147        Rational(RationalImpl) t = { -r.numerator, r.denominator };
     148        return t;
     149} // -?
     150
     151forall( otype RationalImpl | arithmetic( RationalImpl ) )
     152Rational(RationalImpl) ?+?( Rational(RationalImpl) l, Rational(RationalImpl) r ) {
     153        if ( l.denominator == r.denominator ) {                         // special case
     154                Rational(RationalImpl) t = { l.numerator + r.numerator, l.denominator };
     155                return t;
     156        } else {
     157                Rational(RationalImpl) t = { l.numerator * r.denominator + l.denominator * r.numerator, l.denominator * r.denominator };
     158                return t;
     159        } // if
     160} // ?+?
     161
     162forall( otype RationalImpl | arithmetic( RationalImpl ) )
     163Rational(RationalImpl) ?-?( Rational(RationalImpl) l, Rational(RationalImpl) r ) {
     164        if ( l.denominator == r.denominator ) {                         // special case
     165                Rational(RationalImpl) t = { l.numerator - r.numerator, l.denominator };
     166                return t;
     167        } else {
     168                Rational(RationalImpl) t = { l.numerator * r.denominator - l.denominator * r.numerator, l.denominator * r.denominator };
     169                return t;
     170        } // if
     171} // ?-?
     172
     173forall( otype RationalImpl | arithmetic( RationalImpl ) )
     174Rational(RationalImpl) ?*?( Rational(RationalImpl) l, Rational(RationalImpl) r ) {
     175        Rational(RationalImpl) t = { l.numerator * r.numerator, l.denominator * r.denominator };
     176        return t;
     177} // ?*?
     178
     179forall( otype RationalImpl | arithmetic( RationalImpl ) )
     180Rational(RationalImpl) ?/?( Rational(RationalImpl) l, Rational(RationalImpl) r ) {
     181        if ( r.numerator < (RationalImpl){0} ) {
     182                r.numerator = -r.numerator;
     183                r.denominator = -r.denominator;
     184        } // if
     185        Rational(RationalImpl) t = { l.numerator * r.denominator, l.denominator * r.numerator };
     186        return t;
     187} // ?/?
     188
    182189
    183190// conversion
    … …  
    188195} // widen
    189196
     197// http://www.ics.uci.edu/~eppstein/numth/frap.c
    190198forall( otype RationalImpl | arithmetic( RationalImpl ) | { double convert( RationalImpl ); RationalImpl convert( double ); } )
    191199Rational(RationalImpl) narrow( double f, RationalImpl md ) {
    192         // http://www.ics.uci.edu/~eppstein/numth/frap.c
    193200        if ( md <= (RationalImpl){1} ) {                                        // maximum fractional digits too small?
    194201                return (Rational(RationalImpl)){ convert( f ), (RationalImpl){1}}; // truncate fraction
    … …  
    217224} // narrow
    218225
     226
     227// I/O
     228
     229forall( otype RationalImpl | arithmetic( RationalImpl ) )
     230forall( dtype istype | istream( istype ) | { istype & ?|?( istype &, RationalImpl & ); } )
     231istype & ?|?( istype & is, Rational(RationalImpl) & r ) {
     232        RationalImpl t;
     233        is | r.numerator | r.denominator;
     234        t = simplify( r.numerator, r.denominator );
     235        r.numerator /= t;
     236        r.denominator /= t;
     237        return is;
     238} // ?|?
     239
     240forall( otype RationalImpl | arithmetic( RationalImpl ) )
     241forall( dtype ostype | ostream( ostype ) | { ostype & ?|?( ostype &, RationalImpl ); } )
     242ostype & ?|?( ostype & os, Rational(RationalImpl ) r ) {
     243        return os | r.numerator | '/' | r.denominator;
     244} // ?|?
     245
    219246// Local Variables: //
    220247// tab-width: 4 //
Note: See TracChangeset for help on using the changeset viewer.