Changes in src/libcfa/rational.c [6e4b913:a6151ba]

File:

• src/libcfa/rational.c (modified) (9 diffs)

src/libcfa/rational.c

@@ -10 +10 @@
 // Created On       : Wed Apr  6 17:54:28 2016
 // Last Modified By : Peter A. Buhr
-// Last Modified On : Sat Jul  9 11:18:04 2016
-// Update Count     : 40
+// Last Modified On : Tue Jul  5 18:29:12 2016
+// Update Count     : 26
 // 
 
@@ -31 +31 @@
 // alternative: https://en.wikipedia.org/wiki/Binary_GCD_algorithm
 static long int gcd( long int a, long int b ) {
-        for ( ;; ) {                                                                            // Euclid's algorithm
+    for ( ;; ) {                                                                                // Euclid's algorithm
                 long int r = a % b;
           if ( r == 0 ) break;
                 a = b;
                 b = r;
-        } // for
+    } // for
         return b;
 } // gcd
 
 static long int simplify( long int *n, long int *d ) {
-        if ( *d == 0 ) {
+    if ( *d == 0 ) {
                 serr | "Invalid rational number construction: denominator cannot be equal to 0." | endl;
                 exit( EXIT_FAILURE );
-        } // exit
-        if ( *d < 0 ) { *d = -*d; *n = -*n; }                           // move sign to numerator
-        return gcd( abs( *n ), *d );                                            // simplify
+    } // exit
+    if ( *d < 0 ) { *d = -*d; *n = -*n; }                               // move sign to numerator
+    return gcd( abs( *n ), *d );                                                // simplify
 } // Rationalnumber::simplify
 
@@ -53 +53 @@
 
 void ?{}( Rational * r ) {
-        r{ 0, 1 };
+    r{ 0, 1 };
 } // rational
 
 void ?{}( Rational * r, long int n ) {
-        r{ n, 1 };
+    r{ n, 1 };
 } // rational
 
 void ?{}( Rational * r, long int n, long int d ) {
-        long int t = simplify( &n, &d );                                        // simplify
-        r->numerator = n / t;
+    long int t = simplify( &n, &d );                                    // simplify
+    r->numerator = n / t;
         r->denominator = d / t;
 } // rational
@@ -70 +70 @@
 
 long int numerator( Rational r ) {
-        return r.numerator;
+    return r.numerator;
 } // numerator
 
 long int numerator( Rational r, long int n ) {
-        long int prev = r.numerator;
-        long int t = gcd( abs( n ), r.denominator );            // simplify
-        r.numerator = n / t;
-        r.denominator = r.denominator / t;
-        return prev;
+    long int prev = r.numerator;
+    long int t = gcd( abs( n ), r.denominator );                // simplify
+    r.numerator = n / t;
+    r.denominator = r.denominator / t;
+    return prev;
 } // numerator
 
 long int denominator( Rational r ) {
-        return r.denominator;
+    return r.denominator;
 } // denominator
 
 long int denominator( Rational r, long int d ) {
-        long int prev = r.denominator;
-        long int t = simplify( &r.numerator, &d );                      // simplify
-        r.numerator = r.numerator / t;
-        r.denominator = d / t;
-        return prev;
+    long int prev = r.denominator;
+    long int t = simplify( &r.numerator, &d );                  // simplify
+    r.numerator = r.numerator / t;
+    r.denominator = d / t;
+    return prev;
 } // denominator
 
@@ -97 +97 @@
 
 int ?==?( Rational l, Rational r ) {
-        return l.numerator * r.denominator == l.denominator * r.numerator;
+    return l.numerator * r.denominator == l.denominator * r.numerator;
 } // ?==?
 
 int ?!=?( Rational l, Rational r ) {
-        return ! ( l == r );
+    return ! ( l == r );
 } // ?!=?
 
 int ?<?( Rational l, Rational r ) {
-        return l.numerator * r.denominator < l.denominator * r.numerator;
+    return l.numerator * r.denominator < l.denominator * r.numerator;
 } // ?<?
 
 int ?<=?( Rational l, Rational r ) {
-        return l < r || l == r;
+    return l < r || l == r;
 } // ?<=?
 
 int ?>?( Rational l, Rational r ) {
-        return ! ( l <= r );
+    return ! ( l <= r );
 } // ?>?
 
 int ?>=?( Rational l, Rational r ) {
-        return ! ( l < r );
+    return ! ( l < r );
 } // ?>=?
 
