/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

  1. /* evaluate complex polynomial */
  2. #ifndef lint
  3. static 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
  11. pj_zpoly1(COMPLEX z, COMPLEX *C, int n) {
  12. COMPLEX a;
  13. double t;
  14. a = *(C += n);
  15. while (n-- > 0) {
  16. a.r = (--C)->r + z.r * (t = a.r) - z.i * a.i;
  17. a.i = C->i + z.r * a.i + z.i * t;
  18. }
  19. a.r = z.r * (t = a.r) - z.i * a.i;
  20. a.i = z.r * a.i + z.i * t;
  21. return a;
  22. }
  23. /* evaluate complex polynomial and derivative */
  24. COMPLEX
  25. pj_zpolyd1(COMPLEX z, COMPLEX *C, int n, COMPLEX *der) {
  26. COMPLEX a, b;
  27. double t;
  28. int first = 1;
  29. b = a = *(C += n);
  30. while (n-- > 0) {
  31. if (first) {
  32. first = 0;
  33. } else {
  34. b.r = a.r + z.r * (t = b.r) - z.i * b.i;
  35. b.i = a.i + z.r * b.i + z.i * t;
  36. }
  37. a.r = (--C)->r + z.r * (t = a.r) - z.i * a.i;
  38. a.i = C->i + z.r * a.i + z.i * t;
  39. }
  40. b.r = a.r + z.r * (t = b.r) - z.i * b.i;
  41. b.i = a.i + z.r * b.i + z.i * t;
  42. a.r = z.r * (t = a.r) - z.i * a.i;
  43. a.i = z.r * a.i + z.i * t;
  44. *der = b;
  45. return a;
  46. }