diff options
Diffstat (limited to 'libm/double/euclid.c')
| -rw-r--r-- | libm/double/euclid.c | 251 |
1 files changed, 0 insertions, 251 deletions
diff --git a/libm/double/euclid.c b/libm/double/euclid.c deleted file mode 100644 index 3a899a6d2..000000000 --- a/libm/double/euclid.c +++ /dev/null @@ -1,251 +0,0 @@ -/* euclid.c - * - * Rational arithmetic routines - * - * - * - * SYNOPSIS: - * - * - * typedef struct - * { - * double n; numerator - * double d; denominator - * }fract; - * - * radd( a, b, c ) c = b + a - * rsub( a, b, c ) c = b - a - * rmul( a, b, c ) c = b * a - * rdiv( a, b, c ) c = b / a - * euclid( &n, &d ) Reduce n/d to lowest terms, - * return greatest common divisor. - * - * Arguments of the routines are pointers to the structures. - * The double precision numbers are assumed, without checking, - * to be integer valued. Overflow conditions are reported. - */ - - -#include <math.h> -#ifdef ANSIPROT -extern double fabs ( double ); -extern double floor ( double ); -double euclid( double *, double * ); -#else -double fabs(), floor(), euclid(); -#endif - -extern double MACHEP; -#define BIG (1.0/MACHEP) - -typedef struct - { - double n; /* numerator */ - double d; /* denominator */ - }fract; - -/* Add fractions. */ - -void radd( f1, f2, f3 ) -fract *f1, *f2, *f3; -{ -double gcd, d1, d2, gcn, n1, n2; - -n1 = f1->n; -d1 = f1->d; -n2 = f2->n; -d2 = f2->d; -if( n1 == 0.0 ) - { - f3->n = n2; - f3->d = d2; - return; - } -if( n2 == 0.0 ) - { - f3->n = n1; - f3->d = d1; - return; - } - -gcd = euclid( &d1, &d2 ); /* common divisors of denominators */ -gcn = euclid( &n1, &n2 ); /* common divisors of numerators */ -/* Note, factoring the numerators - * makes overflow slightly less likely. - */ -f3->n = ( n1 * d2 + n2 * d1) * gcn; -f3->d = d1 * d2 * gcd; -euclid( &f3->n, &f3->d ); -} - - -/* Subtract fractions. */ - -void rsub( f1, f2, f3 ) -fract *f1, *f2, *f3; -{ -double gcd, d1, d2, gcn, n1, n2; - -n1 = f1->n; -d1 = f1->d; -n2 = f2->n; -d2 = f2->d; -if( n1 == 0.0 ) - { - f3->n = n2; - f3->d = d2; - return; - } -if( n2 == 0.0 ) - { - f3->n = -n1; - f3->d = d1; - return; - } - -gcd = euclid( &d1, &d2 ); -gcn = euclid( &n1, &n2 ); -f3->n = (n2 * d1 - n1 * d2) * gcn; -f3->d = d1 * d2 * gcd; -euclid( &f3->n, &f3->d ); -} - - - - -/* Multiply fractions. */ - -void rmul( ff1, ff2, ff3 ) -fract *ff1, *ff2, *ff3; -{ -double d1, d2, n1, n2; - -n1 = ff1->n; -d1 = ff1->d; -n2 = ff2->n; -d2 = ff2->d; - -if( (n1 == 0.0) || (n2 == 0.0) ) - { - ff3->n = 0.0; - ff3->d = 1.0; - return; - } -euclid( &n1, &d2 ); /* cross cancel common divisors */ -euclid( &n2, &d1 ); -ff3->n = n1 * n2; -ff3->d = d1 * d2; -/* Report overflow. */ -if( (fabs(ff3->n) >= BIG) || (fabs(ff3->d) >= BIG) ) - { - mtherr( "rmul", OVERFLOW ); - return; - } -/* euclid( &ff3->n, &ff3->d );*/ -} - - - -/* Divide fractions. */ - -void rdiv( ff1, ff2, ff3 ) -fract *ff1, *ff2, *ff3; -{ -double d1, d2, n1, n2; - -n1 = ff1->d; /* Invert ff1, then multiply */ -d1 = ff1->n; -if( d1 < 0.0 ) - { /* keep denominator positive */ - n1 = -n1; - d1 = -d1; - } -n2 = ff2->n; -d2 = ff2->d; -if( (n1 == 0.0) || (n2 == 0.0) ) - { - ff3->n = 0.0; - ff3->d = 1.0; - return; - } - -euclid( &n1, &d2 ); /* cross cancel any common divisors */ -euclid( &n2, &d1 ); -ff3->n = n1 * n2; -ff3->d = d1 * d2; -/* Report overflow. */ -if( (fabs(ff3->n) >= BIG) || (fabs(ff3->d) >= BIG) ) - { - mtherr( "rdiv", OVERFLOW ); - return; - } -/* euclid( &ff3->n, &ff3->d );*/ -} - - - - - -/* Euclidean algorithm - * reduces fraction to lowest terms, - * returns greatest common divisor. - */ - - -double euclid( num, den ) -double *num, *den; -{ -double n, d, q, r; - -n = *num; /* Numerator. */ -d = *den; /* Denominator. */ - -/* Make numbers positive, locally. */ -if( n < 0.0 ) - n = -n; -if( d < 0.0 ) - d = -d; - -/* Abort if numbers are too big for integer arithmetic. */ -if( (n >= BIG) || (d >= BIG) ) - { - mtherr( "euclid", OVERFLOW ); - return(1.0); - } - -/* Divide by zero, gcd = 1. */ -if(d == 0.0) - return( 1.0 ); - -/* Zero. Return 0/1, gcd = denominator. */ -if(n == 0.0) - { -/* - if( *den < 0.0 ) - *den = -1.0; - else - *den = 1.0; -*/ - *den = 1.0; - return( d ); - } - -while( d > 0.5 ) - { -/* Find integer part of n divided by d. */ - q = floor( n/d ); -/* Find remainder after dividing n by d. */ - r = n - d * q; -/* The next fraction is d/r. */ - n = d; - d = r; - } - -if( n < 0.0 ) - mtherr( "euclid", UNDERFLOW ); - -*num /= n; -*den /= n; -return( n ); -} - |