Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • libcfa/src/rational.cfa

    r5dc4c7e rfd54fef  
    1010// Created On       : Wed Apr  6 17:54:28 2016
    1111// Last Modified By : Peter A. Buhr
    12 // Last Modified On : Tue Jul 20 16:30:06 2021
    13 // Update Count     : 193
     12// Last Modified On : Sat Feb  8 17:56:36 2020
     13// Update Count     : 187
    1414//
    1515
    … …  
    1818#include "stdlib.hfa"
    1919
    20 forall( T | Arithmetic( T ) ) {
     20forall( RationalImpl | arithmetic( RationalImpl ) ) {
    2121        // helper routines
    2222
    2323        // Calculate greatest common denominator of two numbers, the first of which may be negative. Used to reduce
    2424        // rationals.  alternative: https://en.wikipedia.org/wiki/Binary_GCD_algorithm
    25         static T gcd( T a, T b ) {
     25        static RationalImpl gcd( RationalImpl a, RationalImpl b ) {
    2626                for ( ;; ) {                                                                    // Euclid's algorithm
    27                         T r = a % b;
    28                   if ( r == (T){0} ) break;
     27                        RationalImpl r = a % b;
     28                  if ( r == (RationalImpl){0} ) break;
    2929                        a = b;
    3030                        b = r;
    … …  
    3333        } // gcd
    3434
    35         static T simplify( T & n, T & d ) {
    36                 if ( d == (T){0} ) {
     35        static RationalImpl simplify( RationalImpl & n, RationalImpl & d ) {
     36                if ( d == (RationalImpl){0} ) {
    3737                        abort | "Invalid rational number construction: denominator cannot be equal to 0.";
    3838                } // exit
    39                 if ( d < (T){0} ) { d = -d; n = -n; } // move sign to numerator
     39                if ( d < (RationalImpl){0} ) { d = -d; n = -n; } // move sign to numerator
    4040                return gcd( abs( n ), d );                                              // simplify
    4141        } // Rationalnumber::simplify
    … …  
    4343        // constructors
    4444
    45         void ?{}( Rational(T) & r, zero_t ) {
    46                 r{ (T){0}, (T){1} };
    47         } // rational
    48 
    49         void ?{}( Rational(T) & r, one_t ) {
    50                 r{ (T){1}, (T){1} };
    51         } // rational
    52 
    53         void ?{}( Rational(T) & r ) {
    54                 r{ (T){0}, (T){1} };
    55         } // rational
    56 
    57         void ?{}( Rational(T) & r, T n ) {
    58                 r{ n, (T){1} };
    59         } // rational
    60 
    61         void ?{}( Rational(T) & r, T n, T d ) {
    62                 T t = simplify( n, d );                         // simplify
     45        void ?{}( Rational(RationalImpl) & r ) {
     46                r{ (RationalImpl){0}, (RationalImpl){1} };
     47        } // rational
     48
     49        void ?{}( Rational(RationalImpl) & r, RationalImpl n ) {
     50                r{ n, (RationalImpl){1} };
     51        } // rational
     52
     53        void ?{}( Rational(RationalImpl) & r, RationalImpl n, RationalImpl d ) {
     54                RationalImpl t = simplify( n, d );                              // simplify
    6355                r.[numerator, denominator] = [n / t, d / t];
    6456        } // rational
    6557
     58        void ?{}( Rational(RationalImpl) & r, zero_t ) {
     59                r{ (RationalImpl){0}, (RationalImpl){1} };
     60        } // rational
     61
     62        void ?{}( Rational(RationalImpl) & r, one_t ) {
     63                r{ (RationalImpl){1}, (RationalImpl){1} };
     64        } // rational
     65
    6666        // getter for numerator/denominator
    6767
    68         T numerator( Rational(T) r ) {
     68        RationalImpl numerator( Rational(RationalImpl) r ) {
    6969                return r.numerator;
    7070        } // numerator
    7171
    72         T denominator( Rational(T) r ) {
     72        RationalImpl denominator( Rational(RationalImpl) r ) {
    7373                return r.denominator;
    7474        } // denominator
    7575
    76         [ T, T ] ?=?( & [ T, T ] dest, Rational(T) src ) {
     76        [ RationalImpl, RationalImpl ] ?=?( & [ RationalImpl, RationalImpl ] dest, Rational(RationalImpl) src ) {
    7777                return dest = src.[ numerator, denominator ];
    7878        } // ?=?
    … …  
    8080        // setter for numerator/denominator
    8181
    82         T numerator( Rational(T) r, T n ) {
    83                 T prev = r.numerator;
    84                 T t = gcd( abs( n ), r.denominator ); // simplify
     82        RationalImpl numerator( Rational(RationalImpl) r, RationalImpl n ) {
     83                RationalImpl prev = r.numerator;
     84                RationalImpl t = gcd( abs( n ), r.denominator ); // simplify
    8585                r.[numerator, denominator] = [n / t, r.denominator / t];
    8686                return prev;
    8787        } // numerator
    8888
    89         T denominator( Rational(T) r, T d ) {
    90                 T prev = r.denominator;
    91                 T t = simplify( r.numerator, d );       // simplify
     89        RationalImpl denominator( Rational(RationalImpl) r, RationalImpl d ) {
     90                RationalImpl prev = r.denominator;
     91                RationalImpl t = simplify( r.numerator, d );    // simplify
    9292                r.[numerator, denominator] = [r.numerator / t, d / t];
    9393                return prev;
    … …  
    9696        // comparison
    9797
    98         int ?==?( Rational(T) l, Rational(T) r ) {
     98        int ?==?( Rational(RationalImpl) l, Rational(RationalImpl) r ) {
    9999                return l.numerator * r.denominator == l.denominator * r.numerator;
    100100        } // ?==?
    101101
    102         int ?!=?( Rational(T) l, Rational(T) r ) {
     102        int ?!=?( Rational(RationalImpl) l, Rational(RationalImpl) r ) {
    103103                return ! ( l == r );
    104104        } // ?!=?
    105105
    106         int ?!=?( Rational(T) l, zero_t ) {
    107                 return ! ( l == (Rational(T)){ 0 } );
    108         } // ?!=?
    109 
    110         int ?<?( Rational(T) l, Rational(T) r ) {
     106        int ?<?( Rational(RationalImpl) l, Rational(RationalImpl) r ) {
    111107                return l.numerator * r.denominator < l.denominator * r.numerator;
    112108        } // ?<?
    113109
    114         int ?<=?( Rational(T) l, Rational(T) r ) {
     110        int ?<=?( Rational(RationalImpl) l, Rational(RationalImpl) r ) {
    115111                return l.numerator * r.denominator <= l.denominator * r.numerator;
    116112        } // ?<=?
    117113
    118         int ?>?( Rational(T) l, Rational(T) r ) {
     114        int ?>?( Rational(RationalImpl) l, Rational(RationalImpl) r ) {
    119115                return ! ( l <= r );
    120116        } // ?>?
    121117
    122         int ?>=?( Rational(T) l, Rational(T) r ) {
     118        int ?>=?( Rational(RationalImpl) l, Rational(RationalImpl) r ) {
    123119                return ! ( l < r );
    124120        } // ?>=?
    … …  
    126122        // arithmetic
    127123
    128         Rational(T) +?( Rational(T) r ) {
    129                 return (Rational(T)){ r.numerator, r.denominator };
     124        Rational(RationalImpl) +?( Rational(RationalImpl) r ) {
     125                return (Rational(RationalImpl)){ r.numerator, r.denominator };
    130126        } // +?
    131127
    132         Rational(T) -?( Rational(T) r ) {
    133                 return (Rational(T)){ -r.numerator, r.denominator };
     128        Rational(RationalImpl) -?( Rational(RationalImpl) r ) {
     129                return (Rational(RationalImpl)){ -r.numerator, r.denominator };
    134130        } // -?
    135131
    136         Rational(T) ?+?( Rational(T) l, Rational(T) r ) {
     132        Rational(RationalImpl) ?+?( Rational(RationalImpl) l, Rational(RationalImpl) r ) {
    137133                if ( l.denominator == r.denominator ) {                 // special case
    138                         return (Rational(T)){ l.numerator + r.numerator, l.denominator };
     134                        return (Rational(RationalImpl)){ l.numerator + r.numerator, l.denominator };
    139135                } else {
    140                         return (Rational(T)){ l.numerator * r.denominator + l.denominator * r.numerator, l.denominator * r.denominator };
     136                        return (Rational(RationalImpl)){ l.numerator * r.denominator + l.denominator * r.numerator, l.denominator * r.denominator };
    141137                } // if
    142138        } // ?+?
    143139
    144         Rational(T) ?+=?( Rational(T) & l, Rational(T) r ) {
    145                 l = l + r;
    146                 return l;
    147         } // ?+?
    148 
    149         Rational(T) ?+=?( Rational(T) & l, one_t ) {
    150                 l = l + (Rational(T)){ 1 };
    151                 return l;
    152         } // ?+?
    153 
    154         Rational(T) ?-?( Rational(T) l, Rational(T) r ) {
     140        Rational(RationalImpl) ?-?( Rational(RationalImpl) l, Rational(RationalImpl) r ) {
    155141                if ( l.denominator == r.denominator ) {                 // special case
    156                         return (Rational(T)){ l.numerator - r.numerator, l.denominator };
     142                        return (Rational(RationalImpl)){ l.numerator - r.numerator, l.denominator };
    157143                } else {
    158                         return (Rational(T)){ l.numerator * r.denominator - l.denominator * r.numerator, l.denominator * r.denominator };
     144                        return (Rational(RationalImpl)){ l.numerator * r.denominator - l.denominator * r.numerator, l.denominator * r.denominator };
    159145                } // if
    160146        } // ?-?
    161147
    162         Rational(T) ?-=?( Rational(T) & l, Rational(T) r ) {
    163                 l = l - r;
    164                 return l;
    165         } // ?-?
    166 
    167         Rational(T) ?-=?( Rational(T) & l, one_t ) {
    168                 l = l - (Rational(T)){ 1 };
    169                 return l;
    170         } // ?-?
    171 
    172         Rational(T) ?*?( Rational(T) l, Rational(T) r ) {
    173                 return (Rational(T)){ l.numerator * r.numerator, l.denominator * r.denominator };
     148        Rational(RationalImpl) ?*?( Rational(RationalImpl) l, Rational(RationalImpl) r ) {
     149                return (Rational(RationalImpl)){ l.numerator * r.numerator, l.denominator * r.denominator };
    174150        } // ?*?
    175151
    176         Rational(T) ?*=?( Rational(T) & l, Rational(T) r ) {
    177                 return l = l * r;
    178         } // ?*?
    179 
    180         Rational(T) ?/?( Rational(T) l, Rational(T) r ) {
    181                 if ( r.numerator < (T){0} ) {
     152        Rational(RationalImpl) ?/?( Rational(RationalImpl) l, Rational(RationalImpl) r ) {
     153                if ( r.numerator < (RationalImpl){0} ) {
    182154                        r.[numerator, denominator] = [-r.numerator, -r.denominator];
    183155                } // if
    184                 return (Rational(T)){ l.numerator * r.denominator, l.denominator * r.numerator };
     156                return (Rational(RationalImpl)){ l.numerator * r.denominator, l.denominator * r.numerator };
    185157        } // ?/?
    186158
    187         Rational(T) ?/=?( Rational(T) & l, Rational(T) r ) {
    188                 return l = l / r;
    189         } // ?/?
    190 
    191159        // I/O
    192160
    193         forall( istype & | istream( istype ) | { istype & ?|?( istype &, T & ); } )
    194         istype & ?|?( istype & is, Rational(T) & r ) {
     161        forall( istype & | istream( istype ) | { istype & ?|?( istype &, RationalImpl & ); } )
     162        istype & ?|?( istype & is, Rational(RationalImpl) & r ) {
    195163                is | r.numerator | r.denominator;
    196                 T t = simplify( r.numerator, r.denominator );
     164                RationalImpl t = simplify( r.numerator, r.denominator );
    197165                r.numerator /= t;
    198166                r.denominator /= t;
    … …  
    200168        } // ?|?
    201169
    202         forall( ostype & | ostream( ostype ) | { ostype & ?|?( ostype &, T ); } ) {
    203                 ostype & ?|?( ostype & os, Rational(T) r ) {
     170        forall( ostype & | ostream( ostype ) | { ostype & ?|?( ostype &, RationalImpl ); } ) {
     171                ostype & ?|?( ostype & os, Rational(RationalImpl) r ) {
    204172                        return os | r.numerator | '/' | r.denominator;
    205173                } // ?|?
    206174
    207                 void ?|?( ostype & os, Rational(T) r ) {
     175                void ?|?( ostype & os, Rational(RationalImpl) r ) {
    208176                        (ostype &)(os | r); ends( os );
    209177                } // ?|?
    … …  
    211179} // distribution
    212180
    213 forall( T | Arithmetic( T ) | { T ?\?( T, unsigned long ); } ) {
    214         Rational(T) ?\?( Rational(T) x, long int y ) {
    215                 if ( y < 0 ) {
    216                         return (Rational(T)){ x.denominator \ -y, x.numerator \ -y };
    217                 } else {
    218                         return (Rational(T)){ x.numerator \ y, x.denominator \ y };
    219                 } // if
    220         } // ?\?
    221 
    222         Rational(T) ?\=?( Rational(T) & x, long int y ) {
    223                 return x = x \ y;
    224         } // ?\?
    225 } // distribution
     181forall( RationalImpl | arithmetic( RationalImpl ) | { RationalImpl ?\?( RationalImpl, unsigned long ); } )
     182Rational(RationalImpl) ?\?( Rational(RationalImpl) x, long int y ) {
     183        if ( y < 0 ) {
     184                return (Rational(RationalImpl)){ x.denominator \ -y, x.numerator \ -y };
     185        } else {
     186                return (Rational(RationalImpl)){ x.numerator \ y, x.denominator \ y };
     187        } // if
     188}
    226189
    227190// conversion
    228191
    229 forall( T | Arithmetic( T ) | { double convert( T ); } )
    230 double widen( Rational(T) r ) {
     192forall( RationalImpl | arithmetic( RationalImpl ) | { double convert( RationalImpl ); } )
     193double widen( Rational(RationalImpl) r ) {
    231194        return convert( r.numerator ) / convert( r.denominator );
    232195} // widen
    233196
    234 forall( T | Arithmetic( T ) | { double convert( T ); T convert( double ); } )
    235 Rational(T) narrow( double f, T md ) {
     197forall( RationalImpl | arithmetic( RationalImpl ) | { double convert( RationalImpl ); RationalImpl convert( double ); } )
     198Rational(RationalImpl) narrow( double f, RationalImpl md ) {
    236199        // http://www.ics.uci.edu/~eppstein/numth/frap.c
    237         if ( md <= (T){1} ) {                                   // maximum fractional digits too small?
    238                 return (Rational(T)){ convert( f ), (T){1}}; // truncate fraction
     200        if ( md <= (RationalImpl){1} ) {                                        // maximum fractional digits too small?
     201                return (Rational(RationalImpl)){ convert( f ), (RationalImpl){1}}; // truncate fraction
    239202        } // if
    240203
    241204        // continued fraction coefficients
    242         T m00 = {1}, m11 = { 1 }, m01 = { 0 }, m10 = { 0 };
    243         T ai, t;
     205        RationalImpl m00 = {1}, m11 = { 1 }, m01 = { 0 }, m10 = { 0 };
     206        RationalImpl ai, t;
    244207
    245208        // find terms until denom gets too big
    … …  
    258221          if ( f > (double)0x7FFFFFFF ) break;                          // representation failure
    259222        } // for
    260         return (Rational(T)){ m00, m10 };
     223        return (Rational(RationalImpl)){ m00, m10 };
    261224} // narrow
    262225
Note: See TracChangeset for help on using the changeset viewer.