@@ -125 +125 @@
 Rational -?( Rational r ) {
         Rational t = { -r.numerator, r.denominator };
-        return t;
+    return t;
 } // -?
 
 Rational ?+?( Rational l, Rational r ) {
-        if ( l.denominator == r.denominator ) {                         // special case
+    if ( l.denominator == r.denominator ) {                             // special case
                 Rational t = { l.numerator + r.numerator, l.denominator };
                 return t;
-        } else {
+    } else {
                 Rational t = { l.numerator * r.denominator + l.denominator * r.numerator, l.denominator * r.denominator };
                 return t;
-        } // if
+    } // if
 } // ?+?
 
 Rational ?-?( Rational l, Rational r ) {
-        if ( l.denominator == r.denominator ) {                         // special case
+    if ( l.denominator == r.denominator ) {                             // special case
                 Rational t = { l.numerator - r.numerator, l.denominator };
                 return t;
-        } else {
+    } else {
                 Rational t = { l.numerator * r.denominator - l.denominator * r.numerator, l.denominator * r.denominator };
                 return t;
-        } // if
+    } // if
 } // ?-?
 
 Rational ?*?( Rational l, Rational r ) {
-        Rational t = { l.numerator * r.numerator, l.denominator * r.denominator };
+    Rational t = { l.numerator * r.numerator, l.denominator * r.denominator };
         return t;
 } // ?*?
 
 Rational ?/?( Rational l, Rational r ) {
-        if ( r.numerator < 0 ) {
+    if ( r.numerator < 0 ) {
                 r.numerator = -r.numerator;
                 r.denominator = -r.denominator;
         } // if
         Rational t = { l.numerator * r.denominator, l.denominator * r.numerator };
-        return t;
+    return t;
 } // ?/?
 
@@ -169 +169 @@
 } // widen
 
-// http://www.ics.uci.edu/~eppstein/numth/frap.c
+// https://rosettacode.org/wiki/Convert_decimal_number_to_rational#C
 Rational narrow( double f, long int md ) {
         if ( md <= 1 ) {                                                                        // maximum fractional digits too small?
@@ -176 +176 @@
 
         // continued fraction coefficients
-        long int m00 = 1, m11 = 1, m01 = 0, m10 = 0;
-        long int ai, t;
-
-        // find terms until denom gets too big
-        for ( ;; ) {
-                ai = (long int)f;
-          if ( ! (m10 * ai + m11 <= md) ) break;
-                t = m00 * ai + m01;
-                m01 = m00;
-                m00 = t;
-                t = m10 * ai + m11;
-                m11 = m10;
-                m10 = t;
-                t = (double)ai;
-          if ( f == t ) break;                                                          // prevent division by zero
-                f = 1 / (f - t);
-          if ( f > (double)0x7FFFFFFF ) break;                          // representation failure
-        } 
-        return (Rational){ m00, m10 };
+        long int a, h[3] = { 0, 1, 0 }, k[3] = { 1, 0, 0 };
+        long int x, d, n = 1;
+        int i, neg = 0;
+ 
+        if ( f < 0 ) { neg = 1; f = -f; }
+        while ( f != floor( f ) ) { n <<= 1; f *= 2; }
+        d = f;
+ 
+        // continued fraction and check denominator each step
+        for (i = 0; i < 64; i++) {
+                a = n ? d / n : 0;
+          if (i && !a) break;
+                x = d; d = n; n = x % n;
+                x = a;
+                if (k[1] * a + k[0] >= md) {
+                        x = (md - k[0]) / k[1];
+                  if ( ! (x * 2 >= a || k[1] >= md) ) break;
+                        i = 65;
+                } // if
+                h[2] = x * h[1] + h[0]; h[0] = h[1]; h[1] = h[2];
+                k[2] = x * k[1] + k[0]; k[0] = k[1]; k[1] = k[2];
+        } // for
+        return (Rational){ neg ? -h[1] : h[1], k[1] };
 } // narrow
 
@@ -203 +207 @@
 istype * ?|?( istype *is, Rational *r ) {
         long int t;
-        is | &(r->numerator) | &(r->denominator);
+    is | &(r->numerator) | &(r->denominator);
         t = simplify( &(r->numerator), &(r->denominator) );
-        r->numerator /= t;
-        r->denominator /= t;
-        return is;
+    r->numerator /= t;
+    r->denominator /= t;
+    return is;
 } // ?|?
 
 forall( dtype ostype | ostream( ostype ) )
 ostype * ?|?( ostype *os, Rational r ) {
-        return os | r.numerator | '/' | r.denominator;
+    return os | r.numerator | '/' | r.denominator;
 } // ?|?