1(* $Id: ComplexMath.Mod,v 1.5 1999/09/02 13:05:36 acken Exp $ *) 2MODULE oocComplexMath; 3 4 (* 5 ComplexMath - Mathematical functions for the type COMPLEX. 6 7 Copyright (C) 1995-1996 Michael Griebling 8 9 This module is free software; you can redistribute it and/or modify 10 it under the terms of the GNU Lesser General Public License as 11 published by the Free Software Foundation; either version 2 of the 12 License, or (at your option) any later version. 13 14 This module is distributed in the hope that it will be useful, 15 but WITHOUT ANY WARRANTY; without even the implied warranty of 16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 17 GNU Lesser General Public License for more details. 18 19 You should have received a copy of the GNU Lesser General Public 20 License along with this program; if not, write to the Free Software 21 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 22 23*) 24 25IMPORT m := oocRealMath; 26 27TYPE 28 COMPLEX * = POINTER TO COMPLEXDesc; 29 COMPLEXDesc = RECORD 30 r, i : REAL 31 END; 32 33CONST 34 ZERO=0.0; HALF=0.5; ONE=1.0; TWO=2.0; 35 36VAR 37 i-, one-, zero- : COMPLEX; 38 39PROCEDURE CMPLX * (r, i: REAL): COMPLEX; 40VAR c: COMPLEX; 41BEGIN 42 NEW(c); c.r:=r; c.i:=i; 43 RETURN c 44END CMPLX; 45 46(* 47 NOTE: This function provides the only way 48 of reliably assigning COMPLEX numbers. DO 49 NOT use ` a := b' where a, b are COMPLEX! 50 *) 51PROCEDURE Copy * (z: COMPLEX): COMPLEX; 52BEGIN 53 RETURN CMPLX(z.r, z.i) 54END Copy; 55 56PROCEDURE RealPart * (z: COMPLEX): REAL; 57BEGIN 58 RETURN z.r 59END RealPart; 60 61PROCEDURE ImagPart * (z: COMPLEX): REAL; 62BEGIN 63 RETURN z.i 64END ImagPart; 65 66PROCEDURE add * (z1, z2: COMPLEX): COMPLEX; 67BEGIN 68 RETURN CMPLX(z1.r+z2.r, z1.i+z2.i) 69END add; 70 71PROCEDURE sub * (z1, z2: COMPLEX): COMPLEX; 72BEGIN 73 RETURN CMPLX(z1.r-z2.r, z1.i-z2.i) 74END sub; 75 76PROCEDURE mul * (z1, z2: COMPLEX): COMPLEX; 77BEGIN 78 RETURN CMPLX(z1.r*z2.r-z1.i*z2.i, z1.r*z2.i+z1.i*z2.r) 79END mul; 80 81PROCEDURE div * (z1, z2: COMPLEX): COMPLEX; 82 VAR d, h: REAL; 83BEGIN 84 (* Note: this algorith avoids overflow by avoiding 85 multiplications and using divisions instead so that: 86 87 Re(z1/z2) = (z1.r*z2.r+z1.i*z2.i)/(z2.r^2+z2.i^2) 88 = (z1.r+z1.i*z2.i/z2.r)/(z2.r+z2.i^2/z2.r) 89 = (z1.r+h*z1.i)/(z2.r+h*z2.i) 90 Im(z1/z2) = (z1.i*z2.r-z1.r*z2.i)/(z2.r^2+z2.i^2) 91 = (z1.i-z1.r*z2.i/z2.r)/(z2.r+z2.i^2/z2.r) 92 = (z1.i-h*z1.r)/(z2.r+h*z2.i) 93 94 where h=z2.i/z2.r, provided z2.i<=z2.r and similarly 95 for z2.i>z2.r we have: 96 97 Re(z1/z2) = (h*z1.r+z1.i)/(h*z2.r+z2.i) 98 Im(z1/z2) = (h*z1.i-z1.r)/(h*z2.r+z2.i) 99 100 where h=z2.r/z2.i *) 101 102 (* we always guarantee h<=1 *) 103 IF ABS(z2.r)>ABS(z2.i) THEN 104 h:=z2.i/z2.r; d:=z2.r+h*z2.i; 105 RETURN CMPLX((z1.r+h*z1.i)/d, (z1.i-h*z1.r)/d) 106 ELSE 107 h:=z2.r/z2.i; d:=h*z2.r+z2.i; 108 RETURN CMPLX((h*z1.r+z1.i)/d, (h*z1.i-z1.r)/d) 109 END 110END div; 111 112PROCEDURE abs * (z: COMPLEX): REAL; 113 (* Returns the length of z *) 114VAR 115 r, i, h: REAL; 116BEGIN 117 (* Note: this algorithm avoids overflow by avoiding 118 multiplications and using divisions instead so that: 119 120 abs(z) = sqrt(z.r*z.r+z.i*z.i) 121 = sqrt(z.r^2*(1+(z.i/z.r)^2)) 122 = z.r*sqrt(1+(z.i/z.r)^2) 123 124 where z.i/z.r <= 1.0 by swapping z.r & z.i so that 125 for z.r>z.i we have z.r*sqrt(1+(z.i/z.r)^2) and 126 otherwise we have z.i*sqrt(1+(z.r/z.i)^2) *) 127 r:=ABS(z.r); i:=ABS(z.i); 128 IF i>r THEN h:=i; i:=r; r:=h END; (* guarantees i<=r *) 129 IF i=ZERO THEN RETURN r END; (* i=0, so sqrt(0+r^2)=r *) 130 h:=i/r; 131 RETURN r*m.sqrt(ONE+h*h) (* r*sqrt(1+(i/r)^2) *) 132END abs; 133 134PROCEDURE arg * (z: COMPLEX): REAL; 135 (* Returns the angle that z subtends to the positive real axis, in the range [-pi, pi] *) 136BEGIN 137 RETURN m.arctan2(z.i, z.r) 138END arg; 139 140PROCEDURE conj * (z: COMPLEX): COMPLEX; 141 (* Returns the complex conjugate of z *) 142BEGIN 143 RETURN CMPLX(z.r, -z.i) 144END conj; 145 146PROCEDURE power * (base: COMPLEX; exponent: REAL): COMPLEX; 147 (* Returns the value of the number base raised to the power exponent *) 148VAR c, s, r: REAL; 149BEGIN 150 m.sincos(arg(base)*exponent, s, c); r:=m.power(abs(base), exponent); 151 RETURN CMPLX(c*r, s*r) 152END power; 153 154PROCEDURE sqrt * (z: COMPLEX): COMPLEX; 155 (* Returns the principal square root of z, with arg in the range [-pi/2, pi/2] *) 156VAR u, v: REAL; 157BEGIN 158 (* Note: the following algorithm is more efficient since 159 it doesn't require a sincos or arctan evaluation: 160 161 Re(sqrt(z)) = sqrt((abs(z)+z.r)/2), Im(sqrt(z)) = +/-sqrt((abs(z)-z.r)/2) 162 = u = +/-v 163 164 where z.r >= 0 and z.i = 2*u*v and unknown sign is sign of z.i *) 165 166 (* initially force z.r >= 0 to calculate u, v *) 167 u:=m.sqrt((abs(z)+ABS(z.r))*HALF); 168 IF z.i#ZERO THEN v:=(HALF*z.i)/u ELSE v:=ZERO END; (* slight optimization *) 169 170 (* adjust u, v for the signs of z.r and z.i *) 171 IF z.r>=ZERO THEN RETURN CMPLX(u, v) (* no change *) 172 ELSIF z.i>=ZERO THEN RETURN CMPLX(v, u) (* z.r<0 so swap u, v *) 173 ELSE RETURN CMPLX(-v, -u) (* z.r<0, z.i<0 *) 174 END 175END sqrt; 176 177PROCEDURE exp * (z: COMPLEX): COMPLEX; 178 (* Returns the complex exponential of z *) 179VAR c, s, e: REAL; 180BEGIN 181 m.sincos(z.i, s, c); e:=m.exp(z.r); 182 RETURN CMPLX(e*c, e*s) 183END exp; 184 185PROCEDURE ln * (z: COMPLEX): COMPLEX; 186 (* Returns the principal value of the natural logarithm of z *) 187BEGIN 188 RETURN CMPLX(m.ln(abs(z)), arg(z)) 189END ln; 190 191PROCEDURE sin * (z: COMPLEX): COMPLEX; 192 (* Returns the sine of z *) 193 VAR s, c: REAL; 194BEGIN 195 m.sincos(z.r, s, c); 196 RETURN CMPLX(s*m.cosh(z.i), c*m.sinh(z.i)) 197END sin; 198 199PROCEDURE cos * (z: COMPLEX): COMPLEX; 200 (* Returns the cosine of z *) 201 VAR s, c: REAL; 202BEGIN 203 m.sincos(z.r, s, c); 204 RETURN CMPLX(c*m.cosh(z.i), -s*m.sinh(z.i)) 205END cos; 206 207PROCEDURE tan * (z: COMPLEX): COMPLEX; 208 (* Returns the tangent of z *) 209 VAR s, c, y, d: REAL; 210BEGIN 211 m.sincos(TWO*z.r, s, c); 212 y:=TWO*z.i; d:=c+m.cosh(y); 213 RETURN CMPLX(s/d, m.sinh(y)/d) 214END tan; 215 216PROCEDURE CalcAlphaBeta(z: COMPLEX; VAR a, b: REAL); 217 VAR x, x2, y, r, t: REAL; 218BEGIN x:=z.r+ONE; x:=x*x; y:=z.i*z.i; 219 x2:=z.r-ONE; x2:=x2*x2; 220 r:=m.sqrt(x+y); t:=m.sqrt(x2+y); 221 a:=HALF*(r+t); b:=HALF*(r-t); 222END CalcAlphaBeta; 223 224PROCEDURE arcsin * (z: COMPLEX): COMPLEX; 225 (* Returns the arcsine of z *) 226 VAR a, b: REAL; 227BEGIN 228 CalcAlphaBeta(z, a, b); 229 RETURN CMPLX(m.arcsin(b), m.ln(a+m.sqrt(a*a-1))) 230END arcsin; 231 232PROCEDURE arccos * (z: COMPLEX): COMPLEX; 233 (* Returns the arccosine of z *) 234 VAR a, b: REAL; 235BEGIN 236 CalcAlphaBeta(z, a, b); 237 RETURN CMPLX(m.arccos(b), -m.ln(a+m.sqrt(a*a-1))) 238END arccos; 239 240PROCEDURE arctan * (z: COMPLEX): COMPLEX; 241 (* Returns the arctangent of z *) 242 VAR x, y, yp, x2, y2: REAL; 243BEGIN 244 x:=TWO*z.r; y:=z.i+ONE; y:=y*y; 245 yp:=z.i-ONE; yp:=yp*yp; 246 x2:=z.r*z.r; y2:=z.i*z.i; 247 RETURN CMPLX(HALF*m.arctan(x/(ONE-x2-y2)), 0.25*m.ln((x2+y)/(x2+yp))) 248END arctan; 249 250PROCEDURE polarToComplex * (abs, arg: REAL): COMPLEX; 251 (* Returns the complex number with the specified polar coordinates *) 252BEGIN 253 RETURN CMPLX(abs*m.cos(arg), abs*m.sin(arg)) 254END polarToComplex; 255 256PROCEDURE scalarMult * (scalar: REAL; z: COMPLEX): COMPLEX; 257 (* Returns the scalar product of scalar with z *) 258BEGIN 259 RETURN CMPLX(z.r*scalar, z.i*scalar) 260END scalarMult; 261 262PROCEDURE IsCMathException * (): BOOLEAN; 263 (* Returns TRUE if the current coroutine is in the exceptional execution state 264 because of the ComplexMath exception; otherwise returns FALSE. 265 *) 266BEGIN 267 RETURN FALSE 268END IsCMathException; 269 270BEGIN 271 i:=CMPLX (ZERO, ONE); 272 one:=CMPLX (ONE, ZERO); 273 zero:=CMPLX (ZERO, ZERO) 274END oocComplexMath. 275