diff options
Diffstat (limited to 'libm/double/euclid.c')
| -rw-r--r-- | libm/double/euclid.c | 251 |
1 files changed, 251 insertions, 0 deletions
diff --git a/libm/double/euclid.c b/libm/double/euclid.c new file mode 100644 index 000000000..3a899a6d2 --- /dev/null +++ b/libm/double/euclid.c @@ -0,0 +1,251 @@ +/* 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 ); +} + |