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 }
The D Programming Language's Unofficial Docs
beta