From 7ce331c01ce6eb7b3f5c715a38a24359da9c6ee2 Mon Sep 17 00:00:00 2001 From: Eric Andersen Date: Thu, 22 Nov 2001 14:04:29 +0000 Subject: 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 --- libm/double/polyn.c | 471 ---------------------------------------------------- 1 file changed, 471 deletions(-) delete mode 100644 libm/double/polyn.c (limited to 'libm/double/polyn.c') diff --git a/libm/double/polyn.c b/libm/double/polyn.c deleted file mode 100644 index 2927e77f0..000000000 --- a/libm/double/polyn.c +++ /dev/null @@ -1,471 +0,0 @@ -/* polyn.c - * polyr.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 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 . - * - * - * - * sum = poleva( a, na, x ); 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 -#include -#if ANSIPROT -void exit (int); -extern void * malloc ( long ); -extern void free ( void * ); -void polclr ( double *, int ); -void polmov ( double *, int, double * ); -void polmul ( double *, int, double *, int, double * ); -int poldiv ( double *, int, double *, int, double * ); -#else -void exit(); -void * malloc(); -void free (); -void polclr(), polmov(), poldiv(), polmul(); -#endif -#ifndef NULL -#define NULL 0 -#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 double *pt1 = 0; -static double *pt2 = 0; -static double *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(double); - -/* Release previously allocated memory, if any. */ -if( pt3 ) - free(pt3); -if( pt2 ) - free(pt2); -if( pt1 ) - free(pt1); - -/* Allocate new arrays */ -pt1 = (double * )malloc(psize); /* used by polsbt */ -pt2 = (double * )malloc(psize); /* used by polsbt */ -pt3 = (double * )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 ) -double 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] ); - } -printf( "\n" ); -} - - - -/* Set a = 0. - */ -void polclr( a, n ) -register double *a; -int n; -{ -int i; - -if( n > MAXPOL ) - n = MAXPOL; -for( i=0; i<=n; i++ ) - *a++ = 0.0; -} - - - -/* Set b = a. - */ -void polmov( a, na, b ) -register double *a, *b; -int na; -{ -int i; - -if( na > MAXPOL ) - na = MAXPOL; - -for( i=0; i<= na; i++ ) - { - *b++ = *a++; - } -} - - -/* c = b * a. - */ -void polmul( a, na, b, nb, c ) -double a[], b[], c[]; -int na, nb; -{ -int i, j, k, nc; -double x; - -nc = na + nb; -polclr( pt3, MAXPOL ); - -for( i=0; i<=na; i++ ) - { - x = a[i]; - for( j=0; j<=nb; j++ ) - { - k = i + j; - if( k > MAXPOL ) - break; - pt3[k] += x * b[j]; - } - } - -if( nc > MAXPOL ) - nc = MAXPOL; -for( i=0; i<=nc; i++ ) - c[i] = pt3[i]; -} - - - - -/* c = b + a. - */ -void poladd( a, na, b, nb, c ) -double 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] = b[i]; - else if( i > nb ) - c[i] = a[i]; - else - c[i] = b[i] + a[i]; - } -} - -/* c = b - a. - */ -void polsub( a, na, b, nb, c ) -double 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] = b[i]; - else if( i > nb ) - c[i] = -a[i]; - else - c[i] = b[i] - a[i]; - } -} - - - -/* c = b/a - */ -int poldiv( a, na, b, nb, c ) -double a[], b[], c[]; -int na, nb; -{ -double quot; -double *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 = (double * )malloc( psize ); -polclr( ta, MAXPOL ); -polmov( a, na, ta ); - -tb = (double * )malloc( psize ); -polclr( tb, MAXPOL ); -polmov( b, nb, tb ); - -tq = (double * )malloc( psize ); -polclr( tq, MAXPOL ); - -/* 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( "poldiv", SING ); - goto done; - -nzero: -/* Reduce the degree of the denominator. */ - for( i=0; i MAXPOL ) - break; - tb[k] -= quot * ta[j]; - } - tq[i] = quot; - } -/* 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 ) -double a[], b[], c[]; -int na, nb; -{ -int i, j, k, n2; -double x; - -/* 0th degree term: - */ -polclr( pt1, MAXPOL ); -pt1[0] = b[0]; - -polclr( pt2, MAXPOL ); -pt2[0] = 1.0; -n2 = 0; - -for( i=1; i<=nb; i++ ) - { -/* Form ith power of a. */ - polmul( 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 > MAXPOL ) - break; - pt1[j] += x * pt2[j]; - } - } - -k = n2 + nb; -if( k > MAXPOL ) - k = MAXPOL; -for( i=0; i<=k; i++ ) - c[i] = pt1[i]; -} - - - - -/* Evaluate polynomial a(t) at t = x. - */ -double poleva( a, na, x ) -double a[]; -int na; -double x; -{ -double s; -int i; - -s = a[na]; -for( i=na-1; i>=0; i-- ) - { - s = s * x + a[i]; - } -return(s); -} - -- cgit v1.2.3