redoing dev style to be more test centric, from lessons learned with lisp-matrix.
[CommonLispStat.git] / lib / complex.c
blob4056b2e4092901afa003d7c51a27a894b25005ec
1 /* complex - Complex number functions */
2 /* Copyright (c) 1990, by Luke Tierney */
4 /* patched up and semi-ansified, (c) 2006, AJ Rossini, blindglobe@gmail.com */
6 #include "xmath.h"
7 #include "complex.h"
9 extern void xlfail(char *);
11 /*static */
12 double
13 phase(Complex c)
15 double phi;
17 if (c.real == 0.0) {
18 if (c.imag > 0.0) phi = PI / 2;
19 else if (c.imag == 0.0) phi = 0.0;
20 else phi = -PI / 2;
21 } else {
22 phi = atan(c.imag / c.real);
23 if (c.real < 0.0) {
24 if (c.imag > 0.0) {
25 phi += PI;
26 } else {
27 if (c.imag < 0.0) {
28 phi -= PI;
29 } else {
30 phi = PI;
35 return(phi);
38 double
39 modulus(Complex c)
41 return(sqrt(c.real * c.real + c.imag * c.imag));
44 Complex
45 cart2complex(double real, double imag)
47 Complex val;
48 val.real = real;
49 val.imag = imag;
50 return(val);
53 Complex
54 polar2complex(double mod, double phi)
56 Complex val;
57 double cs, sn;
59 if (phi == 0) {
60 cs = 1.0;
61 sn = 0.0;
62 } else {
63 if (phi == PI / 2) {
64 cs = 0.0;
65 sn = 1.0;
66 } else {
67 if (phi == PI) {
68 cs = -1.0;
69 sn = 0.0;
70 } else {
71 if (phi == -PI / 2) {
72 cs = 0.0;
73 sn = -1.0;
74 } else {
75 cs = cos(phi);
76 sn = sin(phi);
81 val.real = mod * cs;
82 val.imag = mod * sn;
83 return(val);
86 Complex
87 cadd(Complex c1,
88 Complex c2)
90 return(cart2complex(c1.real + c2.real, c1.imag + c2.imag));
93 Complex
94 csub(Complex c1,
95 Complex c2)
97 return(cart2complex(c1.real - c2.real, c1.imag - c2.imag));
100 Complex
101 cmul(Complex c1, Complex c2)
103 double m1, m2, p1, p2;
105 m1 = modulus(c1);
106 p1 = phase(c1);
107 m2 = modulus(c2);
108 p2 = phase(c2);
109 return(polar2complex(m1 * m2, p1 + p2));
112 static void
113 checkfzero(double x)
115 if (x == 0.0) {
116 xlfail("division by zero");
120 Complex
121 cdiv(Complex c1, Complex c2)
123 double m1, m2, p1, p2;
125 m1 = modulus(c1);
126 p1 = phase(c1);
127 m2 = modulus(c2);
128 p2 = phase(c2);
129 checkfzero(m2);
130 return(polar2complex(m1 / m2, p1 - p2));