1 /* polyn.c
2 * polyr.c
3 * Arithmetic operations on polynomials
4 *
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
11 *
12 * polini( maxpol );
13 *
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
17 * by polini().
18 *
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,
23 *
24 * 2 na
25 * a(x) = a[0] + a[1] * x + a[2] * x + ... + a[na] * x .
26 *
27 *
28 *
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
36 *
37 *
38 * Division:
39 *
40 * i = poldiv( a, na, b, nb, c ); c = b / a, nc = MAXPOL
41 *
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.
45 *
46 *
47 * Change of variables:
48 * If a and b are polynomials, and t = a(x), then
49 * c(t) = b(a(x))
50 * is a polynomial found by substituting a(x) for t. The
51 * subroutine call for this is
52 *
53 * polsbt( a, na, b, nb, c );
54 *
55 *
56 * Notes:
57 * poldiv() is an integer routine; poleva() is double.
58 * Any of the arguments a, b, c may refer to the same array.
59 *
60 */
61
62 #include "mconf.h"
63 #include "cephes.h"
64
65 #include <stdio.h>
66 #include <stdlib.h>
67
68 /* Pointers to internal arrays. Note poldiv() allocates
69 * and deallocates some temporary arrays every time it is called.
70 */
71 static double *pt1 = 0;
72 static double *pt2 = 0;
73 static double *pt3 = 0;
74
75 /* Maximum degree of polynomial. */
76 int MAXPOL = 0;
77 extern int MAXPOL;
78
79 /* Number of bytes (chars) in maximum size polynomial. */
80 static int psize = 0;
81
82
83 /* Initialize max degree of polynomials
84 * and allocate temporary storage.
85 */
polini(maxdeg)86 void polini( maxdeg )
87 int maxdeg;
88 {
89
90 MAXPOL = maxdeg;
91 psize = (maxdeg + 1) * sizeof(double);
92
93 /* Release previously allocated memory, if any. */
94 if( pt3 )
95 free(pt3);
96 if( pt2 )
97 free(pt2);
98 if( pt1 )
99 free(pt1);
100
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 */
105
106 /* Report if failure */
107 if( (pt1 == NULL) || (pt2 == NULL) || (pt3 == NULL) )
108 {
109 mtherr( "polini", ERANGE );
110 exit(1);
111 }
112 }
113
114
115 /* Set a = 0.
116 */
polclr(a,n)117 void polclr( a, n )
118 register double *a;
119 int n;
120 {
121 int i;
122
123 if( n > MAXPOL )
124 n = MAXPOL;
125 for( i=0; i<=n; i++ )
126 *a++ = 0.0;
127 }
128
129
130
131 /* Set b = a.
132 */
polmov(a,na,b)133 void polmov( a, na, b )
134 register double *a, *b;
135 int na;
136 {
137 int i;
138
139 if( na > MAXPOL )
140 na = MAXPOL;
141
142 for( i=0; i<= na; i++ )
143 {
144 *b++ = *a++;
145 }
146 }
147
148
149 /* c = b * a.
150 */
polmul(a,na,b,nb,c)151 void polmul( a, na, b, nb, c )
152 double a[], b[], c[];
153 int na, nb;
154 {
155 int i, j, k, nc;
156 double x;
157
158 nc = na + nb;
159 polclr( pt3, MAXPOL );
160
161 for( i=0; i<=na; i++ )
162 {
163 x = a[i];
164 for( j=0; j<=nb; j++ )
165 {
166 k = i + j;
167 if( k > MAXPOL )
168 break;
169 pt3[k] += x * b[j];
170 }
171 }
172
173 if( nc > MAXPOL )
174 nc = MAXPOL;
175 for( i=0; i<=nc; i++ )
176 c[i] = pt3[i];
177 }
178
179
180
181
182 /* c = b + a.
183 */
poladd(a,na,b,nb,c)184 void poladd( a, na, b, nb, c )
185 double a[], b[], c[];
186 int na, nb;
187 {
188 int i, n;
189
190
191 if( na > nb )
192 n = na;
193 else
194 n = nb;
195
196 if( n > MAXPOL )
197 n = MAXPOL;
198
199 for( i=0; i<=n; i++ )
200 {
201 if( i > na )
202 c[i] = b[i];
203 else if( i > nb )
204 c[i] = a[i];
205 else
206 c[i] = b[i] + a[i];
207 }
208 }
209
210 /* c = b - a.
211 */
polsub(a,na,b,nb,c)212 void polsub( a, na, b, nb, c )
213 double a[], b[], c[];
214 int na, nb;
215 {
216 int i, n;
217
218
219 if( na > nb )
220 n = na;
221 else
222 n = nb;
223
224 if( n > MAXPOL )
225 n = MAXPOL;
226
227 for( i=0; i<=n; i++ )
228 {
229 if( i > na )
230 c[i] = b[i];
231 else if( i > nb )
232 c[i] = -a[i];
233 else
234 c[i] = b[i] - a[i];
235 }
236 }
237
238
239
240 /* c = b/a
241 */
poldiv(a,na,b,nb,c)242 int poldiv( a, na, b, nb, c )
243 double a[], b[], c[];
244 int na, nb;
245 {
246 double quot;
247 double *ta, *tb, *tq;
248 int i, j, k, sing;
249
250 sing = 0;
251
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.
255 */
256 ta = (double * )malloc( psize );
257 polclr( ta, MAXPOL );
258 polmov( a, na, ta );
259
260 tb = (double * )malloc( psize );
261 polclr( tb, MAXPOL );
262 polmov( b, nb, tb );
263
264 tq = (double * )malloc( psize );
265 polclr( tq, MAXPOL );
266
267 /* What to do if leading (constant) coefficient
268 * of denominator is zero.
269 */
270 if( a[0] == 0.0 )
271 {
272 for( i=0; i<=na; i++ )
273 {
274 if( ta[i] != 0.0 )
275 goto nzero;
276 }
277 mtherr( "poldiv", SING );
278 goto done;
279
280 nzero:
281 /* Reduce the degree of the denominator. */
282 for( i=0; i<na; i++ )
283 ta[i] = ta[i+1];
284 ta[na] = 0.0;
285
286 if( b[0] != 0.0 )
287 {
288 /* Optional message:
289 printf( "poldiv singularity, divide quotient by x\n" );
290 */
291 sing += 1;
292 }
293 else
294 {
295 /* Reduce degree of numerator. */
296 for( i=0; i<nb; i++ )
297 tb[i] = tb[i+1];
298 tb[nb] = 0.0;
299 }
300 /* Call self, using reduced polynomials. */
301 sing += poldiv( ta, na, tb, nb, c );
302 goto done;
303 }
304
305 /* Long division algorithm. ta[0] is nonzero.
306 */
307 for( i=0; i<=MAXPOL; i++ )
308 {
309 quot = tb[i]/ta[0];
310 for( j=0; j<=MAXPOL; j++ )
311 {
312 k = j + i;
313 if( k > MAXPOL )
314 break;
315 tb[k] -= quot * ta[j];
316 }
317 tq[i] = quot;
318 }
319 /* Send quotient to output array. */
320 polmov( tq, MAXPOL, c );
321
322 done:
323
324 /* Restore allocated memory. */
325 free(tq);
326 free(tb);
327 free(ta);
328 return( sing );
329 }
330
331
332
333
334 /* Change of variables
335 * Substitute a(y) for the variable x in b(x).
336 * x = a(y)
337 * c(x) = b(x) = b(a(y)).
338 */
339
polsbt(a,na,b,nb,c)340 void polsbt( a, na, b, nb, c )
341 double a[], b[], c[];
342 int na, nb;
343 {
344 int i, j, k, n2;
345 double x;
346
347 /* 0th degree term:
348 */
349 polclr( pt1, MAXPOL );
350 pt1[0] = b[0];
351
352 polclr( pt2, MAXPOL );
353 pt2[0] = 1.0;
354 n2 = 0;
355
356 for( i=1; i<=nb; i++ )
357 {
358 /* Form ith power of a. */
359 polmul( a, na, pt2, n2, pt2 );
360 n2 += na;
361 x = b[i];
362 /* Add the ith coefficient of b times the ith power of a. */
363 for( j=0; j<=n2; j++ )
364 {
365 if( j > MAXPOL )
366 break;
367 pt1[j] += x * pt2[j];
368 }
369 }
370
371 k = n2 + nb;
372 if( k > MAXPOL )
373 k = MAXPOL;
374 for( i=0; i<=k; i++ )
375 c[i] = pt1[i];
376 }
377
378
379
380
381 /* Evaluate polynomial a(t) at t = x.
382 */
poleva(a,na,x)383 double poleva( a, na, x )
384 double a[];
385 int na;
386 double x;
387 {
388 double s;
389 int i;
390
391 s = a[na];
392 for( i=na-1; i>=0; i-- )
393 {
394 s = s * x + a[i];
395 }
396 return(s);
397 }
398
399