1 // Written in the D programming language. 2 3 /** This module contains the $(LREF Complex) type, which is used to represent 4 complex numbers, along with related mathematical operations and functions. 5 6 $(LREF Complex) will eventually 7 $(DDLINK deprecate, Deprecated Features, replace) 8 the built-in types `cfloat`, `cdouble`, `creal`, `ifloat`, 9 `idouble`, and `ireal`. 10 11 Macros: 12 TABLE_SV = <table border="1" cellpadding="4" cellspacing="0"> 13 <caption>Special Values</caption> 14 $0</table> 15 PLUSMN = ± 16 NAN = $(RED NAN) 17 INFIN = ∞ 18 PI = π 19 20 Authors: Lars Tandle Kyllingstad, Don Clugston 21 Copyright: Copyright (c) 2010, Lars T. Kyllingstad. 22 License: $(HTTP boost.org/LICENSE_1_0.txt, Boost License 1.0) 23 Source: $(PHOBOSSRC std/complex.d) 24 */ 25 module std.complex; 26 27 import std.traits; 28 29 /** Helper function that returns a complex number with the specified 30 real and imaginary parts. 31 32 Params: 33 R = (template parameter) type of real part of complex number 34 I = (template parameter) type of imaginary part of complex number 35 36 re = real part of complex number to be constructed 37 im = (optional) imaginary part of complex number, 0 if omitted. 38 39 Returns: 40 `Complex` instance with real and imaginary parts set 41 to the values provided as input. If neither `re` nor 42 `im` are floating-point numbers, the return type will 43 be `Complex!double`. Otherwise, the return type is 44 deduced using $(D std.traits.CommonType!(R, I)). 45 */ 46 auto complex(R)(const R re) @safe pure nothrow @nogc 47 if (is(R : double)) 48 { 49 static if (isFloatingPoint!R) 50 return Complex!R(re, 0); 51 else 52 return Complex!double(re, 0); 53 } 54 55 /// ditto 56 auto complex(R, I)(const R re, const I im) @safe pure nothrow @nogc 57 if (is(R : double) && is(I : double)) 58 { 59 static if (isFloatingPoint!R || isFloatingPoint!I) 60 return Complex!(CommonType!(R, I))(re, im); 61 else 62 return Complex!double(re, im); 63 } 64 65 /// 66 @safe pure nothrow unittest 67 { 68 auto a = complex(1.0); 69 static assert(is(typeof(a) == Complex!double)); 70 assert(a.re == 1.0); 71 assert(a.im == 0.0); 72 73 auto b = complex(2.0L); 74 static assert(is(typeof(b) == Complex!real)); 75 assert(b.re == 2.0L); 76 assert(b.im == 0.0L); 77 78 auto c = complex(1.0, 2.0); 79 static assert(is(typeof(c) == Complex!double)); 80 assert(c.re == 1.0); 81 assert(c.im == 2.0); 82 83 auto d = complex(3.0, 4.0L); 84 static assert(is(typeof(d) == Complex!real)); 85 assert(d.re == 3.0); 86 assert(d.im == 4.0L); 87 88 auto e = complex(1); 89 static assert(is(typeof(e) == Complex!double)); 90 assert(e.re == 1); 91 assert(e.im == 0); 92 93 auto f = complex(1L, 2); 94 static assert(is(typeof(f) == Complex!double)); 95 assert(f.re == 1L); 96 assert(f.im == 2); 97 98 auto g = complex(3, 4.0L); 99 static assert(is(typeof(g) == Complex!real)); 100 assert(g.re == 3); 101 assert(g.im == 4.0L); 102 } 103 104 105 /** A complex number parametrised by a type `T`, which must be either 106 `float`, `double` or `real`. 107 */ 108 struct Complex(T) 109 if (isFloatingPoint!T) 110 { 111 import std.format : FormatSpec; 112 import std.range.primitives : isOutputRange; 113 114 /** The real part of the number. */ 115 T re; 116 117 /** The imaginary part of the number. */ 118 T im; 119 120 /** Converts the complex number to a string representation. 121 122 The second form of this function is usually not called directly; 123 instead, it is used via $(REF format, std,string), as shown in the examples 124 below. Supported format characters are 'e', 'f', 'g', 'a', and 's'. 125 126 See the $(MREF std, format) and $(REF format, std,string) 127 documentation for more information. 128 */ 129 string toString() const @safe /* TODO: pure nothrow */ 130 { 131 import std.exception : assumeUnique; 132 char[] buf; 133 buf.reserve(100); 134 auto fmt = FormatSpec!char("%s"); 135 toString((const(char)[] s) { buf ~= s; }, fmt); 136 static trustedAssumeUnique(T)(T t) @trusted { return assumeUnique(t); } 137 return trustedAssumeUnique(buf); 138 } 139 140 static if (is(T == double)) 141 /// 142 @safe unittest 143 { 144 auto c = complex(1.2, 3.4); 145 146 // Vanilla toString formatting: 147 assert(c.toString() == "1.2+3.4i"); 148 149 // Formatting with std.string.format specs: the precision and width 150 // specifiers apply to both the real and imaginary parts of the 151 // complex number. 152 import std.format : format; 153 assert(format("%.2f", c) == "1.20+3.40i"); 154 assert(format("%4.1f", c) == " 1.2+ 3.4i"); 155 } 156 157 /// ditto 158 void toString(Writer, Char)(scope Writer w, scope const ref FormatSpec!Char formatSpec) const 159 if (isOutputRange!(Writer, const(Char)[])) 160 { 161 import std.format : formatValue; 162 import std.math : signbit; 163 import std.range.primitives : put; 164 formatValue(w, re, formatSpec); 165 if (signbit(im) == 0) 166 put(w, "+"); 167 formatValue(w, im, formatSpec); 168 put(w, "i"); 169 } 170 171 @safe pure nothrow @nogc: 172 173 /** Construct a complex number with the specified real and 174 imaginary parts. In the case where a single argument is passed 175 that is not complex, the imaginary part of the result will be 176 zero. 177 */ 178 this(R : T)(Complex!R z) 179 { 180 re = z.re; 181 im = z.im; 182 } 183 184 /// ditto 185 this(Rx : T, Ry : T)(const Rx x, const Ry y) 186 { 187 re = x; 188 im = y; 189 } 190 191 /// ditto 192 this(R : T)(const R r) 193 { 194 re = r; 195 im = 0; 196 } 197 198 // ASSIGNMENT OPERATORS 199 200 // this = complex 201 ref Complex opAssign(R : T)(Complex!R z) 202 { 203 re = z.re; 204 im = z.im; 205 return this; 206 } 207 208 // this = numeric 209 ref Complex opAssign(R : T)(const R r) 210 { 211 re = r; 212 im = 0; 213 return this; 214 } 215 216 // COMPARISON OPERATORS 217 218 // this == complex 219 bool opEquals(R : T)(Complex!R z) const 220 { 221 return re == z.re && im == z.im; 222 } 223 224 // this == numeric 225 bool opEquals(R : T)(const R r) const 226 { 227 return re == r && im == 0; 228 } 229 230 // UNARY OPERATORS 231 232 // +complex 233 Complex opUnary(string op)() const 234 if (op == "+") 235 { 236 return this; 237 } 238 239 // -complex 240 Complex opUnary(string op)() const 241 if (op == "-") 242 { 243 return Complex(-re, -im); 244 } 245 246 // BINARY OPERATORS 247 248 // complex op complex 249 Complex!(CommonType!(T,R)) opBinary(string op, R)(Complex!R z) const 250 { 251 alias C = typeof(return); 252 auto w = C(this.re, this.im); 253 return w.opOpAssign!(op)(z); 254 } 255 256 // complex op numeric 257 Complex!(CommonType!(T,R)) opBinary(string op, R)(const R r) const 258 if (isNumeric!R) 259 { 260 alias C = typeof(return); 261 auto w = C(this.re, this.im); 262 return w.opOpAssign!(op)(r); 263 } 264 265 // numeric + complex, numeric * complex 266 Complex!(CommonType!(T, R)) opBinaryRight(string op, R)(const R r) const 267 if ((op == "+" || op == "*") && (isNumeric!R)) 268 { 269 return opBinary!(op)(r); 270 } 271 272 // numeric - complex 273 Complex!(CommonType!(T, R)) opBinaryRight(string op, R)(const R r) const 274 if (op == "-" && isNumeric!R) 275 { 276 return Complex(r - re, -im); 277 } 278 279 // numeric / complex 280 Complex!(CommonType!(T, R)) opBinaryRight(string op, R)(const R r) const 281 if (op == "/" && isNumeric!R) 282 { 283 import std.math : fabs; 284 typeof(return) w = void; 285 if (fabs(re) < fabs(im)) 286 { 287 immutable ratio = re/im; 288 immutable rdivd = r/(re*ratio + im); 289 290 w.re = rdivd*ratio; 291 w.im = -rdivd; 292 } 293 else 294 { 295 immutable ratio = im/re; 296 immutable rdivd = r/(re + im*ratio); 297 298 w.re = rdivd; 299 w.im = -rdivd*ratio; 300 } 301 302 return w; 303 } 304 305 // numeric ^^ complex 306 Complex!(CommonType!(T, R)) opBinaryRight(string op, R)(const R lhs) const 307 if (op == "^^" && isNumeric!R) 308 { 309 import std.math : cos, exp, log, sin, PI; 310 Unqual!(CommonType!(T, R)) ab = void, ar = void; 311 312 if (lhs >= 0) 313 { 314 // r = lhs 315 // theta = 0 316 ab = lhs ^^ this.re; 317 ar = log(lhs) * this.im; 318 } 319 else 320 { 321 // r = -lhs 322 // theta = PI 323 ab = (-lhs) ^^ this.re * exp(-PI * this.im); 324 ar = PI * this.re + log(-lhs) * this.im; 325 } 326 327 return typeof(return)(ab * cos(ar), ab * sin(ar)); 328 } 329 330 // OP-ASSIGN OPERATORS 331 332 // complex += complex, complex -= complex 333 ref Complex opOpAssign(string op, C)(const C z) 334 if ((op == "+" || op == "-") && is(C R == Complex!R)) 335 { 336 mixin ("re "~op~"= z.re;"); 337 mixin ("im "~op~"= z.im;"); 338 return this; 339 } 340 341 // complex *= complex 342 ref Complex opOpAssign(string op, C)(const C z) 343 if (op == "*" && is(C R == Complex!R)) 344 { 345 auto temp = re*z.re - im*z.im; 346 im = im*z.re + re*z.im; 347 re = temp; 348 return this; 349 } 350 351 // complex /= complex 352 ref Complex opOpAssign(string op, C)(const C z) 353 if (op == "/" && is(C R == Complex!R)) 354 { 355 import std.math : fabs; 356 if (fabs(z.re) < fabs(z.im)) 357 { 358 immutable ratio = z.re/z.im; 359 immutable denom = z.re*ratio + z.im; 360 361 immutable temp = (re*ratio + im)/denom; 362 im = (im*ratio - re)/denom; 363 re = temp; 364 } 365 else 366 { 367 immutable ratio = z.im/z.re; 368 immutable denom = z.re + z.im*ratio; 369 370 immutable temp = (re + im*ratio)/denom; 371 im = (im - re*ratio)/denom; 372 re = temp; 373 } 374 return this; 375 } 376 377 // complex ^^= complex 378 ref Complex opOpAssign(string op, C)(const C z) 379 if (op == "^^" && is(C R == Complex!R)) 380 { 381 import std.math : exp, log, cos, sin; 382 immutable r = abs(this); 383 immutable t = arg(this); 384 immutable ab = r^^z.re * exp(-t*z.im); 385 immutable ar = t*z.re + log(r)*z.im; 386 387 re = ab*cos(ar); 388 im = ab*sin(ar); 389 return this; 390 } 391 392 // complex += numeric, complex -= numeric 393 ref Complex opOpAssign(string op, U : T)(const U a) 394 if (op == "+" || op == "-") 395 { 396 mixin ("re "~op~"= a;"); 397 return this; 398 } 399 400 // complex *= numeric, complex /= numeric 401 ref Complex opOpAssign(string op, U : T)(const U a) 402 if (op == "*" || op == "/") 403 { 404 mixin ("re "~op~"= a;"); 405 mixin ("im "~op~"= a;"); 406 return this; 407 } 408 409 // complex ^^= real 410 ref Complex opOpAssign(string op, R)(const R r) 411 if (op == "^^" && isFloatingPoint!R) 412 { 413 import std.math : cos, sin; 414 immutable ab = abs(this)^^r; 415 immutable ar = arg(this)*r; 416 re = ab*cos(ar); 417 im = ab*sin(ar); 418 return this; 419 } 420 421 // complex ^^= int 422 ref Complex opOpAssign(string op, U)(const U i) 423 if (op == "^^" && isIntegral!U) 424 { 425 switch (i) 426 { 427 case 0: 428 re = 1.0; 429 im = 0.0; 430 break; 431 case 1: 432 // identity; do nothing 433 break; 434 case 2: 435 this *= this; 436 break; 437 case 3: 438 auto z = this; 439 this *= z; 440 this *= z; 441 break; 442 default: 443 this ^^= cast(real) i; 444 } 445 return this; 446 } 447 } 448 449 @safe pure nothrow unittest 450 { 451 import std.complex; 452 import std.math; 453 454 enum EPS = double.epsilon; 455 auto c1 = complex(1.0, 1.0); 456 457 // Check unary operations. 458 auto c2 = Complex!double(0.5, 2.0); 459 460 assert(c2 == +c2); 461 462 assert((-c2).re == -(c2.re)); 463 assert((-c2).im == -(c2.im)); 464 assert(c2 == -(-c2)); 465 466 // Check complex-complex operations. 467 auto cpc = c1 + c2; 468 assert(cpc.re == c1.re + c2.re); 469 assert(cpc.im == c1.im + c2.im); 470 471 auto cmc = c1 - c2; 472 assert(cmc.re == c1.re - c2.re); 473 assert(cmc.im == c1.im - c2.im); 474 475 auto ctc = c1 * c2; 476 assert(approxEqual(abs(ctc), abs(c1)*abs(c2), EPS)); 477 assert(approxEqual(arg(ctc), arg(c1)+arg(c2), EPS)); 478 479 auto cdc = c1 / c2; 480 assert(approxEqual(abs(cdc), abs(c1)/abs(c2), EPS)); 481 assert(approxEqual(arg(cdc), arg(c1)-arg(c2), EPS)); 482 483 auto cec = c1^^c2; 484 assert(approxEqual(cec.re, 0.11524131979943839881, EPS)); 485 assert(approxEqual(cec.im, 0.21870790452746026696, EPS)); 486 487 // Check complex-real operations. 488 double a = 123.456; 489 490 auto cpr = c1 + a; 491 assert(cpr.re == c1.re + a); 492 assert(cpr.im == c1.im); 493 494 auto cmr = c1 - a; 495 assert(cmr.re == c1.re - a); 496 assert(cmr.im == c1.im); 497 498 auto ctr = c1 * a; 499 assert(ctr.re == c1.re*a); 500 assert(ctr.im == c1.im*a); 501 502 auto cdr = c1 / a; 503 assert(approxEqual(abs(cdr), abs(c1)/a, EPS)); 504 assert(approxEqual(arg(cdr), arg(c1), EPS)); 505 506 auto cer = c1^^3.0; 507 assert(approxEqual(abs(cer), abs(c1)^^3, EPS)); 508 assert(approxEqual(arg(cer), arg(c1)*3, EPS)); 509 510 auto rpc = a + c1; 511 assert(rpc == cpr); 512 513 auto rmc = a - c1; 514 assert(rmc.re == a-c1.re); 515 assert(rmc.im == -c1.im); 516 517 auto rtc = a * c1; 518 assert(rtc == ctr); 519 520 auto rdc = a / c1; 521 assert(approxEqual(abs(rdc), a/abs(c1), EPS)); 522 assert(approxEqual(arg(rdc), -arg(c1), EPS)); 523 524 rdc = a / c2; 525 assert(approxEqual(abs(rdc), a/abs(c2), EPS)); 526 assert(approxEqual(arg(rdc), -arg(c2), EPS)); 527 528 auto rec1a = 1.0 ^^ c1; 529 assert(rec1a.re == 1.0); 530 assert(rec1a.im == 0.0); 531 532 auto rec2a = 1.0 ^^ c2; 533 assert(rec2a.re == 1.0); 534 assert(rec2a.im == 0.0); 535 536 auto rec1b = (-1.0) ^^ c1; 537 assert(approxEqual(abs(rec1b), std.math.exp(-PI * c1.im), EPS)); 538 auto arg1b = arg(rec1b); 539 /* The argument _should_ be PI, but floating-point rounding error 540 * means that in fact the imaginary part is very slightly negative. 541 */ 542 assert(approxEqual(arg1b, PI, EPS) || approxEqual(arg1b, -PI, EPS)); 543 544 auto rec2b = (-1.0) ^^ c2; 545 assert(approxEqual(abs(rec2b), std.math.exp(-2 * PI), EPS)); 546 assert(approxEqual(arg(rec2b), PI_2, EPS)); 547 548 auto rec3a = 0.79 ^^ complex(6.8, 5.7); 549 auto rec3b = complex(0.79, 0.0) ^^ complex(6.8, 5.7); 550 assert(approxEqual(rec3a.re, rec3b.re, EPS)); 551 assert(approxEqual(rec3a.im, rec3b.im, EPS)); 552 553 auto rec4a = (-0.79) ^^ complex(6.8, 5.7); 554 auto rec4b = complex(-0.79, 0.0) ^^ complex(6.8, 5.7); 555 assert(approxEqual(rec4a.re, rec4b.re, EPS)); 556 assert(approxEqual(rec4a.im, rec4b.im, EPS)); 557 558 auto rer = a ^^ complex(2.0, 0.0); 559 auto rcheck = a ^^ 2.0; 560 static assert(is(typeof(rcheck) == double)); 561 assert(feqrel(rer.re, rcheck) == double.mant_dig); 562 assert(isIdentical(rer.re, rcheck)); 563 assert(rer.im == 0.0); 564 565 auto rer2 = (-a) ^^ complex(2.0, 0.0); 566 rcheck = (-a) ^^ 2.0; 567 assert(feqrel(rer2.re, rcheck) == double.mant_dig); 568 assert(isIdentical(rer2.re, rcheck)); 569 assert(approxEqual(rer2.im, 0.0, EPS)); 570 571 auto rer3 = (-a) ^^ complex(-2.0, 0.0); 572 rcheck = (-a) ^^ (-2.0); 573 assert(feqrel(rer3.re, rcheck) == double.mant_dig); 574 assert(isIdentical(rer3.re, rcheck)); 575 assert(approxEqual(rer3.im, 0.0, EPS)); 576 577 auto rer4 = a ^^ complex(-2.0, 0.0); 578 rcheck = a ^^ (-2.0); 579 assert(feqrel(rer4.re, rcheck) == double.mant_dig); 580 assert(isIdentical(rer4.re, rcheck)); 581 assert(rer4.im == 0.0); 582 583 // Check Complex-int operations. 584 foreach (i; 0 .. 6) 585 { 586 auto cei = c1^^i; 587 assert(approxEqual(abs(cei), abs(c1)^^i, EPS)); 588 // Use cos() here to deal with arguments that go outside 589 // the (-pi,pi] interval (only an issue for i>3). 590 assert(approxEqual(std.math.cos(arg(cei)), std.math.cos(arg(c1)*i), EPS)); 591 } 592 593 // Check operations between different complex types. 594 auto cf = Complex!float(1.0, 1.0); 595 auto cr = Complex!real(1.0, 1.0); 596 auto c1pcf = c1 + cf; 597 auto c1pcr = c1 + cr; 598 static assert(is(typeof(c1pcf) == Complex!double)); 599 static assert(is(typeof(c1pcr) == Complex!real)); 600 assert(c1pcf.re == c1pcr.re); 601 assert(c1pcf.im == c1pcr.im); 602 603 auto c1c = c1; 604 auto c2c = c2; 605 606 c1c /= c1; 607 assert(approxEqual(c1c.re, 1.0, EPS)); 608 assert(approxEqual(c1c.im, 0.0, EPS)); 609 610 c1c = c1; 611 c1c /= c2; 612 assert(approxEqual(c1c.re, 0.588235, EPS)); 613 assert(approxEqual(c1c.im, -0.352941, EPS)); 614 615 c2c /= c1; 616 assert(approxEqual(c2c.re, 1.25, EPS)); 617 assert(approxEqual(c2c.im, 0.75, EPS)); 618 619 c2c = c2; 620 c2c /= c2; 621 assert(approxEqual(c2c.re, 1.0, EPS)); 622 assert(approxEqual(c2c.im, 0.0, EPS)); 623 } 624 625 @safe pure nothrow unittest 626 { 627 // Initialization 628 Complex!double a = 1; 629 assert(a.re == 1 && a.im == 0); 630 Complex!double b = 1.0; 631 assert(b.re == 1.0 && b.im == 0); 632 Complex!double c = Complex!real(1.0, 2); 633 assert(c.re == 1.0 && c.im == 2); 634 } 635 636 @safe pure nothrow unittest 637 { 638 // Assignments and comparisons 639 Complex!double z; 640 641 z = 1; 642 assert(z == 1); 643 assert(z.re == 1.0 && z.im == 0.0); 644 645 z = 2.0; 646 assert(z == 2.0); 647 assert(z.re == 2.0 && z.im == 0.0); 648 649 z = 1.0L; 650 assert(z == 1.0L); 651 assert(z.re == 1.0 && z.im == 0.0); 652 653 auto w = Complex!real(1.0, 1.0); 654 z = w; 655 assert(z == w); 656 assert(z.re == 1.0 && z.im == 1.0); 657 658 auto c = Complex!float(2.0, 2.0); 659 z = c; 660 assert(z == c); 661 assert(z.re == 2.0 && z.im == 2.0); 662 } 663 664 665 /* Makes Complex!(Complex!T) fold to Complex!T. 666 667 The rationale for this is that just like the real line is a 668 subspace of the complex plane, the complex plane is a subspace 669 of itself. Example of usage: 670 --- 671 Complex!T addI(T)(T x) 672 { 673 return x + Complex!T(0.0, 1.0); 674 } 675 --- 676 The above will work if T is both real and complex. 677 */ 678 template Complex(T) 679 if (is(T R == Complex!R)) 680 { 681 alias Complex = T; 682 } 683 684 @safe pure nothrow unittest 685 { 686 static assert(is(Complex!(Complex!real) == Complex!real)); 687 688 Complex!T addI(T)(T x) 689 { 690 return x + Complex!T(0.0, 1.0); 691 } 692 693 auto z1 = addI(1.0); 694 assert(z1.re == 1.0 && z1.im == 1.0); 695 696 enum one = Complex!double(1.0, 0.0); 697 auto z2 = addI(one); 698 assert(z1 == z2); 699 } 700 701 702 /** 703 Params: z = A complex number. 704 Returns: The absolute value (or modulus) of `z`. 705 */ 706 T abs(T)(Complex!T z) @safe pure nothrow @nogc 707 { 708 import std.math : hypot; 709 return hypot(z.re, z.im); 710 } 711 712 /// 713 @safe pure nothrow unittest 714 { 715 static import std.math; 716 assert(abs(complex(1.0)) == 1.0); 717 assert(abs(complex(0.0, 1.0)) == 1.0); 718 assert(abs(complex(1.0L, -2.0L)) == std.math.sqrt(5.0L)); 719 } 720 721 @safe pure nothrow @nogc unittest 722 { 723 static import std.math; 724 assert(abs(complex(0.0L, -3.2L)) == 3.2L); 725 assert(abs(complex(0.0L, 71.6L)) == 71.6L); 726 assert(abs(complex(-1.0L, 1.0L)) == std.math.sqrt(2.0L)); 727 } 728 729 @safe pure nothrow @nogc unittest 730 { 731 import std.meta : AliasSeq; 732 static foreach (T; AliasSeq!(float, double, real)) 733 {{ 734 static import std.math; 735 Complex!T a = complex(T(-12), T(3)); 736 T b = std.math.hypot(a.re, a.im); 737 assert(std.math.approxEqual(abs(a), b)); 738 assert(std.math.approxEqual(abs(-a), b)); 739 }} 740 } 741 742 /++ 743 Params: 744 z = A complex number. 745 x = A real number. 746 Returns: The squared modulus of `z`. 747 For genericity, if called on a real number, returns its square. 748 +/ 749 T sqAbs(T)(Complex!T z) @safe pure nothrow @nogc 750 { 751 return z.re*z.re + z.im*z.im; 752 } 753 754 /// 755 @safe pure nothrow unittest 756 { 757 import std.math; 758 assert(sqAbs(complex(0.0)) == 0.0); 759 assert(sqAbs(complex(1.0)) == 1.0); 760 assert(sqAbs(complex(0.0, 1.0)) == 1.0); 761 assert(approxEqual(sqAbs(complex(1.0L, -2.0L)), 5.0L)); 762 assert(approxEqual(sqAbs(complex(-3.0L, 1.0L)), 10.0L)); 763 assert(approxEqual(sqAbs(complex(1.0f,-1.0f)), 2.0f)); 764 } 765 766 /// ditto 767 T sqAbs(T)(const T x) @safe pure nothrow @nogc 768 if (isFloatingPoint!T) 769 { 770 return x*x; 771 } 772 773 @safe pure nothrow unittest 774 { 775 import std.math; 776 assert(sqAbs(0.0) == 0.0); 777 assert(sqAbs(-1.0) == 1.0); 778 assert(approxEqual(sqAbs(-3.0L), 9.0L)); 779 assert(approxEqual(sqAbs(-5.0f), 25.0f)); 780 } 781 782 783 /** 784 Params: z = A complex number. 785 Returns: The argument (or phase) of `z`. 786 */ 787 T arg(T)(Complex!T z) @safe pure nothrow @nogc 788 { 789 import std.math : atan2; 790 return atan2(z.im, z.re); 791 } 792 793 /// 794 @safe pure nothrow unittest 795 { 796 import std.math; 797 assert(arg(complex(1.0)) == 0.0); 798 assert(arg(complex(0.0L, 1.0L)) == PI_2); 799 assert(arg(complex(1.0L, 1.0L)) == PI_4); 800 } 801 802 803 /** 804 * Extracts the norm of a complex number. 805 * Params: 806 * z = A complex number 807 * Returns: 808 * The squared magnitude of `z`. 809 */ 810 T norm(T)(Complex!T z) @safe pure nothrow @nogc 811 { 812 return z.re * z.re + z.im * z.im; 813 } 814 815 /// 816 @safe pure nothrow @nogc unittest 817 { 818 import std.math : isClose, PI; 819 assert(norm(complex(3.0, 4.0)) == 25.0); 820 assert(norm(fromPolar(5.0, 0.0)) == 25.0); 821 assert(isClose(norm(fromPolar(5.0L, PI / 6)), 25.0L)); 822 assert(isClose(norm(fromPolar(5.0L, 13 * PI / 6)), 25.0L)); 823 } 824 825 826 /** 827 Params: z = A complex number. 828 Returns: The complex conjugate of `z`. 829 */ 830 Complex!T conj(T)(Complex!T z) @safe pure nothrow @nogc 831 { 832 return Complex!T(z.re, -z.im); 833 } 834 835 /// 836 @safe pure nothrow unittest 837 { 838 assert(conj(complex(1.0)) == complex(1.0)); 839 assert(conj(complex(1.0, 2.0)) == complex(1.0, -2.0)); 840 } 841 842 @safe pure nothrow @nogc unittest 843 { 844 import std.meta : AliasSeq; 845 static foreach (T; AliasSeq!(float, double, real)) 846 {{ 847 auto c = Complex!T(7, 3L); 848 assert(conj(c) == Complex!T(7, -3L)); 849 auto z = Complex!T(0, -3.2L); 850 assert(conj(z) == -z); 851 }} 852 } 853 854 /** 855 * Returns the projection of `z` onto the Riemann sphere. 856 * Params: 857 * z = A complex number 858 * Returns: 859 * The projection of `z` onto the Riemann sphere. 860 */ 861 Complex!T proj(T)(Complex!T z) 862 { 863 static import std.math; 864 865 if (std.math.isInfinity(z.re) || std.math.isInfinity(z.im)) 866 return Complex!T(T.infinity, std.math.copysign(0.0, z.im)); 867 868 return z; 869 } 870 871 /// 872 @safe pure nothrow unittest 873 { 874 assert(proj(complex(1.0)) == complex(1.0)); 875 assert(proj(complex(double.infinity, 5.0)) == complex(double.infinity, 0.0)); 876 assert(proj(complex(5.0, -double.infinity)) == complex(double.infinity, -0.0)); 877 } 878 879 880 /** 881 Constructs a complex number given its absolute value and argument. 882 Params: 883 modulus = The modulus 884 argument = The argument 885 Returns: The complex number with the given modulus and argument. 886 */ 887 Complex!(CommonType!(T, U)) fromPolar(T, U)(const T modulus, const U argument) 888 @safe pure nothrow @nogc 889 { 890 import std.math : sin, cos; 891 return Complex!(CommonType!(T,U)) 892 (modulus*cos(argument), modulus*sin(argument)); 893 } 894 895 /// 896 @safe pure nothrow unittest 897 { 898 import std.math; 899 auto z = fromPolar(std.math.sqrt(2.0), PI_4); 900 assert(approxEqual(z.re, 1.0L, real.epsilon)); 901 assert(approxEqual(z.im, 1.0L, real.epsilon)); 902 } 903 904 905 /** 906 Trigonometric functions on complex numbers. 907 908 Params: z = A complex number. 909 Returns: The sine, cosine and tangent of `z`, respectively. 910 */ 911 Complex!T sin(T)(Complex!T z) @safe pure nothrow @nogc 912 { 913 auto cs = expi(z.re); 914 auto csh = coshisinh(z.im); 915 return typeof(return)(cs.im * csh.re, cs.re * csh.im); 916 } 917 918 /// 919 @safe pure nothrow unittest 920 { 921 static import std.math; 922 assert(sin(complex(0.0)) == 0.0); 923 assert(sin(complex(2.0L, 0)) == std.math.sin(2.0L)); 924 } 925 926 927 /// ditto 928 Complex!T cos(T)(Complex!T z) @safe pure nothrow @nogc 929 { 930 auto cs = expi(z.re); 931 auto csh = coshisinh(z.im); 932 return typeof(return)(cs.re * csh.re, - cs.im * csh.im); 933 } 934 935 /// 936 @safe pure nothrow unittest 937 { 938 static import std.math; 939 assert(cos(complex(0.0)) == 1.0); 940 assert(cos(complex(1.3L, 0.0)) == std.math.cos(1.3L)); 941 assert(cos(complex(0.0L, 5.2L)) == std.math.cosh(5.2L)); 942 } 943 944 version (StdUnittest) 945 { 946 int ceqrel(T)(const Complex!T x, const Complex!T y) @safe pure nothrow @nogc 947 { 948 import std.math : feqrel; 949 const r = feqrel(x.re, y.re); 950 const i = feqrel(x.im, y.im); 951 return r < i ? r : i; 952 } 953 } 954 955 @safe pure nothrow unittest 956 { 957 static import std.math; 958 assert(ceqrel(cos(complex(0, 5.2L)), complex(std.math.cosh(5.2L), 0.0L)) >= real.mant_dig - 1); 959 assert(cos(complex(1.3L)) == std.math.cos(1.3L)); 960 } 961 962 /// ditto 963 Complex!T tan(T)(Complex!T z) @safe pure nothrow @nogc 964 { 965 return sin(z) / cos(z); 966 } 967 968 /// 969 @safe pure nothrow @nogc unittest 970 { 971 static import std.math; 972 assert(ceqrel(tan(complex(1.0, 0.0)), complex(std.math.tan(1.0), 0.0)) >= double.mant_dig - 2); 973 assert(ceqrel(tan(complex(0.0, 1.0)), complex(0.0, std.math.tanh(1.0))) >= double.mant_dig - 2); 974 } 975 976 /** 977 Params: y = A real number. 978 Returns: The value of cos(y) + i sin(y). 979 980 Note: 981 `expi` is included here for convenience and for easy migration of code. 982 */ 983 Complex!real expi(real y) @trusted pure nothrow @nogc 984 { 985 import std.math : cos, sin; 986 return Complex!real(cos(y), sin(y)); 987 } 988 989 /// 990 @safe pure nothrow unittest 991 { 992 import std.math : cos, sin; 993 assert(expi(0.0L) == 1.0L); 994 assert(expi(1.3e5L) == complex(cos(1.3e5L), sin(1.3e5L))); 995 } 996 997 /** 998 Params: y = A real number. 999 Returns: The value of cosh(y) + i sinh(y) 1000 1001 Note: 1002 `coshisinh` is included here for convenience and for easy migration of code. 1003 */ 1004 Complex!real coshisinh(real y) @safe pure nothrow @nogc 1005 { 1006 static import std.math; 1007 if (std.math.fabs(y) <= 0.5) 1008 return Complex!real(std.math.cosh(y), std.math.sinh(y)); 1009 else 1010 { 1011 auto z = std.math.exp(y); 1012 auto zi = 0.5 / z; 1013 z = 0.5 * z; 1014 return Complex!real(z + zi, z - zi); 1015 } 1016 } 1017 1018 /// 1019 @safe pure nothrow @nogc unittest 1020 { 1021 import std.math : cosh, sinh; 1022 assert(coshisinh(3.0L) == complex(cosh(3.0L), sinh(3.0L))); 1023 } 1024 1025 /** 1026 Params: z = A complex number. 1027 Returns: The square root of `z`. 1028 */ 1029 Complex!T sqrt(T)(Complex!T z) @safe pure nothrow @nogc 1030 { 1031 static import std.math; 1032 typeof(return) c; 1033 real x,y,w,r; 1034 1035 if (z == 0) 1036 { 1037 c = typeof(return)(0, 0); 1038 } 1039 else 1040 { 1041 real z_re = z.re; 1042 real z_im = z.im; 1043 1044 x = std.math.fabs(z_re); 1045 y = std.math.fabs(z_im); 1046 if (x >= y) 1047 { 1048 r = y / x; 1049 w = std.math.sqrt(x) 1050 * std.math.sqrt(0.5 * (1 + std.math.sqrt(1 + r * r))); 1051 } 1052 else 1053 { 1054 r = x / y; 1055 w = std.math.sqrt(y) 1056 * std.math.sqrt(0.5 * (r + std.math.sqrt(1 + r * r))); 1057 } 1058 1059 if (z_re >= 0) 1060 { 1061 c = typeof(return)(w, z_im / (w + w)); 1062 } 1063 else 1064 { 1065 if (z_im < 0) 1066 w = -w; 1067 c = typeof(return)(z_im / (w + w), w); 1068 } 1069 } 1070 return c; 1071 } 1072 1073 /// 1074 @safe pure nothrow unittest 1075 { 1076 static import std.math; 1077 assert(sqrt(complex(0.0)) == 0.0); 1078 assert(sqrt(complex(1.0L, 0)) == std.math.sqrt(1.0L)); 1079 assert(sqrt(complex(-1.0L, 0)) == complex(0, 1.0L)); 1080 assert(sqrt(complex(-8.0, -6.0)) == complex(1.0, -3.0)); 1081 } 1082 1083 @safe pure nothrow unittest 1084 { 1085 import std.math : approxEqual; 1086 1087 auto c1 = complex(1.0, 1.0); 1088 auto c2 = Complex!double(0.5, 2.0); 1089 1090 auto c1s = sqrt(c1); 1091 assert(approxEqual(c1s.re, 1.09868411)); 1092 assert(approxEqual(c1s.im, 0.45508986)); 1093 1094 auto c2s = sqrt(c2); 1095 assert(approxEqual(c2s.re, 1.1317134)); 1096 assert(approxEqual(c2s.im, 0.8836155)); 1097 } 1098 1099 // support %f formatting of complex numbers 1100 // https://issues.dlang.org/show_bug.cgi?id=10881 1101 @safe unittest 1102 { 1103 import std.format : format; 1104 1105 auto x = complex(1.2, 3.4); 1106 assert(format("%.2f", x) == "1.20+3.40i"); 1107 1108 auto y = complex(1.2, -3.4); 1109 assert(format("%.2f", y) == "1.20-3.40i"); 1110 } 1111 1112 @safe unittest 1113 { 1114 // Test wide string formatting 1115 import std.format; 1116 wstring wformat(T)(string format, Complex!T c) 1117 { 1118 import std.array : appender; 1119 auto w = appender!wstring(); 1120 auto n = formattedWrite(w, format, c); 1121 return w.data; 1122 } 1123 1124 auto x = complex(1.2, 3.4); 1125 assert(wformat("%.2f", x) == "1.20+3.40i"w); 1126 } 1127 1128 @safe unittest 1129 { 1130 // Test ease of use (vanilla toString() should be supported) 1131 assert(complex(1.2, 3.4).toString() == "1.2+3.4i"); 1132 } 1133 1134 @safe pure nothrow @nogc unittest 1135 { 1136 auto c = complex(3.0L, 4.0L); 1137 c = sqrt(c); 1138 assert(c.re == 2.0L); 1139 assert(c.im == 1.0L); 1140 } 1141 1142 /** 1143 * Calculates e$(SUPERSCRIPT x). 1144 * Params: 1145 * x = A complex number 1146 * Returns: 1147 * The complex base e exponential of `x` 1148 * 1149 * $(TABLE_SV 1150 * $(TR $(TH x) $(TH exp(x))) 1151 * $(TR $(TD ($(PLUSMN)0, +0)) $(TD (1, +0))) 1152 * $(TR $(TD (any, +$(INFIN))) $(TD ($(NAN), $(NAN)))) 1153 * $(TR $(TD (any, $(NAN)) $(TD ($(NAN), $(NAN))))) 1154 * $(TR $(TD (+$(INFIN), +0)) $(TD (+$(INFIN), +0))) 1155 * $(TR $(TD (-$(INFIN), any)) $(TD ($(PLUSMN)0, cis(x.im)))) 1156 * $(TR $(TD (+$(INFIN), any)) $(TD ($(PLUSMN)$(INFIN), cis(x.im)))) 1157 * $(TR $(TD (-$(INFIN), +$(INFIN))) $(TD ($(PLUSMN)0, $(PLUSMN)0))) 1158 * $(TR $(TD (+$(INFIN), +$(INFIN))) $(TD ($(PLUSMN)$(INFIN), $(NAN)))) 1159 * $(TR $(TD (-$(INFIN), $(NAN))) $(TD ($(PLUSMN)0, $(PLUSMN)0))) 1160 * $(TR $(TD (+$(INFIN), $(NAN))) $(TD ($(PLUSMN)$(INFIN), $(NAN)))) 1161 * $(TR $(TD ($(NAN), +0)) $(TD ($(NAN), +0))) 1162 * $(TR $(TD ($(NAN), any)) $(TD ($(NAN), $(NAN)))) 1163 * $(TR $(TD ($(NAN), $(NAN))) $(TD ($(NAN), $(NAN)))) 1164 * ) 1165 */ 1166 Complex!T exp(T)(Complex!T x) @trusted pure nothrow @nogc // TODO: @safe 1167 { 1168 static import std.math; 1169 1170 // Handle special cases explicitly here, as fromPolar will otherwise 1171 // cause them to return Complex!T(NaN, NaN), or with the wrong sign. 1172 if (std.math.isInfinity(x.re)) 1173 { 1174 if (std.math.isNaN(x.im)) 1175 { 1176 if (std.math.signbit(x.re)) 1177 return Complex!T(0, std.math.copysign(0, x.im)); 1178 else 1179 return x; 1180 } 1181 if (std.math.isInfinity(x.im)) 1182 { 1183 if (std.math.signbit(x.re)) 1184 return Complex!T(0, std.math.copysign(0, x.im)); 1185 else 1186 return Complex!T(T.infinity, -T.nan); 1187 } 1188 if (x.im == 0.0) 1189 { 1190 if (std.math.signbit(x.re)) 1191 return Complex!T(0.0); 1192 else 1193 return Complex!T(T.infinity); 1194 } 1195 } 1196 if (std.math.isNaN(x.re)) 1197 { 1198 if (std.math.isNaN(x.im) || std.math.isInfinity(x.im)) 1199 return Complex!T(T.nan, T.nan); 1200 if (x.im == 0.0) 1201 return x; 1202 } 1203 if (x.re == 0.0) 1204 { 1205 if (std.math.isNaN(x.im) || std.math.isInfinity(x.im)) 1206 return Complex!T(T.nan, T.nan); 1207 if (x.im == 0.0) 1208 return Complex!T(1.0, 0.0); 1209 } 1210 1211 return fromPolar!(T, T)(std.math.exp(x.re), x.im); 1212 } 1213 1214 /// 1215 @safe pure nothrow @nogc unittest 1216 { 1217 import std.math : isClose, PI; 1218 1219 assert(exp(complex(0.0, 0.0)) == complex(1.0, 0.0)); 1220 1221 auto a = complex(2.0, 1.0); 1222 assert(exp(conj(a)) == conj(exp(a))); 1223 1224 auto b = exp(complex(0.0L, 1.0L) * PI); 1225 assert(isClose(b, -1.0L, 0.0, 1e-15)); 1226 } 1227 1228 @safe pure nothrow @nogc unittest 1229 { 1230 import std.math : isNaN, isInfinity; 1231 1232 auto a = exp(complex(0.0, double.infinity)); 1233 assert(a.re.isNaN && a.im.isNaN); 1234 auto b = exp(complex(0.0, double.infinity)); 1235 assert(b.re.isNaN && b.im.isNaN); 1236 auto c = exp(complex(0.0, double.nan)); 1237 assert(c.re.isNaN && c.im.isNaN); 1238 1239 auto d = exp(complex(+double.infinity, 0.0)); 1240 assert(d == complex(double.infinity, 0.0)); 1241 auto e = exp(complex(-double.infinity, 0.0)); 1242 assert(e == complex(0.0)); 1243 auto f = exp(complex(-double.infinity, 1.0)); 1244 assert(f == complex(0.0)); 1245 auto g = exp(complex(+double.infinity, 1.0)); 1246 assert(g == complex(double.infinity, double.infinity)); 1247 auto h = exp(complex(-double.infinity, +double.infinity)); 1248 assert(h == complex(0.0)); 1249 auto i = exp(complex(+double.infinity, +double.infinity)); 1250 assert(i.re.isInfinity && i.im.isNaN); 1251 auto j = exp(complex(-double.infinity, double.nan)); 1252 assert(j == complex(0.0)); 1253 auto k = exp(complex(+double.infinity, double.nan)); 1254 assert(k.re.isInfinity && k.im.isNaN); 1255 1256 auto l = exp(complex(double.nan, 0)); 1257 assert(l.re.isNaN && l.im == 0.0); 1258 auto m = exp(complex(double.nan, 1)); 1259 assert(m.re.isNaN && m.im.isNaN); 1260 auto n = exp(complex(double.nan, double.nan)); 1261 assert(n.re.isNaN && n.im.isNaN); 1262 } 1263 1264 @safe pure nothrow @nogc unittest 1265 { 1266 import std.math : PI, isClose; 1267 1268 auto a = exp(complex(0.0, -PI)); 1269 assert(isClose(a, -1.0, 0.0, 1e-15)); 1270 1271 auto b = exp(complex(0.0, -2.0 * PI / 3.0)); 1272 assert(isClose(b, complex(-0.5L, -0.866025403784438646763L))); 1273 1274 auto c = exp(complex(0.0, PI / 3.0)); 1275 assert(isClose(c, complex(0.5L, 0.866025403784438646763L))); 1276 1277 auto d = exp(complex(0.0, 2.0 * PI / 3.0)); 1278 assert(isClose(d, complex(-0.5L, 0.866025403784438646763L))); 1279 1280 auto e = exp(complex(0.0, PI)); 1281 assert(isClose(e, -1.0, 0.0, 1e-15)); 1282 } 1283 1284 /** 1285 * Calculate the natural logarithm of x. 1286 * The branch cut is along the negative axis. 1287 * Params: 1288 * x = A complex number 1289 * Returns: 1290 * The complex natural logarithm of `x` 1291 * 1292 * $(TABLE_SV 1293 * $(TR $(TH x) $(TH log(x))) 1294 * $(TR $(TD (-0, +0)) $(TD (-$(INFIN), $(PI)))) 1295 * $(TR $(TD (+0, +0)) $(TD (-$(INFIN), +0))) 1296 * $(TR $(TD (any, +$(INFIN))) $(TD (+$(INFIN), $(PI)/2))) 1297 * $(TR $(TD (any, $(NAN))) $(TD ($(NAN), $(NAN)))) 1298 * $(TR $(TD (-$(INFIN), any)) $(TD (+$(INFIN), $(PI)))) 1299 * $(TR $(TD (+$(INFIN), any)) $(TD (+$(INFIN), +0))) 1300 * $(TR $(TD (-$(INFIN), +$(INFIN))) $(TD (+$(INFIN), 3$(PI)/4))) 1301 * $(TR $(TD (+$(INFIN), +$(INFIN))) $(TD (+$(INFIN), $(PI)/4))) 1302 * $(TR $(TD ($(PLUSMN)$(INFIN), $(NAN))) $(TD (+$(INFIN), $(NAN)))) 1303 * $(TR $(TD ($(NAN), any)) $(TD ($(NAN), $(NAN)))) 1304 * $(TR $(TD ($(NAN), +$(INFIN))) $(TD (+$(INFIN), $(NAN)))) 1305 * $(TR $(TD ($(NAN), $(NAN))) $(TD ($(NAN), $(NAN)))) 1306 * ) 1307 */ 1308 Complex!T log(T)(Complex!T x) @safe pure nothrow @nogc 1309 { 1310 static import std.math; 1311 1312 // Handle special cases explicitly here for better accuracy. 1313 // The order here is important, so that the correct path is chosen. 1314 if (std.math.isNaN(x.re)) 1315 { 1316 if (std.math.isInfinity(x.im)) 1317 return Complex!T(T.infinity, T.nan); 1318 else 1319 return Complex!T(T.nan, T.nan); 1320 } 1321 if (std.math.isInfinity(x.re)) 1322 { 1323 if (std.math.isNaN(x.im)) 1324 return Complex!T(T.infinity, T.nan); 1325 else if (std.math.isInfinity(x.im)) 1326 { 1327 if (std.math.signbit(x.re)) 1328 return Complex!T(T.infinity, std.math.copysign(3.0 * std.math.PI_4, x.im)); 1329 else 1330 return Complex!T(T.infinity, std.math.copysign(std.math.PI_4, x.im)); 1331 } 1332 else 1333 { 1334 if (std.math.signbit(x.re)) 1335 return Complex!T(T.infinity, std.math.copysign(std.math.PI, x.im)); 1336 else 1337 return Complex!T(T.infinity, std.math.copysign(0.0, x.im)); 1338 } 1339 } 1340 if (std.math.isNaN(x.im)) 1341 return Complex!T(T.nan, T.nan); 1342 if (std.math.isInfinity(x.im)) 1343 return Complex!T(T.infinity, std.math.copysign(std.math.PI_2, x.im)); 1344 if (x.re == 0.0 && x.im == 0.0) 1345 { 1346 if (std.math.signbit(x.re)) 1347 return Complex!T(-T.infinity, std.math.copysign(std.math.PI, x.im)); 1348 else 1349 return Complex!T(-T.infinity, std.math.copysign(0.0, x.im)); 1350 } 1351 1352 return Complex!T(std.math.log(abs(x)), arg(x)); 1353 } 1354 1355 /// 1356 @safe pure nothrow @nogc unittest 1357 { 1358 import std.math : sqrt, PI, isClose; 1359 1360 auto a = complex(2.0, 1.0); 1361 assert(log(conj(a)) == conj(log(a))); 1362 1363 auto b = 2.0 * log10(complex(0.0, 1.0)); 1364 auto c = 4.0 * log10(complex(sqrt(2.0) / 2, sqrt(2.0) / 2)); 1365 assert(isClose(b, c, 0.0, 1e-15)); 1366 1367 assert(log(complex(-1.0L, 0.0L)) == complex(0.0L, PI)); 1368 assert(log(complex(-1.0L, -0.0L)) == complex(0.0L, -PI)); 1369 } 1370 1371 @safe pure nothrow @nogc unittest 1372 { 1373 import std.math : isNaN, isInfinity, PI, PI_2, PI_4; 1374 1375 auto a = log(complex(-0.0L, 0.0L)); 1376 assert(a == complex(-real.infinity, PI)); 1377 auto b = log(complex(0.0L, 0.0L)); 1378 assert(b == complex(-real.infinity, +0.0L)); 1379 auto c = log(complex(1.0L, real.infinity)); 1380 assert(c == complex(real.infinity, PI_2)); 1381 auto d = log(complex(1.0L, real.nan)); 1382 assert(d.re.isNaN && d.im.isNaN); 1383 1384 auto e = log(complex(-real.infinity, 1.0L)); 1385 assert(e == complex(real.infinity, PI)); 1386 auto f = log(complex(real.infinity, 1.0L)); 1387 assert(f == complex(real.infinity, 0.0L)); 1388 auto g = log(complex(-real.infinity, real.infinity)); 1389 assert(g == complex(real.infinity, 3.0 * PI_4)); 1390 auto h = log(complex(real.infinity, real.infinity)); 1391 assert(h == complex(real.infinity, PI_4)); 1392 auto i = log(complex(real.infinity, real.nan)); 1393 assert(i.re.isInfinity && i.im.isNaN); 1394 1395 auto j = log(complex(real.nan, 1.0L)); 1396 assert(j.re.isNaN && j.im.isNaN); 1397 auto k = log(complex(real.nan, real.infinity)); 1398 assert(k.re.isInfinity && k.im.isNaN); 1399 auto l = log(complex(real.nan, real.nan)); 1400 assert(l.re.isNaN && l.im.isNaN); 1401 } 1402 1403 @safe pure nothrow @nogc unittest 1404 { 1405 import std.math : PI, isClose; 1406 1407 auto a = log(fromPolar(1.0, PI / 6.0)); 1408 assert(isClose(a, complex(0.0L, 0.523598775598298873077L), 0.0, 1e-15)); 1409 1410 auto b = log(fromPolar(1.0, PI / 3.0)); 1411 assert(isClose(b, complex(0.0L, 1.04719755119659774615L), 0.0, 1e-15)); 1412 1413 auto c = log(fromPolar(1.0, PI / 2.0)); 1414 assert(isClose(c, complex(0.0L, 1.57079632679489661923L), 0.0, 1e-15)); 1415 1416 auto d = log(fromPolar(1.0, 2.0 * PI / 3.0)); 1417 assert(isClose(d, complex(0.0L, 2.09439510239319549230L), 0.0, 1e-15)); 1418 1419 auto e = log(fromPolar(1.0, 5.0 * PI / 6.0)); 1420 assert(isClose(e, complex(0.0L, 2.61799387799149436538L), 0.0, 1e-15)); 1421 1422 auto f = log(complex(-1.0L, 0.0L)); 1423 assert(isClose(f, complex(0.0L, PI), 0.0, 1e-15)); 1424 } 1425 1426 /** 1427 * Calculate the base-10 logarithm of x. 1428 * Params: 1429 * x = A complex number 1430 * Returns: 1431 * The complex base 10 logarithm of `x` 1432 */ 1433 Complex!T log10(T)(Complex!T x) @safe pure nothrow @nogc 1434 { 1435 static import std.math; 1436 1437 return log(x) / Complex!T(std.math.log(10.0)); 1438 } 1439 1440 /// 1441 @safe pure nothrow @nogc unittest 1442 { 1443 import std.math : LN10, PI, isClose, sqrt; 1444 1445 auto a = complex(2.0, 1.0); 1446 assert(log10(a) == log(a) / log(complex(10.0))); 1447 1448 auto b = log10(complex(0.0, 1.0)) * 2.0; 1449 auto c = log10(complex(sqrt(2.0) / 2, sqrt(2.0) / 2)) * 4.0; 1450 assert(isClose(b, c, 0.0, 1e-15)); 1451 1452 assert(ceqrel(log10(complex(-100.0L, 0.0L)), complex(2.0L, PI / LN10)) >= real.mant_dig - 1); 1453 assert(ceqrel(log10(complex(-100.0L, -0.0L)), complex(2.0L, -PI / LN10)) >= real.mant_dig - 1); 1454 } 1455 1456 @safe pure nothrow @nogc unittest 1457 { 1458 import std.math : PI, isClose; 1459 1460 auto a = log10(fromPolar(1.0, PI / 6.0)); 1461 assert(isClose(a, complex(0.0L, 0.227396058973640224580L), 0.0, 1e-15)); 1462 1463 auto b = log10(fromPolar(1.0, PI / 3.0)); 1464 assert(isClose(b, complex(0.0L, 0.454792117947280449161L), 0.0, 1e-15)); 1465 1466 auto c = log10(fromPolar(1.0, PI / 2.0)); 1467 assert(isClose(c, complex(0.0L, 0.682188176920920673742L), 0.0, 1e-15)); 1468 1469 auto d = log10(fromPolar(1.0, 2.0 * PI / 3.0)); 1470 assert(isClose(d, complex(0.0L, 0.909584235894560898323L), 0.0, 1e-15)); 1471 1472 auto e = log10(fromPolar(1.0, 5.0 * PI / 6.0)); 1473 assert(isClose(e, complex(0.0L, 1.13698029486820112290L), 0.0, 1e-15)); 1474 1475 auto f = log10(complex(-1.0L, 0.0L)); 1476 assert(isClose(f, complex(0.0L, 1.36437635384184134748L), 0.0, 1e-15)); 1477 } 1478 1479 /** 1480 * Calculates x$(SUPERSCRIPT n). 1481 * The branch cut is on the negative axis. 1482 * Params: 1483 * x = base 1484 * n = exponent 1485 * Returns: 1486 * `x` raised to the power of `n` 1487 */ 1488 Complex!T pow(T, Int)(Complex!T x, const Int n) @safe pure nothrow @nogc 1489 if (isIntegral!Int) 1490 { 1491 alias UInt = Unsigned!(Unqual!Int); 1492 1493 UInt m = (n < 0) ? -cast(UInt) n : n; 1494 Complex!T y = (m % 2) ? x : Complex!T(1); 1495 1496 while (m >>= 1) 1497 { 1498 x *= x; 1499 if (m % 2) 1500 y *= x; 1501 } 1502 1503 return (n < 0) ? Complex!T(1) / y : y; 1504 } 1505 1506 /// 1507 @safe pure nothrow @nogc unittest 1508 { 1509 import std.math : isClose; 1510 1511 auto a = complex(1.0, 2.0); 1512 assert(pow(a, 2) == a * a); 1513 assert(pow(a, 3) == a * a * a); 1514 assert(pow(a, -2) == 1.0 / (a * a)); 1515 assert(isClose(pow(a, -3), 1.0 / (a * a * a))); 1516 1517 auto b = complex(2.0); 1518 assert(isClose(pow(b, 3), exp(3 * log(b)))); 1519 } 1520 1521 /// ditto 1522 Complex!T pow(T)(Complex!T x, const T n) @trusted pure nothrow @nogc 1523 { 1524 static import std.math; 1525 1526 if (x == 0.0) 1527 return Complex!T(0.0); 1528 1529 if (x.im == 0 && x.re > 0.0) 1530 return Complex!T(std.math.pow(x.re, n)); 1531 1532 Complex!T t = log(x); 1533 return fromPolar!(T, T)(std.math.exp(n * t.re), n * t.im); 1534 } 1535 1536 /// 1537 @safe pure nothrow @nogc unittest 1538 { 1539 import std.math : isClose; 1540 assert(pow(complex(0.0), 2.0) == complex(0.0)); 1541 assert(pow(complex(5.0), 2.0) == complex(25.0)); 1542 1543 auto a = pow(complex(-1.0, 0.0), 0.5); 1544 assert(isClose(a, complex(0.0, +1.0), 0.0, 1e-16)); 1545 1546 auto b = pow(complex(-1.0, -0.0), 0.5); 1547 assert(isClose(b, complex(0.0, -1.0), 0.0, 1e-16)); 1548 } 1549 1550 /// ditto 1551 Complex!T pow(T)(Complex!T x, Complex!T y) @trusted pure nothrow @nogc 1552 { 1553 return (x == 0) ? Complex!T(0) : exp(y * log(x)); 1554 } 1555 1556 /// 1557 @safe pure nothrow @nogc unittest 1558 { 1559 import std.math : isClose, exp, PI; 1560 auto a = complex(0.0); 1561 auto b = complex(2.0); 1562 assert(pow(a, b) == complex(0.0)); 1563 1564 auto c = pow(complex(0.0, 1.0), complex(0.0, 1.0)); 1565 assert(isClose(c, exp((-PI) / 2))); 1566 } 1567 1568 /// ditto 1569 Complex!T pow(T)(const T x, Complex!T n) @trusted pure nothrow @nogc 1570 { 1571 static import std.math; 1572 1573 return (x > 0.0) 1574 ? fromPolar!(T, T)(std.math.pow(x, n.re), n.im * std.math.log(x)) 1575 : pow(Complex!T(x), n); 1576 } 1577 1578 /// 1579 @safe pure nothrow @nogc unittest 1580 { 1581 import std.math : isClose; 1582 assert(pow(2.0, complex(0.0)) == complex(1.0)); 1583 assert(pow(2.0, complex(5.0)) == complex(32.0)); 1584 1585 auto a = pow(-2.0, complex(-1.0)); 1586 assert(isClose(a, complex(-0.5), 0.0, 1e-16)); 1587 1588 auto b = pow(-0.5, complex(-1.0)); 1589 assert(isClose(b, complex(-2.0), 0.0, 1e-15)); 1590 } 1591 1592 @safe pure nothrow @nogc unittest 1593 { 1594 import std.math : PI, isClose; 1595 1596 auto a = pow(complex(3.0, 4.0), 2); 1597 assert(isClose(a, complex(-7.0, 24.0))); 1598 1599 auto b = pow(complex(3.0, 4.0), PI); 1600 assert(isClose(b, complex(-152.91512205297134, 35.547499631917738))); 1601 1602 auto c = pow(complex(3.0, 4.0), complex(-2.0, 1.0)); 1603 assert(isClose(c, complex(0.015351734187477306, -0.0038407695456661503))); 1604 1605 auto d = pow(PI, complex(2.0, -1.0)); 1606 assert(isClose(d, complex(4.0790296880118296, -8.9872469554541869))); 1607 } 1608 1609 @safe pure nothrow @nogc unittest 1610 { 1611 import std.meta : AliasSeq; 1612 import std.math : RealFormat, floatTraits; 1613 static foreach (T; AliasSeq!(float, double, real)) 1614 {{ 1615 static if (floatTraits!T.realFormat == RealFormat.ibmExtended) 1616 { 1617 /* For IBM real, epsilon is too small (since 1.0 plus any double is 1618 representable) to be able to expect results within epsilon * 100. */ 1619 } 1620 else 1621 { 1622 T eps = T.epsilon * 100; 1623 1624 T a = -1.0; 1625 T b = 0.5; 1626 Complex!T ref1 = pow(complex(a), complex(b)); 1627 Complex!T res1 = pow(a, complex(b)); 1628 Complex!T res2 = pow(complex(a), b); 1629 assert(abs(ref1 - res1) < eps); 1630 assert(abs(ref1 - res2) < eps); 1631 assert(abs(res1 - res2) < eps); 1632 1633 T c = -3.2; 1634 T d = 1.4; 1635 Complex!T ref2 = pow(complex(a), complex(b)); 1636 Complex!T res3 = pow(a, complex(b)); 1637 Complex!T res4 = pow(complex(a), b); 1638 assert(abs(ref2 - res3) < eps); 1639 assert(abs(ref2 - res4) < eps); 1640 assert(abs(res3 - res4) < eps); 1641 } 1642 }} 1643 }