summaryrefslogtreecommitdiff
path: root/libm/double/euclid.c
diff options
context:
space:
mode:
authorEric Andersen <andersen@codepoet.org>2001-05-10 00:40:28 +0000
committerEric Andersen <andersen@codepoet.org>2001-05-10 00:40:28 +0000
commit1077fa4d772832f77a677ce7fb7c2d513b959e3f (patch)
tree579bee13fb0b58d2800206366ec2caecbb15f3fc /libm/double/euclid.c
parent22358dd7ce7bb49792204b698f01a6f69b9c8e08 (diff)
uClibc now has a math library. muahahahaha!
-Erik
Diffstat (limited to 'libm/double/euclid.c')
-rw-r--r--libm/double/euclid.c251
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 );
+}
+