/Proj4/pj_zpoly1.c
http://github.com/route-me/route-me · C · 48 lines · 40 code · 2 blank · 6 comment · 4 complexity · b25165e8fe086229f75fa2d60e081792 MD5 · raw file
- /* 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;
- }