1 files changed, 520 insertions, 0 deletions
diff --git a/libm/float/polynf.c b/libm/float/polynf.c
new file mode 100644
index 000000000..48c6675d4
--- /dev/null
+++ b/libm/float/polynf.c
@@ -0,0 +1,520 @@
+/*							polynf.c
+ *							polyrf.c
+ * Arithmetic operations on polynomials
+ *
+ * In the following descriptions a, b, c are polynomials of degree
+ * na, nb, nc respectively.  The degree of a polynomial cannot
+ * exceed a run-time value MAXPOLF.  An operation that attempts
+ * to use or generate a polynomial of higher degree may produce a
+ * result that suffers truncation at degree MAXPOL.  The value of
+ * MAXPOL is set by calling the function
+ *
+ *     polinif( maxpol );
+ *
+ * where maxpol is the desired maximum degree.  This must be
+ * done prior to calling any of the other functions in this module.
+ * Memory for internal temporary polynomial storage is allocated
+ * by polinif().
+ *
+ * Each polynomial is represented by an array containing its
+ * coefficients, together with a separately declared integer equal
+ * to the degree of the polynomial.  The coefficients appear in
+ * ascending order; that is,
+ *
+ *                                        2                      na
+ * a(x)  =  a[0]  +  a[1] * x  +  a[2] * x   +  ...  +  a[na] * x  .
+ *
+ *
+ *
+ * sum = poleva( a, na, x );	Evaluate polynomial a(t) at t = x.
+ * polprtf( a, na, D );		Print the coefficients of a to D digits.
+ * polclrf( a, na );		Set a identically equal to zero, up to a[na].
+ * polmovf( a, na, b );		Set b = a.
+ * poladdf( a, na, b, nb, c );	c = b + a, nc = max(na,nb)
+ * polsubf( a, na, b, nb, c );	c = b - a, nc = max(na,nb)
+ * polmulf( a, na, b, nb, c );	c = b * a, nc = na+nb
+ *
+ *
+ * Division:
+ *
+ * i = poldivf( a, na, b, nb, c );	c = b / a, nc = MAXPOL
+ *
+ * returns i = the degree of the first nonzero coefficient of a.
+ * The computed quotient c must be divided by x^i.  An error message
+ * is printed if a is identically zero.
+ *
+ *
+ * Change of variables:
+ * If a and b are polynomials, and t = a(x), then
+ *     c(t) = b(a(x))
+ * is a polynomial found by substituting a(x) for t.  The
+ * subroutine call for this is
+ *
+ * polsbtf( a, na, b, nb, c );
+ *
+ *
+ * Notes:
+ * poldivf() is an integer routine; polevaf() is float.
+ * Any of the arguments a, b, c may refer to the same array.
+ *
+ */
+
+#ifndef NULL
+#define NULL 0
+#endif
+#include <math.h>
+
+#ifdef ANSIC
+void printf(), sprintf(), exit();
+void free(void *);
+void *malloc(int);
+#else
+void printf(), sprintf(), free(), exit();
+void *malloc();
+#endif
+/* near pointer version of malloc() */
+/*#define malloc _nmalloc*/
+/*#define free _nfree*/
+
+/* Pointers to internal arrays.  Note poldiv() allocates
+ * and deallocates some temporary arrays every time it is called.
+ */
+static float *pt1 = 0;
+static float *pt2 = 0;
+static float *pt3 = 0;
+
+/* Maximum degree of polynomial. */
+int MAXPOLF = 0;
+extern int MAXPOLF;
+
+/* Number of bytes (chars) in maximum size polynomial. */
+static int psize = 0;
+
+
+/* Initialize max degree of polynomials
+ * and allocate temporary storage.
+ */
+#ifdef ANSIC
+void polinif( int maxdeg )
+#else
+int polinif( maxdeg )
+int maxdeg;
+#endif
+{
+
+MAXPOLF = maxdeg;
+psize = (maxdeg + 1) * sizeof(float);
+
+/* Release previously allocated memory, if any. */
+if( pt3 )
+	free(pt3);
+if( pt2 )
+	free(pt2);
+if( pt1 )
+	free(pt1);
+
+/* Allocate new arrays */
+pt1 = (float * )malloc(psize); /* used by polsbtf */
+pt2 = (float * )malloc(psize); /* used by polsbtf */
+pt3 = (float * )malloc(psize); /* used by polmul */
+
+/* Report if failure */
+if( (pt1 == NULL) || (pt2 == NULL) || (pt3 == NULL) )
+	{
+	mtherr( "polinif", ERANGE );
+	exit(1);
+	}
+#if !ANSIC
+return 0;
+#endif
+}
+
+
+
+/* Print the coefficients of a, with d decimal precision.
+ */
+static char *form = "abcdefghijk";
+
+#ifdef ANSIC
+void polprtf( float *a, int na, int d )
+#else
+int polprtf( a, na, d )
+float a[];
+int na, d;
+#endif
+{
+int i, j, d1;
+char *p;
+
+/* Create format descriptor string for the printout.
+ * Do this partly by hand, since sprintf() may be too
+ * bug-ridden to accomplish this feat by itself.
+ */
+p = form;
+*p++ = '%';
+d1 = d + 8;
+(void )sprintf( p, "%d ", d1 );
+p += 1;
+if( d1 >= 10 )
+	p += 1;
+*p++ = '.';
+(void )sprintf( p, "%d ", d );
+p += 1;
+if( d >= 10 )
+	p += 1;
+*p++ = 'e';
+*p++ = ' ';
+*p++ = '\0';
+
+
+/* Now do the printing.
+ */
+d1 += 1;
+j = 0;
+for( i=0; i<=na; i++ )
+	{
+/* Detect end of available line */
+	j += d1;
+	if( j >= 78 )
+		{
+		printf( "\n" );
+		j = d1;
+		}
+	printf( form, a[i] );
+	}
+printf( "\n" );
+#if !ANSIC
+return 0;
+#endif
+}
+
+
+
+/* Set a = 0.
+ */
+#ifdef ANSIC
+void polclrf( register float *a, int n )
+#else
+int polclrf( a, n )
+register float *a;
+int n;
+#endif
+{
+int i;
+
+if( n > MAXPOLF )
+	n = MAXPOLF;
+for( i=0; i<=n; i++ )
+	*a++ = 0.0;
+#if !ANSIC
+return 0;
+#endif
+}
+
+
+
+/* Set b = a.
+ */
+#ifdef ANSIC
+void polmovf( register float *a, int na, register float *b )
+#else
+int polmovf( a, na, b )
+register float *a, *b;
+int na;
+#endif
+{
+int i;
+
+if( na > MAXPOLF )
+	na = MAXPOLF;
+
+for( i=0; i<= na; i++ )
+	{
+	*b++ = *a++;
+	}
+#if !ANSIC
+return 0;
+#endif
+}
+
+
+/* c = b * a.
+ */
+#ifdef ANSIC
+void polmulf( float a[], int na, float b[], int nb, float c[] )
+#else
+int polmulf( a, na, b, nb, c )
+float a[], b[], c[];
+int na, nb;
+#endif
+{
+int i, j, k, nc;
+float x;
+
+nc = na + nb;
+polclrf( pt3, MAXPOLF );
+
+for( i=0; i<=na; i++ )
+	{
+	x = a[i];
+	for( j=0; j<=nb; j++ )
+		{
+		k = i + j;
+		if( k > MAXPOLF )
+			break;
+		pt3[k] += x * b[j];
+		}
+	}
+
+if( nc > MAXPOLF )
+	nc = MAXPOLF;
+for( i=0; i<=nc; i++ )
+	c[i] = pt3[i];
+#if !ANSIC
+return 0;
+#endif
+}
+
+
+
+ 
+/* c = b + a.
+ */
+#ifdef ANSIC
+void poladdf( float a[], int na, float b[], int nb, float c[] )
+#else
+int poladdf( a, na, b, nb, c )
+float a[], b[], c[];
+int na, nb;
+#endif
+{
+int i, n;
+
+
+if( na > nb )
+	n = na;
+else
+	n = nb;
+
+if( n > MAXPOLF )
+	n = MAXPOLF;
+
+for( i=0; i<=n; i++ )
+	{
+	if( i > na )
+		c[i] = b[i];
+	else if( i > nb )
+		c[i] = a[i];
+	else
+		c[i] = b[i] + a[i];
+	}
+#if !ANSIC
+return 0;
+#endif
+}
+
+/* c = b - a.
+ */
+#ifdef ANSIC
+void polsubf( float a[], int na, float b[], int nb, float c[] )
+#else
+int polsubf( a, na, b, nb, c )
+float a[], b[], c[];
+int na, nb;
+#endif
+{
+int i, n;
+
+
+if( na > nb )
+	n = na;
+else
+	n = nb;
+
+if( n > MAXPOLF )
+	n = MAXPOLF;
+
+for( i=0; i<=n; i++ )
+	{
+	if( i > na )
+		c[i] = b[i];
+	else if( i > nb )
+		c[i] = -a[i];
+	else
+		c[i] = b[i] - a[i];
+	}
+#if !ANSIC
+return 0;
+#endif
+}
+
+
+
+/* c = b/a
+ */
+#ifdef ANSIC
+int poldivf( float a[], int na, float b[], int nb, float c[] )
+#else
+int poldivf( a, na, b, nb, c )
+float a[], b[], c[];
+int na, nb;
+#endif
+{
+float quot;
+float *ta, *tb, *tq;
+int i, j, k, sing;
+
+sing = 0;
+
+/* Allocate temporary arrays.  This would be quicker
+ * if done automatically on the stack, but stack space
+ * may be hard to obtain on a small computer.
+ */
+ta = (float * )malloc( psize );
+polclrf( ta, MAXPOLF );
+polmovf( a, na, ta );
+
+tb = (float * )malloc( psize );
+polclrf( tb, MAXPOLF );
+polmovf( b, nb, tb );
+
+tq = (float * )malloc( psize );
+polclrf( tq, MAXPOLF );
+
+/* What to do if leading (constant) coefficient
+ * of denominator is zero.
+ */
+if( a[0] == 0.0 )
+	{
+	for( i=0; i<=na; i++ )
+		{
+		if( ta[i] != 0.0 )
+			goto nzero;
+		}
+	mtherr( "poldivf", SING );
+	goto done;
+
+nzero:
+/* Reduce the degree of the denominator. */
+	for( i=0; i<na; i++ )
+		ta[i] = ta[i+1];
+	ta[na] = 0.0;
+
+	if( b[0] != 0.0 )
+		{
+/* Optional message:
+		printf( "poldivf singularity, divide quotient by x\n" );
+*/
+		sing += 1;
+		}
+	else
+		{
+/* Reduce degree of numerator. */
+		for( i=0; i<nb; i++ )
+			tb[i] = tb[i+1];
+		tb[nb] = 0.0;
+		}
+/* Call self, using reduced polynomials. */
+	sing += poldivf( ta, na, tb, nb, c );
+	goto done;
+	}
+
+/* Long division algorithm.  ta[0] is nonzero.
+ */
+for( i=0; i<=MAXPOLF; i++ )
+	{
+	quot = tb[i]/ta[0];
+	for( j=0; j<=MAXPOLF; j++ )
+		{
+		k = j + i;
+		if( k > MAXPOLF )
+			break;
+		tb[k] -= quot * ta[j];
+		}
+	tq[i] = quot;
+	}
+/* Send quotient to output array. */
+polmovf( tq, MAXPOLF, c );
+
+done:
+
+/* Restore allocated memory. */
+free(tq);
+free(tb);
+free(ta);
+return( sing );
+}
+
+
+
+
+/* Change of variables
+ * Substitute a(y) for the variable x in b(x).
+ * x = a(y)
+ * c(x) = b(x) = b(a(y)).
+ */
+
+#ifdef ANSIC
+void polsbtf( float a[], int na, float b[], int nb, float c[] )
+#else
+int polsbtf( a, na, b, nb, c )
+float a[], b[], c[];
+int na, nb;
+#endif
+{
+int i, j, k, n2;
+float x;
+
+/* 0th degree term:
+ */
+polclrf( pt1, MAXPOLF );
+pt1[0] = b[0];
+
+polclrf( pt2, MAXPOLF );
+pt2[0] = 1.0;
+n2 = 0;
+
+for( i=1; i<=nb; i++ )
+	{
+/* Form ith power of a. */
+	polmulf( a, na, pt2, n2, pt2 );
+	n2 += na;
+	x = b[i];
+/* Add the ith coefficient of b times the ith power of a. */
+	for( j=0; j<=n2; j++ )
+		{
+		if( j > MAXPOLF )
+			break;
+		pt1[j] += x * pt2[j];
+		}
+	}
+
+k = n2 + nb;
+if( k > MAXPOLF )
+	k = MAXPOLF;
+for( i=0; i<=k; i++ )
+	c[i] = pt1[i];
+#if !ANSIC
+return 0;
+#endif
+}
+
+
+
+
+/* Evaluate polynomial a(t) at t = x.
+ */
+float polevaf( float *a, int na, float xx )
+{
+float x, s;
+int i;
+
+x = xx;
+s = a[na];
+for( i=na-1; i>=0; i-- )
+	{
+	s = s * x + a[i];
+	}
+return(s);
+}
+