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

1 files changed, 0 insertions, 533 deletions

diff --git a/libm/double/polyr.c b/libm/double/polyr.c
deleted file mode 100644
index 81ca817e3..000000000
--- a/libm/double/polyr.c
+++ /dev/null
@@ -1,533 +0,0 @@
-
-/* Arithmetic operations on polynomials with rational coefficients
- *
- * 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 MAXPOL.  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
- *
- *     polini( 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 polini().
- *
- * 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  .
- *
- *
- *
- * `a', `b', `c' are arrays of fracts.
- * poleva( a, na, &x, &sum );	Evaluate polynomial a(t) at t = x.
- * polprt( a, na, D );		Print the coefficients of a to D digits.
- * polclr( a, na );		Set a identically equal to zero, up to a[na].
- * polmov( a, na, b );		Set b = a.
- * poladd( a, na, b, nb, c );	c = b + a, nc = max(na,nb)
- * polsub( a, na, b, nb, c );	c = b - a, nc = max(na,nb)
- * polmul( a, na, b, nb, c );	c = b * a, nc = na+nb
- *
- *
- * Division:
- *
- * i = poldiv( 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
- *
- * polsbt( a, na, b, nb, c );
- *
- *
- * Notes:
- * poldiv() is an integer routine; poleva() is double.
- * Any of the arguments a, b, c may refer to the same array.
- *
- */
-
-#include <stdio.h>
-#include <math.h>
-#ifndef NULL
-#define NULL 0
-#endif
-typedef struct{
-	double n;
-	double d;
-	}fract;
-
-#ifdef ANSIPROT
-extern void radd ( fract *, fract *, fract * );
-extern void rsub ( fract *, fract *, fract * );
-extern void rmul ( fract *, fract *, fract * );
-extern void rdiv ( fract *, fract *, fract * );
-void polmov ( fract *, int, fract * );
-void polmul ( fract *, int, fract *, int, fract * );
-int poldiv ( fract *, int, fract *, int, fract * );
-void * malloc ( long );
-void free ( void * );
-#else
-void radd(), rsub(), rmul(), rdiv();
-void polmov(), polmul();
-int poldiv();
-void * malloc();
-void free ();
-#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 fract *pt1 = 0;
-static fract *pt2 = 0;
-static fract *pt3 = 0;
-
-/* Maximum degree of polynomial. */
-int MAXPOL = 0;
-extern int MAXPOL;
-
-/* Number of bytes (chars) in maximum size polynomial. */
-static int psize = 0;
-
-
-/* Initialize max degree of polynomials
- * and allocate temporary storage.
- */
-void polini( maxdeg )
-int maxdeg;
-{
-
-MAXPOL = maxdeg;
-psize = (maxdeg + 1) * sizeof(fract);
-
-/* Release previously allocated memory, if any. */
-if( pt3 )
-	free(pt3);
-if( pt2 )
-	free(pt2);
-if( pt1 )
-	free(pt1);
-
-/* Allocate new arrays */
-pt1 = (fract * )malloc(psize); /* used by polsbt */
-pt2 = (fract * )malloc(psize); /* used by polsbt */
-pt3 = (fract * )malloc(psize); /* used by polmul */
-
-/* Report if failure */
-if( (pt1 == NULL) || (pt2 == NULL) || (pt3 == NULL) )
-	{
-	mtherr( "polini", ERANGE );
-	exit(1);
-	}
-}
-
-
-
-/* Print the coefficients of a, with d decimal precision.
- */
-static char *form = "abcdefghijk";
-
-void polprt( a, na, d )
-fract a[];
-int na, d;
-{
-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;
-sprintf( p, "%d ", d1 );
-p += 1;
-if( d1 >= 10 )
-	p += 1;
-*p++ = '.';
-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].n );
-	j += d1;
-	if( j >= 78 )
-		{
-		printf( "\n" );
-		j = d1;
-		}
-	printf( form, a[i].d );
-	}
-printf( "\n" );
-}
-
-
-
-/* Set a = 0.
- */
-void polclr( a, n )
-fract a[];
-int n;
-{
-int i;
-
-if( n > MAXPOL )
-	n = MAXPOL;
-for( i=0; i<=n; i++ )
-	{
-	a[i].n = 0.0;
-	a[i].d = 1.0;
-	}
-}
-
-
-
-/* Set b = a.
- */
-void polmov( a, na, b )
-fract a[], b[];
-int na;
-{
-int i;
-
-if( na > MAXPOL )
-	na = MAXPOL;
-
-for( i=0; i<= na; i++ )
-	{
-	b[i].n = a[i].n;
-	b[i].d = a[i].d;
-	}
-}
-
-
-/* c = b * a.
- */
-void polmul( a, na, b, nb, c )
-fract a[], b[], c[];
-int na, nb;
-{
-int i, j, k, nc;
-fract temp;
-fract *p;
-
-nc = na + nb;
-polclr( pt3, MAXPOL );
-
-p = &a[0];
-for( i=0; i<=na; i++ )
-	{
-	for( j=0; j<=nb; j++ )
-		{
-		k = i + j;
-		if( k > MAXPOL )
-			break;
-		rmul( p, &b[j], &temp ); /*pt3[k] += a[i] * b[j];*/
-		radd( &temp, &pt3[k], &pt3[k] );
-		}
-	++p;
-	}
-
-if( nc > MAXPOL )
-	nc = MAXPOL;
-for( i=0; i<=nc; i++ )
-	{
-	c[i].n = pt3[i].n;
-	c[i].d = pt3[i].d;
-	}
-}
-
-
-
- 
-/* c = b + a.
- */
-void poladd( a, na, b, nb, c )
-fract a[], b[], c[];
-int na, nb;
-{
-int i, n;
-
-
-if( na > nb )
-	n = na;
-else
-	n = nb;
-
-if( n > MAXPOL )
-	n = MAXPOL;
-
-for( i=0; i<=n; i++ )
-	{
-	if( i > na )
-		{
-		c[i].n = b[i].n;
-		c[i].d = b[i].d;
-		}
-	else if( i > nb )
-		{
-		c[i].n = a[i].n;
-		c[i].d = a[i].d;
-		}
-	else
-		{
-		radd( &a[i], &b[i], &c[i] ); /*c[i] = b[i] + a[i];*/
-		}
-	}
-}
-
-/* c = b - a.
- */
-void polsub( a, na, b, nb, c )
-fract a[], b[], c[];
-int na, nb;
-{
-int i, n;
-
-
-if( na > nb )
-	n = na;
-else
-	n = nb;
-
-if( n > MAXPOL )
-	n = MAXPOL;
-
-for( i=0; i<=n; i++ )
-	{
-	if( i > na )
-		{
-		c[i].n = b[i].n;
-		c[i].d = b[i].d;
-		}
-	else if( i > nb )
-		{
-		c[i].n = -a[i].n;
-		c[i].d = a[i].d;
-		}
-	else
-		{
-		rsub( &a[i], &b[i], &c[i] ); /*c[i] = b[i] - a[i];*/
-		}
-	}
-}
-
-
-
-/* c = b/a
- */
-int poldiv( a, na, b, nb, c )
-fract a[], b[], c[];
-int na, nb;
-{
-fract *ta, *tb, *tq;
-fract quot;
-fract temp;
-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 = (fract * )malloc( psize );
-polclr( ta, MAXPOL );
-polmov( a, na, ta );
-
-tb = (fract * )malloc( psize );
-polclr( tb, MAXPOL );
-polmov( b, nb, tb );
-
-tq = (fract * )malloc( psize );
-polclr( tq, MAXPOL );
-
-/* What to do if leading (constant) coefficient
- * of denominator is zero.
- */
-if( a[0].n == 0.0 )
-	{
-	for( i=0; i<=na; i++ )
-		{
-		if( ta[i].n != 0.0 )
-			goto nzero;
-		}
-	mtherr( "poldiv", SING );
-	goto done;
-
-nzero:
-/* Reduce the degree of the denominator. */
-	for( i=0; i<na; i++ )
-		{
-		ta[i].n = ta[i+1].n;
-		ta[i].d = ta[i+1].d;
-		}
-	ta[na].n = 0.0;
-	ta[na].d = 1.0;
-
-	if( b[0].n != 0.0 )
-		{
-/* Optional message:
-		printf( "poldiv singularity, divide quotient by x\n" );
-*/
-		sing += 1;
-		}
-	else
-		{
-/* Reduce degree of numerator. */
-		for( i=0; i<nb; i++ )
-			{
-			tb[i].n = tb[i+1].n;
-			tb[i].d = tb[i+1].d;
-			}
-		tb[nb].n = 0.0;
-		tb[nb].d = 1.0;
-		}
-/* Call self, using reduced polynomials. */
-	sing += poldiv( ta, na, tb, nb, c );
-	goto done;
-	}
-
-/* Long division algorithm.  ta[0] is nonzero.
- */
-for( i=0; i<=MAXPOL; i++ )
-	{
-	rdiv( &ta[0], &tb[i], &quot ); /*quot = tb[i]/ta[0];*/
-	for( j=0; j<=MAXPOL; j++ )
-		{
-		k = j + i;
-		if( k > MAXPOL )
-			break;
-
-		rmul( &ta[j], &quot, &temp ); /*tb[k] -= quot * ta[j];*/
-		rsub( &temp, &tb[k], &tb[k] );
-		}
-	tq[i].n = quot.n;
-	tq[i].d = quot.d;
-	}
-/* Send quotient to output array. */
-polmov( tq, MAXPOL, 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)).
- */
-
-void polsbt( a, na, b, nb, c )
-fract a[], b[], c[];
-int na, nb;
-{
-int i, j, k, n2;
-fract temp;
-fract *p;
-
-/* 0th degree term:
- */
-polclr( pt1, MAXPOL );
-pt1[0].n = b[0].n;
-pt1[0].d = b[0].d;
-
-polclr( pt2, MAXPOL );
-pt2[0].n = 1.0;
-pt2[0].d = 1.0;
-n2 = 0;
-p = &b[1];
-
-for( i=1; i<=nb; i++ )
-	{
-/* Form ith power of a. */
-	polmul( a, na, pt2, n2, pt2 );
-	n2 += na;
-/* Add the ith coefficient of b times the ith power of a. */
-	for( j=0; j<=n2; j++ )
-		{
-		if( j > MAXPOL )
-			break;
-		rmul( &pt2[j], p, &temp ); /*pt1[j] += b[i] * pt2[j];*/
-		radd( &temp, &pt1[j], &pt1[j] );
-		}
-	++p;
-	}
-
-k = n2 + nb;
-if( k > MAXPOL )
-	k = MAXPOL;
-for( i=0; i<=k; i++ )
-	{
-	c[i].n = pt1[i].n;
-	c[i].d = pt1[i].d;
-	}
-}
-
-
-
-
-/* Evaluate polynomial a(t) at t = x.
- */
-void poleva( a, na, x, s )
-fract a[];
-int na;
-fract *x;
-fract *s;
-{
-int i;
-fract temp;
-
-s->n = a[na].n;
-s->d = a[na].d;
-for( i=na-1; i>=0; i-- )
-	{
-	rmul( s, x, &temp ); /*s = s * x + a[i];*/
-	radd( &a[i], &temp, s );
-	}
-}
-