1 2 /++ 3 Complex math 4 5 Copyright: Ilia Ki; 2010, Lars T. Kyllingstad (original Phobos code) 6 Authors: Ilia Ki, Lars Tandle Kyllingstad, Don Clugston 7 +/ 8 module mir.complex.math; 9 10 public import mir.complex; 11 12 /++ 13 Params: z = A complex number. 14 Returns: The square root of `z`. 15 +/ 16 Complex!T sqrt(T)(Complex!T z) @safe pure nothrow @nogc 17 { 18 import mir.math.common: fabs, fmin, fmax, sqrt; 19 20 if (z == 0) 21 return typeof(return)(0, 0); 22 auto x = fabs(z.re); 23 auto y = fabs(z.im); 24 auto n = fmin(x, y); 25 auto m = fmax(x, y); 26 auto r = n / m; 27 auto w = sqrt(m) * sqrt(0.5f * ((x >= y ? 1 : r) + sqrt(1 + r * r))); 28 auto s = typeof(return)(w, z.im / (w + w)); 29 if (z.re < 0) 30 { 31 s = typeof(return)(s.im, s.re); 32 if (z.im < 0) 33 s = -s; 34 } 35 return s; 36 } 37 38 /// 39 @safe pure nothrow unittest 40 { 41 assert(sqrt(complex(0.0)) == 0.0); 42 assert(sqrt(complex(1.0, 0)) == 1.0); 43 assert(sqrt(complex(-1.0, 0)) == complex(0, 1.0)); 44 assert(sqrt(complex(-8.0, -6.0)) == complex(1.0, -3.0)); 45 } 46 47 @safe pure nothrow unittest 48 { 49 assert(complex(1.0, 1.0).sqrt.approxEqual(complex(1.098684113467809966, 0.455089860562227341))); 50 assert(complex(0.5, 2.0).sqrt.approxEqual(complex(1.131713924277869410, 0.883615530875513265))); 51 } 52 53 /** 54 * Calculate the natural logarithm of x. 55 * The branch cut is along the negative axis. 56 * Params: 57 * x = A complex number 58 * Returns: 59 * The complex natural logarithm of `x` 60 * 61 * $(TABLE_SV 62 * $(TR $(TH x) $(TH log(x))) 63 * $(TR $(TD (-0, +0)) $(TD (-$(INFIN), $(PI)))) 64 * $(TR $(TD (+0, +0)) $(TD (-$(INFIN), +0))) 65 * $(TR $(TD (any, +$(INFIN))) $(TD (+$(INFIN), $(PI)/2))) 66 * $(TR $(TD (any, $(NAN))) $(TD ($(NAN), $(NAN)))) 67 * $(TR $(TD (-$(INFIN), any)) $(TD (+$(INFIN), $(PI)))) 68 * $(TR $(TD (+$(INFIN), any)) $(TD (+$(INFIN), +0))) 69 * $(TR $(TD (-$(INFIN), +$(INFIN))) $(TD (+$(INFIN), 3$(PI)/4))) 70 * $(TR $(TD (+$(INFIN), +$(INFIN))) $(TD (+$(INFIN), $(PI)/4))) 71 * $(TR $(TD ($(PLUSMN)$(INFIN), $(NAN))) $(TD (+$(INFIN), $(NAN)))) 72 * $(TR $(TD ($(NAN), any)) $(TD ($(NAN), $(NAN)))) 73 * $(TR $(TD ($(NAN), +$(INFIN))) $(TD (+$(INFIN), $(NAN)))) 74 * $(TR $(TD ($(NAN), $(NAN))) $(TD ($(NAN), $(NAN)))) 75 * ) 76 */ 77 Complex!T log(T)(Complex!T x) @safe pure nothrow @nogc 78 { 79 import mir.math.constant: PI, PI_4, PI_2; 80 import mir.math.common: log, fabs, copysign; 81 alias isNaN = x => x != x; 82 alias isInfinity = x => x.fabs == T.infinity; 83 84 // Handle special cases explicitly here for better accuracy. 85 // The order here is important, so that the correct path is chosen. 86 if (isNaN(x.re)) 87 { 88 if (isInfinity(x.im)) 89 return Complex!T(T.infinity, T.nan); 90 else 91 return Complex!T(T.nan, T.nan); 92 } 93 if (isInfinity(x.re)) 94 { 95 if (isNaN(x.im)) 96 return Complex!T(T.infinity, T.nan); 97 else if (isInfinity(x.im)) 98 { 99 if (copysign(1, x.re) < 0) 100 return Complex!T(T.infinity, copysign(3.0 * PI_4, x.im)); 101 else 102 return Complex!T(T.infinity, copysign(PI_4, x.im)); 103 } 104 else 105 { 106 if (copysign(1, x.re) < 0) 107 return Complex!T(T.infinity, copysign(PI, x.im)); 108 else 109 return Complex!T(T.infinity, copysign(0.0, x.im)); 110 } 111 } 112 if (isNaN(x.im)) 113 return Complex!T(T.nan, T.nan); 114 if (isInfinity(x.im)) 115 return Complex!T(T.infinity, copysign(PI_2, x.im)); 116 if (x.re == 0.0 && x.im == 0.0) 117 { 118 if (copysign(1, x.re) < 0) 119 return Complex!T(-T.infinity, copysign(PI, x.im)); 120 else 121 return Complex!T(-T.infinity, copysign(0.0, x.im)); 122 } 123 124 return Complex!T(log(cabs(x)), arg(x)); 125 } 126 127 /// 128 @safe pure nothrow @nogc version(mir_core_test) unittest 129 { 130 import mir.math.common: sqrt; 131 import mir.math.constant: PI; 132 import mir.math.common: approxEqual; 133 134 auto a = complex(2.0, 1.0); 135 assert(log(conj(a)) == conj(log(a))); 136 137 assert(log(complex(-1.0L, 0.0L)) == complex(0.0L, PI)); 138 assert(log(complex(-1.0L, -0.0L)) == complex(0.0L, -PI)); 139 } 140 141 @safe pure nothrow @nogc version(mir_core_test) unittest 142 { 143 import mir.math.common: fabs; 144 import mir.math.constant: PI, PI_2, PI_4; 145 alias isNaN = x => x != x; 146 alias isInfinity = x => x.fabs == x.infinity; 147 148 auto a = log(complex(-0.0L, 0.0L)); 149 assert(a == complex(-real.infinity, PI)); 150 auto b = log(complex(0.0L, 0.0L)); 151 assert(b == complex(-real.infinity, +0.0L)); 152 auto c = log(complex(1.0L, real.infinity)); 153 assert(c == complex(real.infinity, PI_2)); 154 auto d = log(complex(1.0L, real.nan)); 155 assert(isNaN(d.re) && isNaN(d.im)); 156 157 auto e = log(complex(-real.infinity, 1.0L)); 158 assert(e == complex(real.infinity, PI)); 159 auto f = log(complex(real.infinity, 1.0L)); 160 assert(f == complex(real.infinity, 0.0L)); 161 auto g = log(complex(-real.infinity, real.infinity)); 162 assert(g == complex(real.infinity, 3.0 * PI_4)); 163 auto h = log(complex(real.infinity, real.infinity)); 164 assert(h == complex(real.infinity, PI_4)); 165 auto i = log(complex(real.infinity, real.nan)); 166 assert(isInfinity(i.re) && isNaN(i.im)); 167 168 auto j = log(complex(real.nan, 1.0L)); 169 assert(isNaN(j.re) && isNaN(j.im)); 170 auto k = log(complex(real.nan, real.infinity)); 171 assert(isInfinity(k.re) && isNaN(k.im)); 172 auto l = log(complex(real.nan, real.nan)); 173 assert(isNaN(l.re) && isNaN(l.im)); 174 } 175 176 @safe pure nothrow @nogc version(mir_core_test) unittest 177 { 178 import mir.math.constant: PI; 179 180 auto a = log(fromPolar(1.0, PI / 6.0)); 181 assert(approxEqual(a, complex(0.0L, 0.523598775598298873077L), 0.0, 1e-15)); 182 183 auto b = log(fromPolar(1.0, PI / 3.0)); 184 assert(approxEqual(b, complex(0.0L, 1.04719755119659774615L), 0.0, 1e-15)); 185 186 auto c = log(fromPolar(1.0, PI / 2.0)); 187 assert(approxEqual(c, complex(0.0L, 1.57079632679489661923L), 0.0, 1e-15)); 188 189 auto d = log(fromPolar(1.0, 2.0 * PI / 3.0)); 190 assert(approxEqual(d, complex(0.0L, 2.09439510239319549230L), 0.0, 1e-15)); 191 192 auto e = log(fromPolar(1.0, 5.0 * PI / 6.0)); 193 assert(approxEqual(e, complex(0.0L, 2.61799387799149436538L), 0.0, 1e-15)); 194 195 auto f = log(complex(-1.0L, 0.0L)); 196 assert(approxEqual(f, complex(0.0L, PI), 0.0, 1e-15)); 197 } 198 199 /++ 200 Calculates e$(SUPERSCRIPT x). 201 Params: 202 x = A complex number 203 Returns: 204 The complex base e exponential of `x` 205 $(TABLE_SV 206 $(TR $(TH x) $(TH exp(x))) 207 $(TR $(TD ($(PLUSMN)0, +0)) $(TD (1, +0))) 208 $(TR $(TD (any, +$(INFIN))) $(TD ($(NAN), $(NAN)))) 209 $(TR $(TD (any, $(NAN)) $(TD ($(NAN), $(NAN))))) 210 $(TR $(TD (+$(INFIN), +0)) $(TD (+$(INFIN), +0))) 211 $(TR $(TD (-$(INFIN), any)) $(TD ($(PLUSMN)0, cis(x.im)))) 212 $(TR $(TD (+$(INFIN), any)) $(TD ($(PLUSMN)$(INFIN), cis(x.im)))) 213 $(TR $(TD (-$(INFIN), +$(INFIN))) $(TD ($(PLUSMN)0, $(PLUSMN)0))) 214 $(TR $(TD (+$(INFIN), +$(INFIN))) $(TD ($(PLUSMN)$(INFIN), $(NAN)))) 215 $(TR $(TD (-$(INFIN), $(NAN))) $(TD ($(PLUSMN)0, $(PLUSMN)0))) 216 $(TR $(TD (+$(INFIN), $(NAN))) $(TD ($(PLUSMN)$(INFIN), $(NAN)))) 217 $(TR $(TD ($(NAN), +0)) $(TD ($(NAN), +0))) 218 $(TR $(TD ($(NAN), any)) $(TD ($(NAN), $(NAN)))) 219 $(TR $(TD ($(NAN), $(NAN))) $(TD ($(NAN), $(NAN)))) 220 ) 221 +/ 222 Complex!T exp(T)(Complex!T x) @trusted pure nothrow @nogc // TODO: @safe 223 { 224 import mir.math.common: exp, fabs, copysign; 225 alias isNaN = x => x != x; 226 alias isInfinity = x => x.fabs == T.infinity; 227 228 // Handle special cases explicitly here, as fromPolar will otherwise 229 // cause them to return Complex!T(NaN, NaN), or with the wrong sign. 230 if (isInfinity(x.re)) 231 { 232 if (isNaN(x.im)) 233 { 234 if (copysign(1, x.re) < 0) 235 return Complex!T(0, copysign(0, x.im)); 236 else 237 return x; 238 } 239 if (isInfinity(x.im)) 240 { 241 if (copysign(1, x.re) < 0) 242 return Complex!T(0, copysign(0, x.im)); 243 else 244 return Complex!T(T.infinity, -T.nan); 245 } 246 if (x.im == 0) 247 { 248 if (copysign(1, x.re) < 0) 249 return Complex!T(0); 250 else 251 return Complex!T(T.infinity); 252 } 253 } 254 if (isNaN(x.re)) 255 { 256 if (isNaN(x.im) || isInfinity(x.im)) 257 return Complex!T(T.nan, T.nan); 258 if (x.im == 0) 259 return x; 260 } 261 if (x.re == 0) 262 { 263 if (isNaN(x.im) || isInfinity(x.im)) 264 return Complex!T(T.nan, T.nan); 265 if (x.im == 0) 266 return Complex!T(1, 0); 267 } 268 269 return fromPolar!T(exp(x.re), x.im); 270 } 271 272 /// 273 @safe pure nothrow @nogc version(mir_core_test) unittest 274 { 275 import mir.math.constant: PI; 276 277 assert(exp(complex(0.0, 0.0)) == complex(1.0, 0.0)); 278 279 auto a = complex(2.0, 1.0); 280 assert(exp(conj(a)) == conj(exp(a))); 281 282 auto b = exp(complex(0.0, 1.0) * double(PI)); 283 assert(approxEqual(b, complex(-1.0), 0.0, 1e-15)); 284 } 285 286 @safe pure nothrow @nogc version(mir_core_test) unittest 287 { 288 import mir.math.common: fabs; 289 290 alias isNaN = x => x != x; 291 alias isInfinity = x => x.fabs == x.infinity; 292 293 auto a = exp(complex(0.0, double.infinity)); 294 assert(isNaN(a.re) && isNaN(a.im)); 295 auto b = exp(complex(0.0, double.infinity)); 296 assert(isNaN(b.re) && isNaN(b.im)); 297 auto c = exp(complex(0.0, double.nan)); 298 assert(isNaN(c.re) && isNaN(c.im)); 299 300 auto d = exp(complex(+double.infinity, 0.0)); 301 assert(d == complex(double.infinity, 0.0)); 302 auto e = exp(complex(-double.infinity, 0.0)); 303 assert(e == complex(0.0)); 304 auto f = exp(complex(-double.infinity, 1.0)); 305 assert(f == complex(0.0)); 306 auto g = exp(complex(+double.infinity, 1.0)); 307 assert(g == complex(double.infinity, double.infinity)); 308 auto h = exp(complex(-double.infinity, +double.infinity)); 309 assert(h == complex(0.0)); 310 auto i = exp(complex(+double.infinity, +double.infinity)); 311 assert(isInfinity(i.re) && isNaN(i.im)); 312 auto j = exp(complex(-double.infinity, double.nan)); 313 assert(j == complex(0.0)); 314 auto k = exp(complex(+double.infinity, double.nan)); 315 assert(isInfinity(k.re) && isNaN(k.im)); 316 317 auto l = exp(complex(double.nan, 0)); 318 assert(isNaN(l.re) && l.im == 0.0); 319 auto m = exp(complex(double.nan, 1)); 320 assert(isNaN(m.re) && isNaN(m.im)); 321 auto n = exp(complex(double.nan, double.nan)); 322 assert(isNaN(n.re) && isNaN(n.im)); 323 } 324 325 @safe pure nothrow @nogc version(mir_core_test) unittest 326 { 327 import mir.math.constant : PI; 328 329 auto a = exp(complex(0.0, -PI)); 330 assert(approxEqual(a, complex(-1.0L), 0.0, 1e-15)); 331 332 auto b = exp(complex(0.0, -2.0 * PI / 3.0)); 333 assert(approxEqual(b, complex(-0.5L, -0.866025403784438646763L))); 334 335 auto c = exp(complex(0.0, PI / 3.0)); 336 assert(approxEqual(c, complex(0.5L, 0.866025403784438646763L))); 337 338 auto d = exp(complex(0.0, 2.0 * PI / 3.0)); 339 assert(approxEqual(d, complex(-0.5L, 0.866025403784438646763L))); 340 341 auto e = exp(complex(0.0, PI)); 342 assert(approxEqual(e, complex(-1.0L), 0.0, 1e-15)); 343 } 344 345 /++ 346 Computes whether two values are approximately equal, admitting a maximum 347 relative difference, and a maximum absolute difference. 348 Params: 349 lhs = First item to compare. 350 rhs = Second item to compare. 351 maxRelDiff = Maximum allowable difference relative to `rhs`. Defaults to `0.5 ^^ 20`. 352 maxAbsDiff = Maximum absolute difference. Defaults to `0.5 ^^ 20`. 353 354 Returns: 355 `true` if the two items are equal or approximately equal under either criterium. 356 +/ 357 bool approxEqual(T)(Complex!T lhs, Complex!T rhs, const T maxRelDiff = 0x1p-20f, const T maxAbsDiff = 0x1p-20f) 358 { 359 import mir.math.common: approxEqual; 360 return approxEqual(lhs.re, rhs.re, maxRelDiff, maxAbsDiff) 361 && approxEqual(lhs.im, rhs.im, maxRelDiff, maxAbsDiff); 362 } 363 364 /// Complex types works as `approxEqual(l.re, r.re) && approxEqual(l.im, r.im)` 365 @safe pure nothrow @nogc version(mir_core_test) unittest 366 { 367 assert(approxEqual(complex(1.0, 1), complex(1.0000001, 1), 1.0000001)); 368 assert(!approxEqual(complex(100000.0, 0), complex(100001.0, 0))); 369 }