uClibc now has a math library. muahahahaha!

 -Erik

1 files changed, 309 insertions, 0 deletions

diff --git a/libm/double/polmisc.c b/libm/double/polmisc.c
new file mode 100644
index 000000000..7d517ae69
--- /dev/null
+++ b/libm/double/polmisc.c
@@ -0,0 +1,309 @@
+
+/* Square root, sine, cosine, and arctangent of polynomial.
+ * See polyn.c for data structures and discussion.
+ */
+
+#include <stdio.h>
+#include <math.h>
+#ifdef ANSIPROT
+extern double atan2 ( double, double );
+extern double sqrt ( double );
+extern double fabs ( double );
+extern double sin ( double );
+extern double cos ( double );
+extern void polclr ( double *a, int n );
+extern void polmov ( double *a, int na, double *b );
+extern void polmul ( double a[], int na, double b[], int nb, double c[] );
+extern void poladd ( double a[], int na, double b[], int nb, double c[] );
+extern void polsub ( double a[], int na, double b[], int nb, double c[] );
+extern int poldiv ( double a[], int na, double b[], int nb, double c[] );
+extern void polsbt ( double a[], int na, double b[], int nb, double c[] );
+extern void * malloc ( long );
+extern void free ( void * );
+#else
+double atan2(), sqrt(), fabs(), sin(), cos();
+void polclr(), polmov(), polsbt(), poladd(), polsub(), polmul();
+int poldiv();
+void * malloc();
+void free ();
+#endif
+
+/* Highest degree of polynomial to be handled
+   by the polyn.c subroutine package.  */
+#define N 16
+/* Highest degree actually initialized at runtime.  */
+extern int MAXPOL;
+
+/* Taylor series coefficients for various functions
+ */
+double patan[N+1] = {
+  0.0,     1.0,      0.0, -1.0/3.0,     0.0,
+  1.0/5.0, 0.0, -1.0/7.0,      0.0, 1.0/9.0, 0.0, -1.0/11.0,
+  0.0, 1.0/13.0, 0.0, -1.0/15.0, 0.0 };
+
+double psin[N+1] = {
+  0.0, 1.0, 0.0,   -1.0/6.0,  0.0, 1.0/120.0,  0.0,
+  -1.0/5040.0, 0.0, 1.0/362880.0, 0.0, -1.0/39916800.0,
+  0.0, 1.0/6227020800.0, 0.0, -1.0/1.307674368e12, 0.0};
+
+double pcos[N+1] = {
+  1.0, 0.0,   -1.0/2.0,  0.0, 1.0/24.0,  0.0,
+  -1.0/720.0, 0.0, 1.0/40320.0, 0.0, -1.0/3628800.0, 0.0,
+  1.0/479001600.0, 0.0, -1.0/8.7179291e10, 0.0, 1.0/2.0922789888e13};
+
+double pasin[N+1] = {
+  0.0,     1.0,  0.0, 1.0/6.0,  0.0,
+  3.0/40.0, 0.0, 15.0/336.0, 0.0, 105.0/3456.0, 0.0, 945.0/42240.0,
+  0.0, 10395.0/599040.0 , 0.0, 135135.0/9676800.0 , 0.0
+};
+
+/* Square root of 1 + x.  */
+double psqrt[N+1] = {
+  1.0, 1./2., -1./8., 1./16., -5./128., 7./256., -21./1024., 33./2048.,
+  -429./32768., 715./65536., -2431./262144., 4199./524288., -29393./4194304.,
+  52003./8388608., -185725./33554432., 334305./67108864.,
+  -9694845./2147483648.};
+
+/* Arctangent of the ratio num/den of two polynomials.
+ */
+void
+polatn( num, den, ans, nn )
+     double num[], den[], ans[];
+     int nn;
+{
+  double a, t;
+  double *polq, *polu, *polt;
+  int i;
+
+  if (nn > N)
+    {
+      mtherr ("polatn", OVERFLOW);
+      return;
+    }
+  /* arctan( a + b ) = arctan(a) + arctan( b/(1 + ab + a**2) ) */
+  t = num[0];
+  a = den[0];
+  if( (t == 0.0) && (a == 0.0 ) )
+    {
+      t = num[1];
+      a = den[1];
+    }
+  t = atan2( t, a );  /* arctan(num/den), the ANSI argument order */
+  polq = (double * )malloc( (MAXPOL+1) * sizeof (double) );
+  polu = (double * )malloc( (MAXPOL+1) * sizeof (double) );
+  polt = (double * )malloc( (MAXPOL+1) * sizeof (double) );
+  polclr( polq, MAXPOL );
+  i = poldiv( den, nn, num, nn, polq );
+  a = polq[0]; /* a */
+  polq[0] = 0.0; /* b */
+  polmov( polq, nn, polu ); /* b */
+  /* Form the polynomial
+     1 + ab + a**2
+     where a is a scalar.  */
+  for( i=0; i<=nn; i++ )
+    polu[i] *= a;
+  polu[0] += 1.0 + a * a;
+  poldiv( polu, nn, polq, nn, polt ); /* divide into b */
+  polsbt( polt, nn, patan, nn, polu ); /* arctan(b)  */
+  polu[0] += t; /* plus arctan(a) */
+  polmov( polu, nn, ans );
+  free( polt );
+  free( polu );
+  free( polq );
+}
+
+
+
+/* Square root of a polynomial.
+ * Assumes the lowest degree nonzero term is dominant
+ * and of even degree.  An error message is given
+ * if the Newton iteration does not converge.
+ */
+void
+polsqt( pol, ans, nn )
+     double pol[], ans[];
+     int nn;
+{
+  double t;
+  double *x, *y;
+  int i, n;
+#if 0
+  double z[N+1];
+  double u;
+#endif
+
+  if (nn > N)
+    {
+      mtherr ("polatn", OVERFLOW);
+      return;
+    }
+  x = (double * )malloc( (MAXPOL+1) * sizeof (double) );
+  y = (double * )malloc( (MAXPOL+1) * sizeof (double) );
+  polmov( pol, nn, x );
+  polclr( y, MAXPOL );
+
+  /* Find lowest degree nonzero term.  */
+  t = 0.0;
+  for( n=0; n<nn; n++ )
+    {
+      if( x[n] != 0.0 )
+	goto nzero;
+    }
+  polmov( y, nn, ans );
+  return;
+
+nzero:
+
+  if( n > 0 )
+    {
+      if (n & 1)
+        {
+	  printf("error, sqrt of odd polynomial\n");
+	  return;
+	}
+      /* Divide by x^n.  */
+      y[n] = x[n];
+      poldiv (y, nn, pol, N, x);
+    }
+
+  t = x[0];
+  for( i=1; i<=nn; i++ )
+    x[i] /= t;
+  x[0] = 0.0;
+  /* series development sqrt(1+x) = 1  +  x / 2  -  x**2 / 8  +  x**3 / 16
+     hopes that first (constant) term is greater than what follows   */
+  polsbt( x, nn, psqrt, nn, y);
+  t = sqrt( t );
+  for( i=0; i<=nn; i++ )
+    y[i] *= t;
+
+  /* If first nonzero coefficient was at degree n > 0, multiply by
+     x^(n/2).  */
+  if (n > 0)
+    {
+      polclr (x, MAXPOL);
+      x[n/2] = 1.0;
+      polmul (x, nn, y, nn, y);
+    }
+#if 0
+/* Newton iterations */
+for( n=0; n<10; n++ )
+	{
+	poldiv( y, nn, pol, nn, z );
+	poladd( y, nn, z, nn, y );
+	for( i=0; i<=nn; i++ )
+		y[i] *= 0.5;
+	for( i=0; i<=nn; i++ )
+		{
+		u = fabs( y[i] - z[i] );
+		if( u > 1.0e-15 )
+			goto more;
+		}
+	goto done;
+more:	;
+	}
+printf( "square root did not converge\n" );
+done:
+#endif /* 0 */
+
+polmov( y, nn, ans );
+free( y );
+free( x );
+}
+
+
+
+/* Sine of a polynomial.
+ * The computation uses
+ *     sin(a+b) = sin(a) cos(b) + cos(a) sin(b)
+ * where a is the constant term of the polynomial and
+ * b is the sum of the rest of the terms.
+ * Since sin(b) and cos(b) are computed by series expansions,
+ * the value of b should be small.
+ */
+void
+polsin( x, y, nn )
+     double x[], y[];
+     int nn;
+{
+  double a, sc;
+  double *w, *c;
+  int i;
+
+  if (nn > N)
+    {
+      mtherr ("polatn", OVERFLOW);
+      return;
+    }
+  w = (double * )malloc( (MAXPOL+1) * sizeof (double) );
+  c = (double * )malloc( (MAXPOL+1) * sizeof (double) );
+  polmov( x, nn, w );
+  polclr( c, MAXPOL );
+  polclr( y, nn );
+  /* a, in the description, is x[0].  b is the polynomial x - x[0].  */
+  a = w[0];
+  /* c = cos (b) */
+  w[0] = 0.0;
+  polsbt( w, nn, pcos, nn, c );
+  sc = sin(a);
+  /* sin(a) cos (b) */
+  for( i=0; i<=nn; i++ )
+    c[i] *= sc;
+  /* y = sin (b)  */
+  polsbt( w, nn, psin, nn, y );
+  sc = cos(a);
+  /* cos(a) sin(b) */
+  for( i=0; i<=nn; i++ )
+    y[i] *= sc;
+  poladd( c, nn, y, nn, y );
+  free( c );
+  free( w );
+}
+
+
+/* Cosine of a polynomial.
+ * The computation uses
+ *     cos(a+b) = cos(a) cos(b) - sin(a) sin(b)
+ * where a is the constant term of the polynomial and
+ * b is the sum of the rest of the terms.
+ * Since sin(b) and cos(b) are computed by series expansions,
+ * the value of b should be small.
+ */
+void
+polcos( x, y, nn )
+     double x[], y[];
+     int nn;
+{
+  double a, sc;
+  double *w, *c;
+  int i;
+  double sin(), cos();
+
+  if (nn > N)
+    {
+      mtherr ("polatn", OVERFLOW);
+      return;
+    }
+  w = (double * )malloc( (MAXPOL+1) * sizeof (double) );
+  c = (double * )malloc( (MAXPOL+1) * sizeof (double) );
+  polmov( x, nn, w );
+  polclr( c, MAXPOL );
+  polclr( y, nn );
+  a = w[0];
+  w[0] = 0.0;
+  /* c = cos(b)  */
+  polsbt( w, nn, pcos, nn, c );
+  sc = cos(a);
+  /* cos(a) cos(b)  */
+  for( i=0; i<=nn; i++ )
+    c[i] *= sc;
+  /* y = sin(b) */
+  polsbt( w, nn, psin, nn, y );
+  sc = sin(a);
+  /* sin(a) sin(b) */
+  for( i=0; i<=nn; i++ )
+    y[i] *= sc;
+  polsub( y, nn, c, nn, y );
+  free( c );
+  free( w );
+}