1 /++ 2 Complex numbers 3 4 Copyright: Ilia Ki; 2010, Lars T. Kyllingstad (original Phobos code) 5 Authors: Ilia Ki, Lars Tandle Kyllingstad, Don Clugston 6 +/ 7 module mir.complex; 8 9 import mir.math.common: optmath; 10 11 private alias CommonType(A, B) = typeof(A.init + B.init); 12 13 @optmath: 14 15 /++ 16 Generic complex number type 17 +/ 18 struct Complex(T) 19 if (is(T == float) || is(T == double) || is(T == real)) 20 { 21 import mir.internal.utility: isComplex; 22 import std.traits: isNumeric; 23 24 @optmath: 25 26 /++ 27 Real part. Default value is zero. 28 +/ 29 T re = 0; 30 /++ 31 Imaginary part. Default value is zero. 32 +/ 33 T im = 0; 34 35 /// 36 ref Complex opAssign(R)(Complex!R rhs) 37 if (!is(R == T)) 38 { 39 this.re = rhs.re; 40 this.im = rhs.im; 41 return this; 42 } 43 44 /// 45 ref Complex opAssign(F)(const F rhs) 46 if (isNumeric!F) 47 { 48 this.re = rhs; 49 this.im = 0; 50 return this; 51 } 52 53 /// 54 ref Complex opOpAssign(string op : "+", R)(Complex!R rhs) return 55 { 56 re += rhs.re; 57 im += rhs.im; 58 return this; 59 } 60 61 /// 62 ref Complex opOpAssign(string op : "-", R)(Complex!R rhs) return 63 { 64 re -= rhs.re; 65 im -= rhs.im; 66 return this; 67 } 68 69 /// 70 ref Complex opOpAssign(string op, R)(Complex!R rhs) return 71 if (op == "*" || op == "/") 72 { 73 return this = this.opBinary!op(rhs); 74 } 75 76 /// 77 ref Complex opOpAssign(string op : "+", R)(const R rhs) return 78 if (isNumeric!R) 79 { 80 re += rhs; 81 return this; 82 } 83 84 /// 85 ref Complex opOpAssign(string op : "-", R)(const R rhs) return 86 if (isNumeric!R) 87 { 88 re -= rhs; 89 return this; 90 } 91 92 /// 93 ref Complex opOpAssign(string op : "*", R)(const R rhs) return 94 if (isNumeric!R) 95 { 96 re *= rhs; 97 return this; 98 } 99 100 /// 101 ref Complex opOpAssign(string op : "/", R)(const R rhs) return 102 if (isNumeric!R) 103 { 104 re /= rhs; 105 return this; 106 } 107 108 scope const: 109 110 /// 111 bool opEquals(const Complex rhs) 112 { 113 return re == rhs.re && im == rhs.im; 114 } 115 116 /// 117 size_t toHash() 118 { 119 T[2] val = [re, im]; 120 return hashOf(val) ; 121 } 122 123 /// 124 bool opEquals(R)(Complex!R rhs) 125 if (!is(R == T)) 126 { 127 return re == rhs.re && im == rhs.im; 128 } 129 130 /// 131 bool opEquals(F)(const F rhs) 132 if (isNumeric!F) 133 { 134 return re == rhs && im == 0; 135 } 136 137 /// 138 Complex opUnary(string op : "+")() 139 { 140 return this; 141 } 142 143 /// 144 Complex opUnary(string op : "-")() 145 { 146 return typeof(return)(-re, -im); 147 } 148 149 /// 150 Complex!(CommonType!(T, R)) opBinary(string op : "+", R)(Complex!R rhs) 151 { 152 return typeof(return)(re + rhs.re, im + rhs.im); 153 } 154 155 /// 156 Complex!(CommonType!(T, R)) opBinary(string op : "-", R)(Complex!R rhs) 157 { 158 return typeof(return)(re - rhs.re, im - rhs.im); 159 } 160 161 /// 162 Complex!(CommonType!(T, R)) opBinary(string op : "*", R)(Complex!R rhs) 163 { 164 return typeof(return)(re * rhs.re - im * rhs.im, re * rhs.im + im * rhs.re); 165 } 166 167 /// 168 Complex!(CommonType!(T, R)) opBinary(string op : "/", R)(Complex!R rhs) 169 { 170 // TODO: use more precise algorithm 171 auto norm = rhs.re * rhs.re + rhs.im * rhs.im; 172 return typeof(return)( 173 (re * rhs.re + im * rhs.im) / norm, 174 (im * rhs.re - re * rhs.im) / norm, 175 ); 176 } 177 178 /// 179 Complex!(CommonType!(T, R)) opBinary(string op : "+", R)(const R rhs) 180 if (isNumeric!R) 181 { 182 return typeof(return)(re + rhs, im); 183 } 184 185 /// 186 Complex!(CommonType!(T, R)) opBinary(string op : "-", R)(const R rhs) 187 if (isNumeric!R) 188 { 189 return typeof(return)(re - rhs, im); 190 } 191 192 /// 193 Complex!(CommonType!(T, R)) opBinary(string op : "*", R)(const R rhs) 194 if (isNumeric!R) 195 { 196 return typeof(return)(re * rhs, im * rhs); 197 } 198 199 /// 200 Complex!(CommonType!(T, R)) opBinary(string op : "/", R)(const R rhs) 201 if (isNumeric!R) 202 { 203 return typeof(return)(re / rhs, im / rhs); 204 } 205 206 207 /// 208 Complex!(CommonType!(T, R)) opBinaryRight(string op : "+", R)(const R rhs) 209 if (isNumeric!R) 210 { 211 return typeof(return)(rhs + re, im); 212 } 213 214 /// 215 Complex!(CommonType!(T, R)) opBinaryRight(string op : "-", R)(const R rhs) 216 if (isNumeric!R) 217 { 218 return typeof(return)(rhs - re, -im); 219 } 220 221 /// 222 Complex!(CommonType!(T, R)) opBinaryRight(string op : "*", R)(const R rhs) 223 if (isNumeric!R) 224 { 225 return typeof(return)(rhs * re, rhs * im); 226 } 227 228 /// 229 Complex!(CommonType!(T, R)) opBinaryRight(string op : "/", R)(const R rhs) 230 if (isNumeric!R) 231 { 232 // TODO: use more precise algorithm 233 auto norm = this.re * this.re + this.im * this.im; 234 return typeof(return)( 235 rhs * (this.re / norm), 236 -rhs * (this.im / norm), 237 ); 238 } 239 240 /// 241 R opCast(R)() 242 if (isNumeric!R || isComplex!R) 243 { 244 static if (isNumeric!R) 245 return cast(R) re; 246 else 247 return R(re, im); 248 } 249 } 250 251 /// ditto 252 Complex!T complex(T)(const T re, const T im = 0) 253 if (is(T == float) || is(T == double) || is(T == real)) 254 { 255 return typeof(return)(re, im); 256 } 257 258 private alias _cdouble_ = Complex!double; 259 private alias _cfloat_ = Complex!float; 260 private alias _creal_ = Complex!real; 261 262 /// 263 unittest 264 { 265 auto a = complex(1.0, 3); 266 auto b = a; 267 b.re += 3; 268 a = b; 269 assert(a == b); 270 271 a = Complex!float(5, 6); 272 assert(a == Complex!real(5, 6)); 273 274 a += b; 275 a -= b; 276 a *= b; 277 a /= b; 278 279 a = a + b; 280 a = a - b; 281 a = a * b; 282 a = a / b; 283 284 a += 2; 285 a -= 2; 286 a *= 2; 287 a /= 2; 288 289 a = a + 2; 290 a = a - 2; 291 a = a * 2; 292 a = a / 2; 293 294 a = 2 + a; 295 a = 2 - a; 296 a = 2 * a; 297 a = 2 / a; 298 299 a = -a; 300 a = +a; 301 302 assert(a != 4.0); 303 a = 4; 304 assert(a == 4); 305 assert(cast(int)a == 4); 306 assert(cast(Complex!float)a == 4); 307 308 import std.complex : StdComplex = Complex; 309 assert(cast(StdComplex!double)a == StdComplex!double(4, 0)); 310 } 311 312 /** 313 Constructs a complex number given its absolute value and argument. 314 Params: 315 modulus = The modulus 316 argument = The argument 317 Returns: The complex number with the given modulus and argument. 318 */ 319 Complex!T fromPolar(T)(const T modulus, const T argument) 320 @safe pure nothrow @nogc 321 if (__traits(isFloating, T)) 322 { 323 import mir.math.common: sin, cos; 324 return typeof(return)(modulus * cos(argument), modulus * sin(argument)); 325 } 326 327 /// 328 @safe pure nothrow version(mir_core_test) unittest 329 { 330 import mir.math : approxEqual, PI, sqrt; 331 auto z = fromPolar(sqrt(2.0), double(PI / 4)); 332 assert(approxEqual(z.re, 1.0)); 333 assert(approxEqual(z.im, 1.0)); 334 } 335 336 /++ 337 Params: z = A complex number. 338 Returns: The complex conjugate of `z`. 339 +/ 340 Complex!T conj(T)(Complex!T z) @safe pure nothrow @nogc 341 { 342 return Complex!T(z.re, -z.im); 343 } 344 345 /// 346 @safe pure nothrow version(mir_core_test) unittest 347 { 348 assert(conj(complex(1.0)) == complex(1.0)); 349 assert(conj(complex(1.0, 2.0)) == complex(1.0, -2.0)); 350 } 351 352 /++ 353 Params: z = A complex number. 354 Returns: The argument (or phase) of `z`. 355 +/ 356 T arg(T)(Complex!T z) @safe pure nothrow @nogc 357 { 358 import std.math.trigonometry : atan2; 359 return atan2(z.im, z.re); 360 } 361 362 /// 363 @safe pure nothrow version(mir_core_test) unittest 364 { 365 import mir.math.constant: PI_2, PI_4; 366 assert(arg(complex(1.0)) == 0.0); 367 assert(arg(complex(0.0L, 1.0L)) == PI_2); 368 assert(arg(complex(1.0L, 1.0L)) == PI_4); 369 } 370 371 372 /** 373 Params: z = A complex number. 374 Returns: The absolute value (or modulus) of `z`. 375 */ 376 T cabs(T)(Complex!T z) @safe pure nothrow @nogc 377 { 378 import std.math.algebraic : hypot; 379 return hypot(z.re, z.im); 380 } 381 382 /// 383 @safe pure nothrow version(mir_core_test) unittest 384 { 385 import mir.math.common: sqrt; 386 assert(cabs(complex(1.0)) == 1.0); 387 assert(cabs(complex(0.0, 1.0)) == 1.0); 388 assert(cabs(complex(1.0L, -2.0L)) == sqrt(5.0L)); 389 } 390 391 @safe pure nothrow @nogc version(mir_core_test) unittest 392 { 393 import mir.math.common: sqrt; 394 assert(cabs(complex(0.0L, -3.2L)) == 3.2L); 395 assert(cabs(complex(0.0L, 71.6L)) == 71.6L); 396 assert(cabs(complex(-1.0L, 1.0L)) == sqrt(2.0L)); 397 }