uClibc now has a math library. muahahahaha!

 -Erik

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 );
+}
+