summaryrefslogtreecommitdiff
path: root/libm/float/polynf.c
diff options
context:
space:
mode:
Diffstat (limited to 'libm/float/polynf.c')
-rw-r--r--libm/float/polynf.c520
1 files changed, 0 insertions, 520 deletions
diff --git a/libm/float/polynf.c b/libm/float/polynf.c
deleted file mode 100644
index 48c6675d4..000000000
--- a/libm/float/polynf.c
+++ /dev/null
@@ -1,520 +0,0 @@
-/* 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);
-}
-