3
* Arithmetic operations on polynomials
5
* In the following descriptions a, b, c are polynomials of degree
6
* na, nb, nc respectively. The degree of a polynomial cannot
7
* exceed a run-time value MAXPOL. An operation that attempts
8
* to use or generate a polynomial of higher degree may produce a
9
* result that suffers truncation at degree MAXPOL. The value of
10
* MAXPOL is set by calling the function
14
* where maxpol is the desired maximum degree. This must be
15
* done prior to calling any of the other functions in this module.
16
* Memory for internal temporary polynomial storage is allocated
19
* Each polynomial is represented by an array containing its
20
* coefficients, together with a separately declared integer equal
21
* to the degree of the polynomial. The coefficients appear in
22
* ascending order; that is,
25
* a(x) = a[0] + a[1] * x + a[2] * x + ... + a[na] * x .
29
* sum = poleva( a, na, x ); Evaluate polynomial a(t) at t = x.
30
* polprt( a, na, D ); Print the coefficients of a to D digits.
31
* polclr( a, na ); Set a identically equal to zero, up to a[na].
32
* polmov( a, na, b ); Set b = a.
33
* poladd( a, na, b, nb, c ); c = b + a, nc = max(na,nb)
34
* polsub( a, na, b, nb, c ); c = b - a, nc = max(na,nb)
35
* polmul( a, na, b, nb, c ); c = b * a, nc = na+nb
40
* i = poldiv( a, na, b, nb, c ); c = b / a, nc = MAXPOL
42
* returns i = the degree of the first nonzero coefficient of a.
43
* The computed quotient c must be divided by x^i. An error message
44
* is printed if a is identically zero.
47
* Change of variables:
48
* If a and b are polynomials, and t = a(x), then
50
* is a polynomial found by substituting a(x) for t. The
51
* subroutine call for this is
53
* polsbt( a, na, b, nb, c );
57
* poldiv() is an integer routine; poleva() is double.
58
* Any of the arguments a, b, c may refer to the same array.
68
/* Pointers to internal arrays. Note poldiv() allocates
69
* and deallocates some temporary arrays every time it is called.
71
static double *pt1 = 0;
72
static double *pt2 = 0;
73
static double *pt3 = 0;
75
/* Maximum degree of polynomial. */
79
/* Number of bytes (chars) in maximum size polynomial. */
83
/* Initialize max degree of polynomials
84
* and allocate temporary storage.
91
psize = (maxdeg + 1) * sizeof(double);
93
/* Release previously allocated memory, if any. */
101
/* Allocate new arrays */
102
pt1 = (double * )malloc(psize); /* used by polsbt */
103
pt2 = (double * )malloc(psize); /* used by polsbt */
104
pt3 = (double * )malloc(psize); /* used by polmul */
106
/* Report if failure */
107
if( (pt1 == NULL) || (pt2 == NULL) || (pt3 == NULL) )
109
mtherr( "polini", ERANGE );
125
for( i=0; i<=n; i++ )
133
void polmov( a, na, b )
134
register double *a, *b;
142
for( i=0; i<= na; i++ )
151
void polmul( a, na, b, nb, c )
152
double a[], b[], c[];
159
polclr( pt3, MAXPOL );
161
for( i=0; i<=na; i++ )
164
for( j=0; j<=nb; j++ )
175
for( i=0; i<=nc; i++ )
184
void poladd( a, na, b, nb, c )
185
double a[], b[], c[];
199
for( i=0; i<=n; i++ )
212
void polsub( a, na, b, nb, c )
213
double a[], b[], c[];
227
for( i=0; i<=n; i++ )
242
int poldiv( a, na, b, nb, c )
243
double a[], b[], c[];
247
double *ta, *tb, *tq;
252
/* Allocate temporary arrays. This would be quicker
253
* if done automatically on the stack, but stack space
254
* may be hard to obtain on a small computer.
256
ta = (double * )malloc( psize );
257
polclr( ta, MAXPOL );
260
tb = (double * )malloc( psize );
261
polclr( tb, MAXPOL );
264
tq = (double * )malloc( psize );
265
polclr( tq, MAXPOL );
267
/* What to do if leading (constant) coefficient
268
* of denominator is zero.
272
for( i=0; i<=na; i++ )
277
mtherr( "poldiv", SING );
281
/* Reduce the degree of the denominator. */
282
for( i=0; i<na; i++ )
289
printf( "poldiv singularity, divide quotient by x\n" );
295
/* Reduce degree of numerator. */
296
for( i=0; i<nb; i++ )
300
/* Call self, using reduced polynomials. */
301
sing += poldiv( ta, na, tb, nb, c );
305
/* Long division algorithm. ta[0] is nonzero.
307
for( i=0; i<=MAXPOL; i++ )
310
for( j=0; j<=MAXPOL; j++ )
315
tb[k] -= quot * ta[j];
319
/* Send quotient to output array. */
320
polmov( tq, MAXPOL, c );
324
/* Restore allocated memory. */
334
/* Change of variables
335
* Substitute a(y) for the variable x in b(x).
337
* c(x) = b(x) = b(a(y)).
340
void polsbt( a, na, b, nb, c )
341
double a[], b[], c[];
349
polclr( pt1, MAXPOL );
352
polclr( pt2, MAXPOL );
356
for( i=1; i<=nb; i++ )
358
/* Form ith power of a. */
359
polmul( a, na, pt2, n2, pt2 );
362
/* Add the ith coefficient of b times the ith power of a. */
363
for( j=0; j<=n2; j++ )
367
pt1[j] += x * pt2[j];
374
for( i=0; i<=k; i++ )
381
/* Evaluate polynomial a(t) at t = x.
383
double poleva( a, na, x )
392
for( i=na-1; i>=0; i-- )
added
removed
cephes/polyn.c
