diff options
Diffstat (limited to 'libm/double/polyr.c')
| -rw-r--r-- | libm/double/polyr.c | 533 |
1 files changed, 533 insertions, 0 deletions
diff --git a/libm/double/polyr.c b/libm/double/polyr.c new file mode 100644 index 000000000..81ca817e3 --- /dev/null +++ b/libm/double/polyr.c @@ -0,0 +1,533 @@ + +/* 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 = tb[i]/ta[0];*/ + for( j=0; j<=MAXPOL; j++ ) + { + k = j + i; + if( k > MAXPOL ) + break; + + rmul( &ta[j], ", &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 ); + } +} + |