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*/ 25module std.complex; 26 27import 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*/ 46auto complex(R)(const R re) @safe pure nothrow @nogc 47if (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 56auto complex(R, I)(const R re, const I im) @safe pure nothrow @nogc 57if (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*/ 108struct Complex(T) 109if (isFloatingPoint!T) 110{ 111 import std.format.spec : 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.write : formatValue; 162 import std.math.traits : 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 version (FastMath) 284 { 285 // Compute norm(this) 286 immutable norm = re * re + im * im; 287 // Compute r * conj(this) 288 immutable prod_re = r * re; 289 immutable prod_im = r * -im; 290 // Divide the product by the norm 291 typeof(return) w = void; 292 w.re = prod_re / norm; 293 w.im = prod_im / norm; 294 return w; 295 } 296 else 297 { 298 import core.math : fabs; 299 typeof(return) w = void; 300 if (fabs(re) < fabs(im)) 301 { 302 immutable ratio = re/im; 303 immutable rdivd = r/(re*ratio + im); 304 305 w.re = rdivd*ratio; 306 w.im = -rdivd; 307 } 308 else 309 { 310 immutable ratio = im/re; 311 immutable rdivd = r/(re + im*ratio); 312 313 w.re = rdivd; 314 w.im = -rdivd*ratio; 315 } 316 317 return w; 318 } 319 } 320 321 // numeric ^^ complex 322 Complex!(CommonType!(T, R)) opBinaryRight(string op, R)(const R lhs) const 323 if (op == "^^" && isNumeric!R) 324 { 325 import core.math : cos, sin; 326 import std.math.exponential : exp, log; 327 import std.math.constants : PI; 328 Unqual!(CommonType!(T, R)) ab = void, ar = void; 329 330 if (lhs >= 0) 331 { 332 // r = lhs 333 // theta = 0 334 ab = lhs ^^ this.re; 335 ar = log(lhs) * this.im; 336 } 337 else 338 { 339 // r = -lhs 340 // theta = PI 341 ab = (-lhs) ^^ this.re * exp(-PI * this.im); 342 ar = PI * this.re + log(-lhs) * this.im; 343 } 344 345 return typeof(return)(ab * cos(ar), ab * sin(ar)); 346 } 347 348 // OP-ASSIGN OPERATORS 349 350 // complex += complex, complex -= complex 351 ref Complex opOpAssign(string op, C)(const C z) 352 if ((op == "+" || op == "-") && is(C R == Complex!R)) 353 { 354 mixin ("re "~op~"= z.re;"); 355 mixin ("im "~op~"= z.im;"); 356 return this; 357 } 358 359 // complex *= complex 360 ref Complex opOpAssign(string op, C)(const C z) 361 if (op == "*" && is(C R == Complex!R)) 362 { 363 auto temp = re*z.re - im*z.im; 364 im = im*z.re + re*z.im; 365 re = temp; 366 return this; 367 } 368 369 // complex /= complex 370 ref Complex opOpAssign(string op, C)(const C z) 371 if (op == "/" && is(C R == Complex!R)) 372 { 373 version (FastMath) 374 { 375 // Compute norm(z) 376 immutable norm = z.re * z.re + z.im * z.im; 377 // Compute this * conj(z) 378 immutable prod_re = re * z.re - im * -z.im; 379 immutable prod_im = im * z.re + re * -z.im; 380 // Divide the product by the norm 381 re = prod_re / norm; 382 im = prod_im / norm; 383 return this; 384 } 385 else 386 { 387 import core.math : fabs; 388 if (fabs(z.re) < fabs(z.im)) 389 { 390 immutable ratio = z.re/z.im; 391 immutable denom = z.re*ratio + z.im; 392 393 immutable temp = (re*ratio + im)/denom; 394 im = (im*ratio - re)/denom; 395 re = temp; 396 } 397 else 398 { 399 immutable ratio = z.im/z.re; 400 immutable denom = z.re + z.im*ratio; 401 402 immutable temp = (re + im*ratio)/denom; 403 im = (im - re*ratio)/denom; 404 re = temp; 405 } 406 return this; 407 } 408 } 409 410 // complex ^^= complex 411 ref Complex opOpAssign(string op, C)(const C z) 412 if (op == "^^" && is(C R == Complex!R)) 413 { 414 import core.math : cos, sin; 415 import std.math.exponential : exp, log; 416 immutable r = abs(this); 417 immutable t = arg(this); 418 immutable ab = r^^z.re * exp(-t*z.im); 419 immutable ar = t*z.re + log(r)*z.im; 420 421 re = ab*cos(ar); 422 im = ab*sin(ar); 423 return this; 424 } 425 426 // complex += numeric, complex -= numeric 427 ref Complex opOpAssign(string op, U : T)(const U a) 428 if (op == "+" || op == "-") 429 { 430 mixin ("re "~op~"= a;"); 431 return this; 432 } 433 434 // complex *= numeric, complex /= numeric 435 ref Complex opOpAssign(string op, U : T)(const U a) 436 if (op == "*" || op == "/") 437 { 438 mixin ("re "~op~"= a;"); 439 mixin ("im "~op~"= a;"); 440 return this; 441 } 442 443 // complex ^^= real 444 ref Complex opOpAssign(string op, R)(const R r) 445 if (op == "^^" && isFloatingPoint!R) 446 { 447 import core.math : cos, sin; 448 immutable ab = abs(this)^^r; 449 immutable ar = arg(this)*r; 450 re = ab*cos(ar); 451 im = ab*sin(ar); 452 return this; 453 } 454 455 // complex ^^= int 456 ref Complex opOpAssign(string op, U)(const U i) 457 if (op == "^^" && isIntegral!U) 458 { 459 switch (i) 460 { 461 case 0: 462 re = 1.0; 463 im = 0.0; 464 break; 465 case 1: 466 // identity; do nothing 467 break; 468 case 2: 469 this *= this; 470 break; 471 case 3: 472 auto z = this; 473 this *= z; 474 this *= z; 475 break; 476 default: 477 this ^^= cast(real) i; 478 } 479 return this; 480 } 481} 482 483@safe pure nothrow unittest 484{ 485 import std.complex; 486 static import core.math; 487 import std.math; 488 489 enum EPS = double.epsilon; 490 auto c1 = complex(1.0, 1.0); 491 492 // Check unary operations. 493 auto c2 = Complex!double(0.5, 2.0); 494 495 assert(c2 == +c2); 496 497 assert((-c2).re == -(c2.re)); 498 assert((-c2).im == -(c2.im)); 499 assert(c2 == -(-c2)); 500 501 // Check complex-complex operations. 502 auto cpc = c1 + c2; 503 assert(cpc.re == c1.re + c2.re); 504 assert(cpc.im == c1.im + c2.im); 505 506 auto cmc = c1 - c2; 507 assert(cmc.re == c1.re - c2.re); 508 assert(cmc.im == c1.im - c2.im); 509 510 auto ctc = c1 * c2; 511 assert(isClose(abs(ctc), abs(c1)*abs(c2), EPS)); 512 assert(isClose(arg(ctc), arg(c1)+arg(c2), EPS)); 513 514 auto cdc = c1 / c2; 515 assert(isClose(abs(cdc), abs(c1)/abs(c2), EPS)); 516 assert(isClose(arg(cdc), arg(c1)-arg(c2), EPS)); 517 518 auto cec = c1^^c2; 519 assert(isClose(cec.re, 0.1152413197994, 1e-12)); 520 assert(isClose(cec.im, 0.2187079045274, 1e-12)); 521 522 // Check complex-real operations. 523 double a = 123.456; 524 525 auto cpr = c1 + a; 526 assert(cpr.re == c1.re + a); 527 assert(cpr.im == c1.im); 528 529 auto cmr = c1 - a; 530 assert(cmr.re == c1.re - a); 531 assert(cmr.im == c1.im); 532 533 auto ctr = c1 * a; 534 assert(ctr.re == c1.re*a); 535 assert(ctr.im == c1.im*a); 536 537 auto cdr = c1 / a; 538 assert(isClose(abs(cdr), abs(c1)/a, EPS)); 539 assert(isClose(arg(cdr), arg(c1), EPS)); 540 541 auto cer = c1^^3.0; 542 assert(isClose(abs(cer), abs(c1)^^3, EPS)); 543 assert(isClose(arg(cer), arg(c1)*3, EPS)); 544 545 auto rpc = a + c1; 546 assert(rpc == cpr); 547 548 auto rmc = a - c1; 549 assert(rmc.re == a-c1.re); 550 assert(rmc.im == -c1.im); 551 552 auto rtc = a * c1; 553 assert(rtc == ctr); 554 555 auto rdc = a / c1; 556 assert(isClose(abs(rdc), a/abs(c1), EPS)); 557 assert(isClose(arg(rdc), -arg(c1), EPS)); 558 559 rdc = a / c2; 560 assert(isClose(abs(rdc), a/abs(c2), EPS)); 561 assert(isClose(arg(rdc), -arg(c2), EPS)); 562 563 auto rec1a = 1.0 ^^ c1; 564 assert(rec1a.re == 1.0); 565 assert(rec1a.im == 0.0); 566 567 auto rec2a = 1.0 ^^ c2; 568 assert(rec2a.re == 1.0); 569 assert(rec2a.im == 0.0); 570 571 auto rec1b = (-1.0) ^^ c1; 572 assert(isClose(abs(rec1b), std.math.exp(-PI * c1.im), EPS)); 573 auto arg1b = arg(rec1b); 574 /* The argument _should_ be PI, but floating-point rounding error 575 * means that in fact the imaginary part is very slightly negative. 576 */ 577 assert(isClose(arg1b, PI, EPS) || isClose(arg1b, -PI, EPS)); 578 579 auto rec2b = (-1.0) ^^ c2; 580 assert(isClose(abs(rec2b), std.math.exp(-2 * PI), EPS)); 581 assert(isClose(arg(rec2b), PI_2, EPS)); 582 583 auto rec3a = 0.79 ^^ complex(6.8, 5.7); 584 auto rec3b = complex(0.79, 0.0) ^^ complex(6.8, 5.7); 585 assert(isClose(rec3a.re, rec3b.re, 1e-14)); 586 assert(isClose(rec3a.im, rec3b.im, 1e-14)); 587 588 auto rec4a = (-0.79) ^^ complex(6.8, 5.7); 589 auto rec4b = complex(-0.79, 0.0) ^^ complex(6.8, 5.7); 590 assert(isClose(rec4a.re, rec4b.re, 1e-14)); 591 assert(isClose(rec4a.im, rec4b.im, 1e-14)); 592 593 auto rer = a ^^ complex(2.0, 0.0); 594 auto rcheck = a ^^ 2.0; 595 static assert(is(typeof(rcheck) == double)); 596 assert(feqrel(rer.re, rcheck) == double.mant_dig); 597 assert(isIdentical(rer.re, rcheck)); 598 assert(rer.im == 0.0); 599 600 auto rer2 = (-a) ^^ complex(2.0, 0.0); 601 rcheck = (-a) ^^ 2.0; 602 assert(feqrel(rer2.re, rcheck) == double.mant_dig); 603 assert(isIdentical(rer2.re, rcheck)); 604 assert(isClose(rer2.im, 0.0, 0.0, 1e-10)); 605 606 auto rer3 = (-a) ^^ complex(-2.0, 0.0); 607 rcheck = (-a) ^^ (-2.0); 608 assert(feqrel(rer3.re, rcheck) == double.mant_dig); 609 assert(isIdentical(rer3.re, rcheck)); 610 assert(isClose(rer3.im, 0.0, 0.0, EPS)); 611 612 auto rer4 = a ^^ complex(-2.0, 0.0); 613 rcheck = a ^^ (-2.0); 614 assert(feqrel(rer4.re, rcheck) == double.mant_dig); 615 assert(isIdentical(rer4.re, rcheck)); 616 assert(rer4.im == 0.0); 617 618 // Check Complex-int operations. 619 foreach (i; 0 .. 6) 620 { 621 auto cei = c1^^i; 622 assert(isClose(abs(cei), abs(c1)^^i, 1e-14)); 623 // Use cos() here to deal with arguments that go outside 624 // the (-pi,pi] interval (only an issue for i>3). 625 assert(isClose(core.math.cos(arg(cei)), core.math.cos(arg(c1)*i), 1e-14)); 626 } 627 628 // Check operations between different complex types. 629 auto cf = Complex!float(1.0, 1.0); 630 auto cr = Complex!real(1.0, 1.0); 631 auto c1pcf = c1 + cf; 632 auto c1pcr = c1 + cr; 633 static assert(is(typeof(c1pcf) == Complex!double)); 634 static assert(is(typeof(c1pcr) == Complex!real)); 635 assert(c1pcf.re == c1pcr.re); 636 assert(c1pcf.im == c1pcr.im); 637 638 auto c1c = c1; 639 auto c2c = c2; 640 641 c1c /= c1; 642 assert(isClose(c1c.re, 1.0, EPS)); 643 assert(isClose(c1c.im, 0.0, 0.0, EPS)); 644 645 c1c = c1; 646 c1c /= c2; 647 assert(isClose(c1c.re, 0.5882352941177, 1e-12)); 648 assert(isClose(c1c.im, -0.3529411764706, 1e-12)); 649 650 c2c /= c1; 651 assert(isClose(c2c.re, 1.25, EPS)); 652 assert(isClose(c2c.im, 0.75, EPS)); 653 654 c2c = c2; 655 c2c /= c2; 656 assert(isClose(c2c.re, 1.0, EPS)); 657 assert(isClose(c2c.im, 0.0, 0.0, EPS)); 658} 659 660@safe pure nothrow unittest 661{ 662 // Initialization 663 Complex!double a = 1; 664 assert(a.re == 1 && a.im == 0); 665 Complex!double b = 1.0; 666 assert(b.re == 1.0 && b.im == 0); 667 Complex!double c = Complex!real(1.0, 2); 668 assert(c.re == 1.0 && c.im == 2); 669} 670 671@safe pure nothrow unittest 672{ 673 // Assignments and comparisons 674 Complex!double z; 675 676 z = 1; 677 assert(z == 1); 678 assert(z.re == 1.0 && z.im == 0.0); 679 680 z = 2.0; 681 assert(z == 2.0); 682 assert(z.re == 2.0 && z.im == 0.0); 683 684 z = 1.0L; 685 assert(z == 1.0L); 686 assert(z.re == 1.0 && z.im == 0.0); 687 688 auto w = Complex!real(1.0, 1.0); 689 z = w; 690 assert(z == w); 691 assert(z.re == 1.0 && z.im == 1.0); 692 693 auto c = Complex!float(2.0, 2.0); 694 z = c; 695 assert(z == c); 696 assert(z.re == 2.0 && z.im == 2.0); 697} 698 699 700/* Makes Complex!(Complex!T) fold to Complex!T. 701 702 The rationale for this is that just like the real line is a 703 subspace of the complex plane, the complex plane is a subspace 704 of itself. Example of usage: 705 --- 706 Complex!T addI(T)(T x) 707 { 708 return x + Complex!T(0.0, 1.0); 709 } 710 --- 711 The above will work if T is both real and complex. 712*/ 713template Complex(T) 714if (is(T R == Complex!R)) 715{ 716 alias Complex = T; 717} 718 719@safe pure nothrow unittest 720{ 721 static assert(is(Complex!(Complex!real) == Complex!real)); 722 723 Complex!T addI(T)(T x) 724 { 725 return x + Complex!T(0.0, 1.0); 726 } 727 728 auto z1 = addI(1.0); 729 assert(z1.re == 1.0 && z1.im == 1.0); 730 731 enum one = Complex!double(1.0, 0.0); 732 auto z2 = addI(one); 733 assert(z1 == z2); 734} 735 736 737/** 738 Params: z = A complex number. 739 Returns: The absolute value (or modulus) of `z`. 740*/ 741T abs(T)(Complex!T z) @safe pure nothrow @nogc 742{ 743 import std.math.algebraic : hypot; 744 return hypot(z.re, z.im); 745} 746 747/// 748@safe pure nothrow unittest 749{ 750 static import core.math; 751 assert(abs(complex(1.0)) == 1.0); 752 assert(abs(complex(0.0, 1.0)) == 1.0); 753 assert(abs(complex(1.0L, -2.0L)) == core.math.sqrt(5.0L)); 754} 755 756@safe pure nothrow @nogc unittest 757{ 758 static import core.math; 759 assert(abs(complex(0.0L, -3.2L)) == 3.2L); 760 assert(abs(complex(0.0L, 71.6L)) == 71.6L); 761 assert(abs(complex(-1.0L, 1.0L)) == core.math.sqrt(2.0L)); 762} 763 764@safe pure nothrow @nogc unittest 765{ 766 import std.meta : AliasSeq; 767 static foreach (T; AliasSeq!(float, double, real)) 768 {{ 769 static import std.math; 770 Complex!T a = complex(T(-12), T(3)); 771 T b = std.math.hypot(a.re, a.im); 772 assert(std.math.isClose(abs(a), b)); 773 assert(std.math.isClose(abs(-a), b)); 774 }} 775} 776 777/++ 778 Params: 779 z = A complex number. 780 x = A real number. 781 Returns: The squared modulus of `z`. 782 For genericity, if called on a real number, returns its square. 783+/ 784T sqAbs(T)(Complex!T z) @safe pure nothrow @nogc 785{ 786 return z.re*z.re + z.im*z.im; 787} 788 789/// 790@safe pure nothrow unittest 791{ 792 import std.math.operations : isClose; 793 assert(sqAbs(complex(0.0)) == 0.0); 794 assert(sqAbs(complex(1.0)) == 1.0); 795 assert(sqAbs(complex(0.0, 1.0)) == 1.0); 796 assert(isClose(sqAbs(complex(1.0L, -2.0L)), 5.0L)); 797 assert(isClose(sqAbs(complex(-3.0L, 1.0L)), 10.0L)); 798 assert(isClose(sqAbs(complex(1.0f,-1.0f)), 2.0f)); 799} 800 801/// ditto 802T sqAbs(T)(const T x) @safe pure nothrow @nogc 803if (isFloatingPoint!T) 804{ 805 return x*x; 806} 807 808@safe pure nothrow unittest 809{ 810 import std.math.operations : isClose; 811 assert(sqAbs(0.0) == 0.0); 812 assert(sqAbs(-1.0) == 1.0); 813 assert(isClose(sqAbs(-3.0L), 9.0L)); 814 assert(isClose(sqAbs(-5.0f), 25.0f)); 815} 816 817 818/** 819 Params: z = A complex number. 820 Returns: The argument (or phase) of `z`. 821 */ 822T arg(T)(Complex!T z) @safe pure nothrow @nogc 823{ 824 import std.math.trigonometry : atan2; 825 return atan2(z.im, z.re); 826} 827 828/// 829@safe pure nothrow unittest 830{ 831 import std.math.constants : PI_2, PI_4; 832 assert(arg(complex(1.0)) == 0.0); 833 assert(arg(complex(0.0L, 1.0L)) == PI_2); 834 assert(arg(complex(1.0L, 1.0L)) == PI_4); 835} 836 837 838/** 839 * Extracts the norm of a complex number. 840 * Params: 841 * z = A complex number 842 * Returns: 843 * The squared magnitude of `z`. 844 */ 845T norm(T)(Complex!T z) @safe pure nothrow @nogc 846{ 847 return z.re * z.re + z.im * z.im; 848} 849 850/// 851@safe pure nothrow @nogc unittest 852{ 853 import std.math.operations : isClose; 854 import std.math.constants : PI; 855 assert(norm(complex(3.0, 4.0)) == 25.0); 856 assert(norm(fromPolar(5.0, 0.0)) == 25.0); 857 assert(isClose(norm(fromPolar(5.0L, PI / 6)), 25.0L)); 858 assert(isClose(norm(fromPolar(5.0L, 13 * PI / 6)), 25.0L)); 859} 860 861 862/** 863 Params: z = A complex number. 864 Returns: The complex conjugate of `z`. 865*/ 866Complex!T conj(T)(Complex!T z) @safe pure nothrow @nogc 867{ 868 return Complex!T(z.re, -z.im); 869} 870 871/// 872@safe pure nothrow unittest 873{ 874 assert(conj(complex(1.0)) == complex(1.0)); 875 assert(conj(complex(1.0, 2.0)) == complex(1.0, -2.0)); 876} 877 878@safe pure nothrow @nogc unittest 879{ 880 import std.meta : AliasSeq; 881 static foreach (T; AliasSeq!(float, double, real)) 882 {{ 883 auto c = Complex!T(7, 3L); 884 assert(conj(c) == Complex!T(7, -3L)); 885 auto z = Complex!T(0, -3.2L); 886 assert(conj(z) == -z); 887 }} 888} 889 890/** 891 * Returns the projection of `z` onto the Riemann sphere. 892 * Params: 893 * z = A complex number 894 * Returns: 895 * The projection of `z` onto the Riemann sphere. 896 */ 897Complex!T proj(T)(Complex!T z) 898{ 899 static import std.math; 900 901 if (std.math.isInfinity(z.re) || std.math.isInfinity(z.im)) 902 return Complex!T(T.infinity, std.math.copysign(0.0, z.im)); 903 904 return z; 905} 906 907/// 908@safe pure nothrow unittest 909{ 910 assert(proj(complex(1.0)) == complex(1.0)); 911 assert(proj(complex(double.infinity, 5.0)) == complex(double.infinity, 0.0)); 912 assert(proj(complex(5.0, -double.infinity)) == complex(double.infinity, -0.0)); 913} 914 915 916/** 917 Constructs a complex number given its absolute value and argument. 918 Params: 919 modulus = The modulus 920 argument = The argument 921 Returns: The complex number with the given modulus and argument. 922*/ 923Complex!(CommonType!(T, U)) fromPolar(T, U)(const T modulus, const U argument) 924 @safe pure nothrow @nogc 925{ 926 import core.math : sin, cos; 927 return Complex!(CommonType!(T,U)) 928 (modulus*cos(argument), modulus*sin(argument)); 929} 930 931/// 932@safe pure nothrow unittest 933{ 934 import core.math; 935 import std.math.operations : isClose; 936 import std.math.algebraic : sqrt; 937 import std.math.constants : PI_4; 938 auto z = fromPolar(core.math.sqrt(2.0), PI_4); 939 assert(isClose(z.re, 1.0L)); 940 assert(isClose(z.im, 1.0L)); 941} 942 943version (StdUnittest) 944{ 945 // Helper function for comparing two Complex numbers. 946 int ceqrel(T)(const Complex!T x, const Complex!T y) @safe pure nothrow @nogc 947 { 948 import std.math.operations : 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/** 956 Trigonometric functions on complex numbers. 957 958 Params: z = A complex number. 959 Returns: The sine, cosine and tangent of `z`, respectively. 960*/ 961Complex!T sin(T)(Complex!T z) @safe pure nothrow @nogc 962{ 963 auto cs = expi(z.re); 964 auto csh = coshisinh(z.im); 965 return typeof(return)(cs.im * csh.re, cs.re * csh.im); 966} 967 968/// 969@safe pure nothrow unittest 970{ 971 static import core.math; 972 assert(sin(complex(0.0)) == 0.0); 973 assert(sin(complex(2.0, 0)) == core.math.sin(2.0)); 974} 975 976@safe pure nothrow unittest 977{ 978 static import core.math; 979 assert(ceqrel(sin(complex(2.0L, 0)), complex(core.math.sin(2.0L))) >= real.mant_dig - 1); 980} 981 982/// ditto 983Complex!T cos(T)(Complex!T z) @safe pure nothrow @nogc 984{ 985 auto cs = expi(z.re); 986 auto csh = coshisinh(z.im); 987 return typeof(return)(cs.re * csh.re, - cs.im * csh.im); 988} 989 990/// 991@safe pure nothrow unittest 992{ 993 static import core.math; 994 static import std.math; 995 assert(cos(complex(0.0)) == 1.0); 996 assert(cos(complex(1.3, 0.0)) == core.math.cos(1.3)); 997 assert(cos(complex(0.0, 5.2)) == std.math.cosh(5.2)); 998} 999 1000@safe pure nothrow unittest 1001{ 1002 static import core.math; 1003 static import std.math; 1004 assert(ceqrel(cos(complex(0, 5.2L)), complex(std.math.cosh(5.2L), 0.0L)) >= real.mant_dig - 1); 1005 assert(ceqrel(cos(complex(1.3L)), complex(core.math.cos(1.3L))) >= real.mant_dig - 1); 1006} 1007 1008/// ditto 1009Complex!T tan(T)(Complex!T z) @safe pure nothrow @nogc 1010{ 1011 return sin(z) / cos(z); 1012} 1013 1014/// 1015@safe pure nothrow @nogc unittest 1016{ 1017 static import std.math; 1018 1019 int ceqrel(T)(const Complex!T x, const Complex!T y) @safe pure nothrow @nogc 1020 { 1021 import std.math.operations : feqrel; 1022 const r = feqrel(x.re, y.re); 1023 const i = feqrel(x.im, y.im); 1024 return r < i ? r : i; 1025 } 1026 assert(ceqrel(tan(complex(1.0, 0.0)), complex(std.math.tan(1.0), 0.0)) >= double.mant_dig - 2); 1027 assert(ceqrel(tan(complex(0.0, 1.0)), complex(0.0, std.math.tanh(1.0))) >= double.mant_dig - 2); 1028} 1029 1030/** 1031 Inverse trigonometric functions on complex numbers. 1032 1033 Params: z = A complex number. 1034 Returns: The arcsine, arccosine and arctangent of `z`, respectively. 1035*/ 1036Complex!T asin(T)(Complex!T z) @safe pure nothrow @nogc 1037{ 1038 auto ash = asinh(Complex!T(-z.im, z.re)); 1039 return Complex!T(ash.im, -ash.re); 1040} 1041 1042/// 1043@safe pure nothrow unittest 1044{ 1045 import std.math.operations : isClose; 1046 import std.math.constants : PI; 1047 assert(asin(complex(0.0)) == 0.0); 1048 assert(isClose(asin(complex(0.5L)), PI / 6)); 1049} 1050 1051@safe pure nothrow unittest 1052{ 1053 import std.math.operations : isClose; 1054 import std.math.constants : PI; 1055 version (DigitalMars) {} else // Disabled because of issue 21376 1056 assert(isClose(asin(complex(0.5f)), float(PI) / 6)); 1057} 1058 1059/// ditto 1060Complex!T acos(T)(Complex!T z) @safe pure nothrow @nogc 1061{ 1062 static import std.math; 1063 auto as = asin(z); 1064 return Complex!T(T(std.math.PI_2) - as.re, as.im); 1065} 1066 1067/// 1068@safe pure nothrow unittest 1069{ 1070 import std.math.operations : isClose; 1071 import std.math.constants : PI; 1072 import std.math.trigonometry : std_math_acos = acos; 1073 assert(acos(complex(0.0)) == std_math_acos(0.0)); 1074 assert(isClose(acos(complex(0.5L)), PI / 3)); 1075} 1076 1077@safe pure nothrow unittest 1078{ 1079 import std.math.operations : isClose; 1080 import std.math.constants : PI; 1081 version (DigitalMars) {} else // Disabled because of issue 21376 1082 assert(isClose(acos(complex(0.5f)), float(PI) / 3)); 1083} 1084 1085/// ditto 1086Complex!T atan(T)(Complex!T z) @safe pure nothrow @nogc 1087{ 1088 static import std.math; 1089 const T re2 = z.re * z.re; 1090 const T x = 1 - re2 - z.im * z.im; 1091 1092 T num = z.im + 1; 1093 T den = z.im - 1; 1094 1095 num = re2 + num * num; 1096 den = re2 + den * den; 1097 1098 return Complex!T(T(0.5) * std.math.atan2(2 * z.re, x), 1099 T(0.25) * std.math.log(num / den)); 1100} 1101 1102/// 1103@safe pure nothrow @nogc unittest 1104{ 1105 import std.math.operations : isClose; 1106 import std.math.constants : PI; 1107 assert(atan(complex(0.0)) == 0.0); 1108 assert(isClose(atan(sqrt(complex(3.0L))), PI / 3)); 1109 assert(isClose(atan(sqrt(complex(3.0f))), float(PI) / 3)); 1110} 1111 1112/** 1113 Hyperbolic trigonometric functions on complex numbers. 1114 1115 Params: z = A complex number. 1116 Returns: The hyperbolic sine, cosine and tangent of `z`, respectively. 1117*/ 1118Complex!T sinh(T)(Complex!T z) @safe pure nothrow @nogc 1119{ 1120 static import core.math, std.math; 1121 return Complex!T(std.math.sinh(z.re) * core.math.cos(z.im), 1122 std.math.cosh(z.re) * core.math.sin(z.im)); 1123} 1124 1125/// 1126@safe pure nothrow unittest 1127{ 1128 static import std.math; 1129 assert(sinh(complex(0.0)) == 0.0); 1130 assert(sinh(complex(1.0L)) == std.math.sinh(1.0L)); 1131 assert(sinh(complex(1.0f)) == std.math.sinh(1.0f)); 1132} 1133 1134/// ditto 1135Complex!T cosh(T)(Complex!T z) @safe pure nothrow @nogc 1136{ 1137 static import core.math, std.math; 1138 return Complex!T(std.math.cosh(z.re) * core.math.cos(z.im), 1139 std.math.sinh(z.re) * core.math.sin(z.im)); 1140} 1141 1142/// 1143@safe pure nothrow unittest 1144{ 1145 static import std.math; 1146 assert(cosh(complex(0.0)) == 1.0); 1147 assert(cosh(complex(1.0L)) == std.math.cosh(1.0L)); 1148 assert(cosh(complex(1.0f)) == std.math.cosh(1.0f)); 1149} 1150 1151/// ditto 1152Complex!T tanh(T)(Complex!T z) @safe pure nothrow @nogc 1153{ 1154 return sinh(z) / cosh(z); 1155} 1156 1157/// 1158@safe pure nothrow @nogc unittest 1159{ 1160 import std.math.operations : isClose; 1161 import std.math.trigonometry : std_math_tanh = tanh; 1162 assert(tanh(complex(0.0)) == 0.0); 1163 assert(isClose(tanh(complex(1.0L)), std_math_tanh(1.0L))); 1164 assert(isClose(tanh(complex(1.0f)), std_math_tanh(1.0f))); 1165} 1166 1167/** 1168 Inverse hyperbolic trigonometric functions on complex numbers. 1169 1170 Params: z = A complex number. 1171 Returns: The hyperbolic arcsine, arccosine and arctangent of `z`, respectively. 1172*/ 1173Complex!T asinh(T)(Complex!T z) @safe pure nothrow @nogc 1174{ 1175 auto t = Complex!T((z.re - z.im) * (z.re + z.im) + 1, 2 * z.re * z.im); 1176 return log(sqrt(t) + z); 1177} 1178 1179/// 1180@safe pure nothrow unittest 1181{ 1182 import std.math.operations : isClose; 1183 import std.math.trigonometry : std_math_asinh = asinh; 1184 assert(asinh(complex(0.0)) == 0.0); 1185 assert(isClose(asinh(complex(1.0L)), std_math_asinh(1.0L))); 1186 assert(isClose(asinh(complex(1.0f)), std_math_asinh(1.0f))); 1187} 1188 1189/// ditto 1190Complex!T acosh(T)(Complex!T z) @safe pure nothrow @nogc 1191{ 1192 return 2 * log(sqrt(T(0.5) * (z + 1)) + sqrt(T(0.5) * (z - 1))); 1193} 1194 1195/// 1196@safe pure nothrow unittest 1197{ 1198 import std.math.operations : isClose; 1199 import std.math.trigonometry : std_math_acosh = acosh; 1200 assert(acosh(complex(1.0)) == 0.0); 1201 assert(isClose(acosh(complex(3.0L)), std_math_acosh(3.0L))); 1202 assert(isClose(acosh(complex(3.0f)), std_math_acosh(3.0f))); 1203} 1204 1205/// ditto 1206Complex!T atanh(T)(Complex!T z) @safe pure nothrow @nogc 1207{ 1208 static import std.math; 1209 const T im2 = z.im * z.im; 1210 const T x = 1 - im2 - z.re * z.re; 1211 1212 T num = 1 + z.re; 1213 T den = 1 - z.re; 1214 1215 num = im2 + num * num; 1216 den = im2 + den * den; 1217 1218 return Complex!T(T(0.25) * (std.math.log(num) - std.math.log(den)), 1219 T(0.5) * std.math.atan2(2 * z.im, x)); 1220} 1221 1222/// 1223@safe pure nothrow @nogc unittest 1224{ 1225 import std.math.operations : isClose; 1226 import std.math.trigonometry : std_math_atanh = atanh; 1227 assert(atanh(complex(0.0)) == 0.0); 1228 assert(isClose(atanh(complex(0.5L)), std_math_atanh(0.5L))); 1229 assert(isClose(atanh(complex(0.5f)), std_math_atanh(0.5f))); 1230} 1231 1232/** 1233 Params: y = A real number. 1234 Returns: The value of cos(y) + i sin(y). 1235 1236 Note: 1237 `expi` is included here for convenience and for easy migration of code. 1238*/ 1239Complex!real expi(real y) @trusted pure nothrow @nogc 1240{ 1241 import core.math : cos, sin; 1242 return Complex!real(cos(y), sin(y)); 1243} 1244 1245/// 1246@safe pure nothrow unittest 1247{ 1248 import core.math : cos, sin; 1249 assert(expi(0.0L) == 1.0L); 1250 assert(expi(1.3e5L) == complex(cos(1.3e5L), sin(1.3e5L))); 1251} 1252 1253/** 1254 Params: y = A real number. 1255 Returns: The value of cosh(y) + i sinh(y) 1256 1257 Note: 1258 `coshisinh` is included here for convenience and for easy migration of code. 1259*/ 1260Complex!real coshisinh(real y) @safe pure nothrow @nogc 1261{ 1262 static import core.math; 1263 static import std.math; 1264 if (core.math.fabs(y) <= 0.5) 1265 return Complex!real(std.math.cosh(y), std.math.sinh(y)); 1266 else 1267 { 1268 auto z = std.math.exp(y); 1269 auto zi = 0.5 / z; 1270 z = 0.5 * z; 1271 return Complex!real(z + zi, z - zi); 1272 } 1273} 1274 1275/// 1276@safe pure nothrow @nogc unittest 1277{ 1278 import std.math.trigonometry : cosh, sinh; 1279 assert(coshisinh(3.0L) == complex(cosh(3.0L), sinh(3.0L))); 1280} 1281 1282/** 1283 Params: z = A complex number. 1284 Returns: The square root of `z`. 1285*/ 1286Complex!T sqrt(T)(Complex!T z) @safe pure nothrow @nogc 1287{ 1288 static import core.math; 1289 typeof(return) c; 1290 real x,y,w,r; 1291 1292 if (z == 0) 1293 { 1294 c = typeof(return)(0, 0); 1295 } 1296 else 1297 { 1298 real z_re = z.re; 1299 real z_im = z.im; 1300 1301 x = core.math.fabs(z_re); 1302 y = core.math.fabs(z_im); 1303 if (x >= y) 1304 { 1305 r = y / x; 1306 w = core.math.sqrt(x) 1307 * core.math.sqrt(0.5 * (1 + core.math.sqrt(1 + r * r))); 1308 } 1309 else 1310 { 1311 r = x / y; 1312 w = core.math.sqrt(y) 1313 * core.math.sqrt(0.5 * (r + core.math.sqrt(1 + r * r))); 1314 } 1315 1316 if (z_re >= 0) 1317 { 1318 c = typeof(return)(w, z_im / (w + w)); 1319 } 1320 else 1321 { 1322 if (z_im < 0) 1323 w = -w; 1324 c = typeof(return)(z_im / (w + w), w); 1325 } 1326 } 1327 return c; 1328} 1329 1330/// 1331@safe pure nothrow unittest 1332{ 1333 static import core.math; 1334 assert(sqrt(complex(0.0)) == 0.0); 1335 assert(sqrt(complex(1.0L, 0)) == core.math.sqrt(1.0L)); 1336 assert(sqrt(complex(-1.0L, 0)) == complex(0, 1.0L)); 1337 assert(sqrt(complex(-8.0, -6.0)) == complex(1.0, -3.0)); 1338} 1339 1340@safe pure nothrow unittest 1341{ 1342 import std.math.operations : isClose; 1343 1344 auto c1 = complex(1.0, 1.0); 1345 auto c2 = Complex!double(0.5, 2.0); 1346 1347 auto c1s = sqrt(c1); 1348 assert(isClose(c1s.re, 1.09868411347)); 1349 assert(isClose(c1s.im, 0.455089860562)); 1350 1351 auto c2s = sqrt(c2); 1352 assert(isClose(c2s.re, 1.13171392428)); 1353 assert(isClose(c2s.im, 0.883615530876)); 1354} 1355 1356// support %f formatting of complex numbers 1357// https://issues.dlang.org/show_bug.cgi?id=10881 1358@safe unittest 1359{ 1360 import std.format : format; 1361 1362 auto x = complex(1.2, 3.4); 1363 assert(format("%.2f", x) == "1.20+3.40i"); 1364 1365 auto y = complex(1.2, -3.4); 1366 assert(format("%.2f", y) == "1.20-3.40i"); 1367} 1368 1369@safe unittest 1370{ 1371 // Test wide string formatting 1372 import std.format.write : formattedWrite; 1373 wstring wformat(T)(string format, Complex!T c) 1374 { 1375 import std.array : appender; 1376 auto w = appender!wstring(); 1377 auto n = formattedWrite(w, format, c); 1378 return w.data; 1379 } 1380 1381 auto x = complex(1.2, 3.4); 1382 assert(wformat("%.2f", x) == "1.20+3.40i"w); 1383} 1384 1385@safe unittest 1386{ 1387 // Test ease of use (vanilla toString() should be supported) 1388 assert(complex(1.2, 3.4).toString() == "1.2+3.4i"); 1389} 1390 1391@safe pure nothrow @nogc unittest 1392{ 1393 auto c = complex(3.0L, 4.0L); 1394 c = sqrt(c); 1395 assert(c.re == 2.0L); 1396 assert(c.im == 1.0L); 1397} 1398 1399/** 1400 * Calculates e$(SUPERSCRIPT x). 1401 * Params: 1402 * x = A complex number 1403 * Returns: 1404 * The complex base e exponential of `x` 1405 * 1406 * $(TABLE_SV 1407 * $(TR $(TH x) $(TH exp(x))) 1408 * $(TR $(TD ($(PLUSMN)0, +0)) $(TD (1, +0))) 1409 * $(TR $(TD (any, +$(INFIN))) $(TD ($(NAN), $(NAN)))) 1410 * $(TR $(TD (any, $(NAN)) $(TD ($(NAN), $(NAN))))) 1411 * $(TR $(TD (+$(INFIN), +0)) $(TD (+$(INFIN), +0))) 1412 * $(TR $(TD (-$(INFIN), any)) $(TD ($(PLUSMN)0, cis(x.im)))) 1413 * $(TR $(TD (+$(INFIN), any)) $(TD ($(PLUSMN)$(INFIN), cis(x.im)))) 1414 * $(TR $(TD (-$(INFIN), +$(INFIN))) $(TD ($(PLUSMN)0, $(PLUSMN)0))) 1415 * $(TR $(TD (+$(INFIN), +$(INFIN))) $(TD ($(PLUSMN)$(INFIN), $(NAN)))) 1416 * $(TR $(TD (-$(INFIN), $(NAN))) $(TD ($(PLUSMN)0, $(PLUSMN)0))) 1417 * $(TR $(TD (+$(INFIN), $(NAN))) $(TD ($(PLUSMN)$(INFIN), $(NAN)))) 1418 * $(TR $(TD ($(NAN), +0)) $(TD ($(NAN), +0))) 1419 * $(TR $(TD ($(NAN), any)) $(TD ($(NAN), $(NAN)))) 1420 * $(TR $(TD ($(NAN), $(NAN))) $(TD ($(NAN), $(NAN)))) 1421 * ) 1422 */ 1423Complex!T exp(T)(Complex!T x) @trusted pure nothrow @nogc // TODO: @safe 1424{ 1425 static import std.math; 1426 1427 // Handle special cases explicitly here, as fromPolar will otherwise 1428 // cause them to return Complex!T(NaN, NaN), or with the wrong sign. 1429 if (std.math.isInfinity(x.re)) 1430 { 1431 if (std.math.isNaN(x.im)) 1432 { 1433 if (std.math.signbit(x.re)) 1434 return Complex!T(0, std.math.copysign(0, x.im)); 1435 else 1436 return x; 1437 } 1438 if (std.math.isInfinity(x.im)) 1439 { 1440 if (std.math.signbit(x.re)) 1441 return Complex!T(0, std.math.copysign(0, x.im)); 1442 else 1443 return Complex!T(T.infinity, -T.nan); 1444 } 1445 if (x.im == 0.0) 1446 { 1447 if (std.math.signbit(x.re)) 1448 return Complex!T(0.0); 1449 else 1450 return Complex!T(T.infinity); 1451 } 1452 } 1453 if (std.math.isNaN(x.re)) 1454 { 1455 if (std.math.isNaN(x.im) || std.math.isInfinity(x.im)) 1456 return Complex!T(T.nan, T.nan); 1457 if (x.im == 0.0) 1458 return x; 1459 } 1460 if (x.re == 0.0) 1461 { 1462 if (std.math.isNaN(x.im) || std.math.isInfinity(x.im)) 1463 return Complex!T(T.nan, T.nan); 1464 if (x.im == 0.0) 1465 return Complex!T(1.0, 0.0); 1466 } 1467 1468 return fromPolar!(T, T)(std.math.exp(x.re), x.im); 1469} 1470 1471/// 1472@safe pure nothrow @nogc unittest 1473{ 1474 import std.math.operations : isClose; 1475 import std.math.constants : PI; 1476 1477 assert(exp(complex(0.0, 0.0)) == complex(1.0, 0.0)); 1478 1479 auto a = complex(2.0, 1.0); 1480 assert(exp(conj(a)) == conj(exp(a))); 1481 1482 auto b = exp(complex(0.0L, 1.0L) * PI); 1483 assert(isClose(b, -1.0L, 0.0, 1e-15)); 1484} 1485 1486@safe pure nothrow @nogc unittest 1487{ 1488 import std.math.traits : isNaN, isInfinity; 1489 1490 auto a = exp(complex(0.0, double.infinity)); 1491 assert(a.re.isNaN && a.im.isNaN); 1492 auto b = exp(complex(0.0, double.infinity)); 1493 assert(b.re.isNaN && b.im.isNaN); 1494 auto c = exp(complex(0.0, double.nan)); 1495 assert(c.re.isNaN && c.im.isNaN); 1496 1497 auto d = exp(complex(+double.infinity, 0.0)); 1498 assert(d == complex(double.infinity, 0.0)); 1499 auto e = exp(complex(-double.infinity, 0.0)); 1500 assert(e == complex(0.0)); 1501 auto f = exp(complex(-double.infinity, 1.0)); 1502 assert(f == complex(0.0)); 1503 auto g = exp(complex(+double.infinity, 1.0)); 1504 assert(g == complex(double.infinity, double.infinity)); 1505 auto h = exp(complex(-double.infinity, +double.infinity)); 1506 assert(h == complex(0.0)); 1507 auto i = exp(complex(+double.infinity, +double.infinity)); 1508 assert(i.re.isInfinity && i.im.isNaN); 1509 auto j = exp(complex(-double.infinity, double.nan)); 1510 assert(j == complex(0.0)); 1511 auto k = exp(complex(+double.infinity, double.nan)); 1512 assert(k.re.isInfinity && k.im.isNaN); 1513 1514 auto l = exp(complex(double.nan, 0)); 1515 assert(l.re.isNaN && l.im == 0.0); 1516 auto m = exp(complex(double.nan, 1)); 1517 assert(m.re.isNaN && m.im.isNaN); 1518 auto n = exp(complex(double.nan, double.nan)); 1519 assert(n.re.isNaN && n.im.isNaN); 1520} 1521 1522@safe pure nothrow @nogc unittest 1523{ 1524 import std.math.constants : PI; 1525 import std.math.operations : isClose; 1526 1527 auto a = exp(complex(0.0, -PI)); 1528 assert(isClose(a, -1.0, 0.0, 1e-15)); 1529 1530 auto b = exp(complex(0.0, -2.0 * PI / 3.0)); 1531 assert(isClose(b, complex(-0.5L, -0.866025403784438646763L))); 1532 1533 auto c = exp(complex(0.0, PI / 3.0)); 1534 assert(isClose(c, complex(0.5L, 0.866025403784438646763L))); 1535 1536 auto d = exp(complex(0.0, 2.0 * PI / 3.0)); 1537 assert(isClose(d, complex(-0.5L, 0.866025403784438646763L))); 1538 1539 auto e = exp(complex(0.0, PI)); 1540 assert(isClose(e, -1.0, 0.0, 1e-15)); 1541} 1542 1543/** 1544 * Calculate the natural logarithm of x. 1545 * The branch cut is along the negative axis. 1546 * Params: 1547 * x = A complex number 1548 * Returns: 1549 * The complex natural logarithm of `x` 1550 * 1551 * $(TABLE_SV 1552 * $(TR $(TH x) $(TH log(x))) 1553 * $(TR $(TD (-0, +0)) $(TD (-$(INFIN), $(PI)))) 1554 * $(TR $(TD (+0, +0)) $(TD (-$(INFIN), +0))) 1555 * $(TR $(TD (any, +$(INFIN))) $(TD (+$(INFIN), $(PI)/2))) 1556 * $(TR $(TD (any, $(NAN))) $(TD ($(NAN), $(NAN)))) 1557 * $(TR $(TD (-$(INFIN), any)) $(TD (+$(INFIN), $(PI)))) 1558 * $(TR $(TD (+$(INFIN), any)) $(TD (+$(INFIN), +0))) 1559 * $(TR $(TD (-$(INFIN), +$(INFIN))) $(TD (+$(INFIN), 3$(PI)/4))) 1560 * $(TR $(TD (+$(INFIN), +$(INFIN))) $(TD (+$(INFIN), $(PI)/4))) 1561 * $(TR $(TD ($(PLUSMN)$(INFIN), $(NAN))) $(TD (+$(INFIN), $(NAN)))) 1562 * $(TR $(TD ($(NAN), any)) $(TD ($(NAN), $(NAN)))) 1563 * $(TR $(TD ($(NAN), +$(INFIN))) $(TD (+$(INFIN), $(NAN)))) 1564 * $(TR $(TD ($(NAN), $(NAN))) $(TD ($(NAN), $(NAN)))) 1565 * ) 1566 */ 1567Complex!T log(T)(Complex!T x) @safe pure nothrow @nogc 1568{ 1569 static import std.math; 1570 1571 // Handle special cases explicitly here for better accuracy. 1572 // The order here is important, so that the correct path is chosen. 1573 if (std.math.isNaN(x.re)) 1574 { 1575 if (std.math.isInfinity(x.im)) 1576 return Complex!T(T.infinity, T.nan); 1577 else 1578 return Complex!T(T.nan, T.nan); 1579 } 1580 if (std.math.isInfinity(x.re)) 1581 { 1582 if (std.math.isNaN(x.im)) 1583 return Complex!T(T.infinity, T.nan); 1584 else if (std.math.isInfinity(x.im)) 1585 { 1586 if (std.math.signbit(x.re)) 1587 return Complex!T(T.infinity, std.math.copysign(3.0 * std.math.PI_4, x.im)); 1588 else 1589 return Complex!T(T.infinity, std.math.copysign(std.math.PI_4, x.im)); 1590 } 1591 else 1592 { 1593 if (std.math.signbit(x.re)) 1594 return Complex!T(T.infinity, std.math.copysign(std.math.PI, x.im)); 1595 else 1596 return Complex!T(T.infinity, std.math.copysign(0.0, x.im)); 1597 } 1598 } 1599 if (std.math.isNaN(x.im)) 1600 return Complex!T(T.nan, T.nan); 1601 if (std.math.isInfinity(x.im)) 1602 return Complex!T(T.infinity, std.math.copysign(std.math.PI_2, x.im)); 1603 if (x.re == 0.0 && x.im == 0.0) 1604 { 1605 if (std.math.signbit(x.re)) 1606 return Complex!T(-T.infinity, std.math.copysign(std.math.PI, x.im)); 1607 else 1608 return Complex!T(-T.infinity, std.math.copysign(0.0, x.im)); 1609 } 1610 1611 return Complex!T(std.math.log(abs(x)), arg(x)); 1612} 1613 1614/// 1615@safe pure nothrow @nogc unittest 1616{ 1617 import core.math : sqrt; 1618 import std.math.constants : PI; 1619 import std.math.operations : isClose; 1620 1621 auto a = complex(2.0, 1.0); 1622 assert(log(conj(a)) == conj(log(a))); 1623 1624 auto b = 2.0 * log10(complex(0.0, 1.0)); 1625 auto c = 4.0 * log10(complex(sqrt(2.0) / 2, sqrt(2.0) / 2)); 1626 assert(isClose(b, c, 0.0, 1e-15)); 1627 1628 assert(log(complex(-1.0L, 0.0L)) == complex(0.0L, PI)); 1629 assert(log(complex(-1.0L, -0.0L)) == complex(0.0L, -PI)); 1630} 1631 1632@safe pure nothrow @nogc unittest 1633{ 1634 import std.math.traits : isNaN, isInfinity; 1635 import std.math.constants : PI, PI_2, PI_4; 1636 1637 auto a = log(complex(-0.0L, 0.0L)); 1638 assert(a == complex(-real.infinity, PI)); 1639 auto b = log(complex(0.0L, 0.0L)); 1640 assert(b == complex(-real.infinity, +0.0L)); 1641 auto c = log(complex(1.0L, real.infinity)); 1642 assert(c == complex(real.infinity, PI_2)); 1643 auto d = log(complex(1.0L, real.nan)); 1644 assert(d.re.isNaN && d.im.isNaN); 1645 1646 auto e = log(complex(-real.infinity, 1.0L)); 1647 assert(e == complex(real.infinity, PI)); 1648 auto f = log(complex(real.infinity, 1.0L)); 1649 assert(f == complex(real.infinity, 0.0L)); 1650 auto g = log(complex(-real.infinity, real.infinity)); 1651 assert(g == complex(real.infinity, 3.0 * PI_4)); 1652 auto h = log(complex(real.infinity, real.infinity)); 1653 assert(h == complex(real.infinity, PI_4)); 1654 auto i = log(complex(real.infinity, real.nan)); 1655 assert(i.re.isInfinity && i.im.isNaN); 1656 1657 auto j = log(complex(real.nan, 1.0L)); 1658 assert(j.re.isNaN && j.im.isNaN); 1659 auto k = log(complex(real.nan, real.infinity)); 1660 assert(k.re.isInfinity && k.im.isNaN); 1661 auto l = log(complex(real.nan, real.nan)); 1662 assert(l.re.isNaN && l.im.isNaN); 1663} 1664 1665@safe pure nothrow @nogc unittest 1666{ 1667 import std.math.constants : PI; 1668 import std.math.operations : isClose; 1669 1670 auto a = log(fromPolar(1.0, PI / 6.0)); 1671 assert(isClose(a, complex(0.0L, 0.523598775598298873077L), 0.0, 1e-15)); 1672 1673 auto b = log(fromPolar(1.0, PI / 3.0)); 1674 assert(isClose(b, complex(0.0L, 1.04719755119659774615L), 0.0, 1e-15)); 1675 1676 auto c = log(fromPolar(1.0, PI / 2.0)); 1677 assert(isClose(c, complex(0.0L, 1.57079632679489661923L), 0.0, 1e-15)); 1678 1679 auto d = log(fromPolar(1.0, 2.0 * PI / 3.0)); 1680 assert(isClose(d, complex(0.0L, 2.09439510239319549230L), 0.0, 1e-15)); 1681 1682 auto e = log(fromPolar(1.0, 5.0 * PI / 6.0)); 1683 assert(isClose(e, complex(0.0L, 2.61799387799149436538L), 0.0, 1e-15)); 1684 1685 auto f = log(complex(-1.0L, 0.0L)); 1686 assert(isClose(f, complex(0.0L, PI), 0.0, 1e-15)); 1687} 1688 1689/** 1690 * Calculate the base-10 logarithm of x. 1691 * Params: 1692 * x = A complex number 1693 * Returns: 1694 * The complex base 10 logarithm of `x` 1695 */ 1696Complex!T log10(T)(Complex!T x) @safe pure nothrow @nogc 1697{ 1698 static import std.math; 1699 1700 return log(x) / Complex!T(std.math.log(10.0)); 1701} 1702 1703/// 1704@safe pure nothrow @nogc unittest 1705{ 1706 import core.math : sqrt; 1707 import std.math.constants : LN10, PI; 1708 import std.math.operations : isClose; 1709 1710 auto a = complex(2.0, 1.0); 1711 assert(log10(a) == log(a) / log(complex(10.0))); 1712 1713 auto b = log10(complex(0.0, 1.0)) * 2.0; 1714 auto c = log10(complex(sqrt(2.0) / 2, sqrt(2.0) / 2)) * 4.0; 1715 assert(isClose(b, c, 0.0, 1e-15)); 1716} 1717 1718@safe pure nothrow @nogc unittest 1719{ 1720 import std.math.constants : LN10, PI; 1721 import std.math.operations : isClose; 1722 1723 auto a = log10(fromPolar(1.0, PI / 6.0)); 1724 assert(isClose(a, complex(0.0L, 0.227396058973640224580L), 0.0, 1e-15)); 1725 1726 auto b = log10(fromPolar(1.0, PI / 3.0)); 1727 assert(isClose(b, complex(0.0L, 0.454792117947280449161L), 0.0, 1e-15)); 1728 1729 auto c = log10(fromPolar(1.0, PI / 2.0)); 1730 assert(isClose(c, complex(0.0L, 0.682188176920920673742L), 0.0, 1e-15)); 1731 1732 auto d = log10(fromPolar(1.0, 2.0 * PI / 3.0)); 1733 assert(isClose(d, complex(0.0L, 0.909584235894560898323L), 0.0, 1e-15)); 1734 1735 auto e = log10(fromPolar(1.0, 5.0 * PI / 6.0)); 1736 assert(isClose(e, complex(0.0L, 1.13698029486820112290L), 0.0, 1e-15)); 1737 1738 auto f = log10(complex(-1.0L, 0.0L)); 1739 assert(isClose(f, complex(0.0L, 1.36437635384184134748L), 0.0, 1e-15)); 1740 1741 assert(ceqrel(log10(complex(-100.0L, 0.0L)), complex(2.0L, PI / LN10)) >= real.mant_dig - 1); 1742 assert(ceqrel(log10(complex(-100.0L, -0.0L)), complex(2.0L, -PI / LN10)) >= real.mant_dig - 1); 1743} 1744 1745/** 1746 * Calculates x$(SUPERSCRIPT n). 1747 * The branch cut is on the negative axis. 1748 * Params: 1749 * x = base 1750 * n = exponent 1751 * Returns: 1752 * `x` raised to the power of `n` 1753 */ 1754Complex!T pow(T, Int)(Complex!T x, const Int n) @safe pure nothrow @nogc 1755if (isIntegral!Int) 1756{ 1757 alias UInt = Unsigned!(Unqual!Int); 1758 1759 UInt m = (n < 0) ? -cast(UInt) n : n; 1760 Complex!T y = (m % 2) ? x : Complex!T(1); 1761 1762 while (m >>= 1) 1763 { 1764 x *= x; 1765 if (m % 2) 1766 y *= x; 1767 } 1768 1769 return (n < 0) ? Complex!T(1) / y : y; 1770} 1771 1772/// 1773@safe pure nothrow @nogc unittest 1774{ 1775 import std.math.operations : isClose; 1776 1777 auto a = complex(1.0, 2.0); 1778 assert(pow(a, 2) == a * a); 1779 assert(pow(a, 3) == a * a * a); 1780 assert(pow(a, -2) == 1.0 / (a * a)); 1781 assert(isClose(pow(a, -3), 1.0 / (a * a * a))); 1782} 1783 1784/// ditto 1785Complex!T pow(T)(Complex!T x, const T n) @trusted pure nothrow @nogc 1786{ 1787 static import std.math; 1788 1789 if (x == 0.0) 1790 return Complex!T(0.0); 1791 1792 if (x.im == 0 && x.re > 0.0) 1793 return Complex!T(std.math.pow(x.re, n)); 1794 1795 Complex!T t = log(x); 1796 return fromPolar!(T, T)(std.math.exp(n * t.re), n * t.im); 1797} 1798 1799/// 1800@safe pure nothrow @nogc unittest 1801{ 1802 import std.math.operations : isClose; 1803 assert(pow(complex(0.0), 2.0) == complex(0.0)); 1804 assert(pow(complex(5.0), 2.0) == complex(25.0)); 1805 1806 auto a = pow(complex(-1.0, 0.0), 0.5); 1807 assert(isClose(a, complex(0.0, +1.0), 0.0, 1e-16)); 1808 1809 auto b = pow(complex(-1.0, -0.0), 0.5); 1810 assert(isClose(b, complex(0.0, -1.0), 0.0, 1e-16)); 1811} 1812 1813/// ditto 1814Complex!T pow(T)(Complex!T x, Complex!T y) @trusted pure nothrow @nogc 1815{ 1816 return (x == 0) ? Complex!T(0) : exp(y * log(x)); 1817} 1818 1819/// 1820@safe pure nothrow @nogc unittest 1821{ 1822 import std.math.operations : isClose; 1823 import std.math.exponential : exp; 1824 import std.math.constants : PI; 1825 auto a = complex(0.0); 1826 auto b = complex(2.0); 1827 assert(pow(a, b) == complex(0.0)); 1828 1829 auto c = complex(0.0L, 1.0L); 1830 assert(isClose(pow(c, c), exp((-PI) / 2))); 1831} 1832 1833/// ditto 1834Complex!T pow(T)(const T x, Complex!T n) @trusted pure nothrow @nogc 1835{ 1836 static import std.math; 1837 1838 return (x > 0.0) 1839 ? fromPolar!(T, T)(std.math.pow(x, n.re), n.im * std.math.log(x)) 1840 : pow(Complex!T(x), n); 1841} 1842 1843/// 1844@safe pure nothrow @nogc unittest 1845{ 1846 import std.math.operations : isClose; 1847 assert(pow(2.0, complex(0.0)) == complex(1.0)); 1848 assert(pow(2.0, complex(5.0)) == complex(32.0)); 1849 1850 auto a = pow(-2.0, complex(-1.0)); 1851 assert(isClose(a, complex(-0.5), 0.0, 1e-16)); 1852 1853 auto b = pow(-0.5, complex(-1.0)); 1854 assert(isClose(b, complex(-2.0), 0.0, 1e-15)); 1855} 1856 1857@safe pure nothrow @nogc unittest 1858{ 1859 import std.math.constants : PI; 1860 import std.math.operations : isClose; 1861 1862 auto a = pow(complex(3.0, 4.0), 2); 1863 assert(isClose(a, complex(-7.0, 24.0))); 1864 1865 auto b = pow(complex(3.0, 4.0), PI); 1866 assert(ceqrel(b, complex(-152.91512205297134, 35.547499631917738)) >= double.mant_dig - 3); 1867 1868 auto c = pow(complex(3.0, 4.0), complex(-2.0, 1.0)); 1869 assert(ceqrel(c, complex(0.015351734187477306, -0.0038407695456661503)) >= double.mant_dig - 3); 1870 1871 auto d = pow(PI, complex(2.0, -1.0)); 1872 assert(ceqrel(d, complex(4.0790296880118296, -8.9872469554541869)) >= double.mant_dig - 1); 1873 1874 auto e = complex(2.0); 1875 assert(ceqrel(pow(e, 3), exp(3 * log(e))) >= double.mant_dig - 1); 1876} 1877 1878@safe pure nothrow @nogc unittest 1879{ 1880 import std.meta : AliasSeq; 1881 import std.math : RealFormat, floatTraits; 1882 static foreach (T; AliasSeq!(float, double, real)) 1883 {{ 1884 static if (floatTraits!T.realFormat == RealFormat.ibmExtended) 1885 { 1886 /* For IBM real, epsilon is too small (since 1.0 plus any double is 1887 representable) to be able to expect results within epsilon * 100. */ 1888 } 1889 else 1890 { 1891 T eps = T.epsilon * 100; 1892 1893 T a = -1.0; 1894 T b = 0.5; 1895 Complex!T ref1 = pow(complex(a), complex(b)); 1896 Complex!T res1 = pow(a, complex(b)); 1897 Complex!T res2 = pow(complex(a), b); 1898 assert(abs(ref1 - res1) < eps); 1899 assert(abs(ref1 - res2) < eps); 1900 assert(abs(res1 - res2) < eps); 1901 1902 T c = -3.2; 1903 T d = 1.4; 1904 Complex!T ref2 = pow(complex(a), complex(b)); 1905 Complex!T res3 = pow(a, complex(b)); 1906 Complex!T res4 = pow(complex(a), b); 1907 assert(abs(ref2 - res3) < eps); 1908 assert(abs(ref2 - res4) < eps); 1909 assert(abs(res3 - res4) < eps); 1910 } 1911 }} 1912} 1913