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