/Proj4/pj_zpoly1.c

 1/* evaluate complex polynomial */
 2#ifndef lint
 3static const char SCCSID[]="@(#)pj_zpoly1.c	4.3	93/06/12	GIE	REL";
 4#endif
 5#include "projects.h"
 6/* note: coefficients are always from C_1 to C_n
 7**	i.e. C_0 == (0., 0)
 8**	n should always be >= 1 though no checks are made
 9*/
10	COMPLEX
11pj_zpoly1(COMPLEX z, COMPLEX *C, int n) {
12	COMPLEX a;
13	double t;
14
15	a = *(C += n);
16	while (n-- > 0) {
17		a.r = (--C)->r + z.r * (t = a.r) - z.i * a.i;
18		a.i = C->i + z.r * a.i + z.i * t;
19	}
20	a.r = z.r * (t = a.r) - z.i * a.i;
21	a.i = z.r * a.i + z.i * t;
22	return a;
23}
24/* evaluate complex polynomial and derivative */
25	COMPLEX
26pj_zpolyd1(COMPLEX z, COMPLEX *C, int n, COMPLEX *der) {
27	COMPLEX a, b;
28	double t;
29	int first = 1;
30
31	b = a = *(C += n);
32	while (n-- > 0) {
33		if (first) {
34			first = 0;
35		} else {
36			b.r = a.r + z.r * (t = b.r) - z.i * b.i;
37			b.i = a.i + z.r * b.i + z.i * t;
38		}
39		a.r = (--C)->r + z.r * (t = a.r) - z.i * a.i;
40		a.i = C->i + z.r * a.i + z.i * t;
41	}
42	b.r = a.r + z.r * (t = b.r) - z.i * b.i;
43	b.i = a.i + z.r * b.i + z.i * t;
44	a.r = z.r * (t = a.r) - z.i * a.i;
45	a.i = z.r * a.i + z.i * t;
46	*der = b;
47	return a;
48}