summaryrefslogtreecommitdiff
path: root/libm/double/euclid.c
diff options
context:
space:
mode:
authorEric Andersen <andersen@codepoet.org>2001-11-22 14:04:29 +0000
committerEric Andersen <andersen@codepoet.org>2001-11-22 14:04:29 +0000
commit7ce331c01ce6eb7b3f5c715a38a24359da9c6ee2 (patch)
tree3a7e8476e868ae15f4da1b7ce26b2db6f434468c /libm/double/euclid.c
parentc117dd5fb183afb1a4790a6f6110d88704be6bf8 (diff)
Totally rework the math library, this time based on the MacOs X
math library (which is itself based on the math lib from FreeBSD). -Erik
Diffstat (limited to 'libm/double/euclid.c')
-rw-r--r--libm/double/euclid.c251
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 );
-}
-