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 $(D cfloat), $(D cdouble), $(D creal), $(D ifloat),
9 $(D idouble), and $(D ireal).
10
11 Authors: Lars Tandle Kyllingstad, Don Clugston
12 Copyright: Copyright (c) 2010, Lars T. Kyllingstad.
13 License: $(HTTP boost.org/LICENSE_1_0.txt, Boost License 1.0)
14 Source: $(PHOBOSSRC std/_complex.d)
15 */
16 module std.complex;
17
18 import std.traits;
19
20 /** Helper function that returns a _complex number with the specified
21 real and imaginary parts.
22
23 Params:
24 R = (template parameter) type of real part of complex number
25 I = (template parameter) type of imaginary part of complex number
26
27 re = real part of complex number to be constructed
28 im = (optional) imaginary part of complex number, 0 if omitted.
29
30 Returns:
31 $(D Complex) instance with real and imaginary parts set
32 to the values provided as input. If neither $(D re) nor
33 $(D im) are floating-point numbers, the return type will
34 be $(D Complex!double). Otherwise, the return type is
35 deduced using $(D std.traits.CommonType!(R, I)).
36 */
37 auto complex(R)(R re) @safe pure nothrow @nogc
38 if (is(R : double))
39 {
40 static if (isFloatingPoint!R)
41 return Complex!R(re, 0);
42 else
43 return Complex!double(re, 0);
44 }
45
46 /// ditto
47 auto complex(R, I)(R re, I im) @safe pure nothrow @nogc
48 if (is(R : double) && is(I : double))
49 {
50 static if (isFloatingPoint!R || isFloatingPoint!I)
51 return Complex!(CommonType!(R, I))(re, im);
52 else
53 return Complex!double(re, im);
54 }
55
56 ///
57 @safe pure nothrow unittest
58 {
59 auto a = complex(1.0);
60 static assert(is(typeof(a) == Complex!double));
61 assert(a.re == 1.0);
62 assert(a.im == 0.0);
63
64 auto b = complex(2.0L);
65 static assert(is(typeof(b) == Complex!real));
66 assert(b.re == 2.0L);
67 assert(b.im == 0.0L);
68
69 auto c = complex(1.0, 2.0);
70 static assert(is(typeof(c) == Complex!double));
71 assert(c.re == 1.0);
72 assert(c.im == 2.0);
73
74 auto d = complex(3.0, 4.0L);
75 static assert(is(typeof(d) == Complex!real));
76 assert(d.re == 3.0);
77 assert(d.im == 4.0L);
78
79 auto e = complex(1);
80 static assert(is(typeof(e) == Complex!double));
81 assert(e.re == 1);
82 assert(e.im == 0);
83
84 auto f = complex(1L, 2);
85 static assert(is(typeof(f) == Complex!double));
86 assert(f.re == 1L);
87 assert(f.im == 2);
88
89 auto g = complex(3, 4.0L);
90 static assert(is(typeof(g) == Complex!real));
91 assert(g.re == 3);
92 assert(g.im == 4.0L);
93 }
94
95
96 /** A complex number parametrised by a type $(D T), which must be either
97 $(D float), $(D double) or $(D real).
98 */
99 struct Complex(T)
100 if (isFloatingPoint!T)
101 {
102 import std.format : FormatSpec;
103 import std.range.primitives : isOutputRange;
104
105 /** The real part of the number. */
106 T re;
107
108 /** The imaginary part of the number. */
109 T im;
110
111 /** Converts the complex number to a string representation.
112
113 The second form of this function is usually not called directly;
114 instead, it is used via $(REF format, std,string), as shown in the examples
115 below. Supported format characters are 'e', 'f', 'g', 'a', and 's'.
116
117 See the $(MREF std, format) and $(REF format, std,string)
118 documentation for more information.
119 */
toString()120 string toString() const @safe /* TODO: pure nothrow */
121 {
122 import std.exception : assumeUnique;
123 char[] buf;
124 buf.reserve(100);
125 auto fmt = FormatSpec!char("%s");
126 toString((const(char)[] s) { buf ~= s; }, fmt);
127 static trustedAssumeUnique(T)(T t) @trusted { return assumeUnique(t); }
128 return trustedAssumeUnique(buf);
129 }
130
131 static if (is(T == double))
132 ///
133 @safe unittest
134 {
135 auto c = complex(1.2, 3.4);
136
137 // Vanilla toString formatting:
138 assert(c.toString() == "1.2+3.4i");
139
140 // Formatting with std.string.format specs: the precision and width
141 // specifiers apply to both the real and imaginary parts of the
142 // complex number.
143 import std.format : format;
144 assert(format("%.2f", c) == "1.20+3.40i");
145 assert(format("%4.1f", c) == " 1.2+ 3.4i");
146 }
147
148 /// ditto
149 void toString(Writer, Char)(scope Writer w,
150 FormatSpec!Char formatSpec) const
151 if (isOutputRange!(Writer, const(Char)[]))
152 {
153 import std.format : formatValue;
154 import std.math : signbit;
155 import std.range.primitives : put;
156 formatValue(w, re, formatSpec);
157 if (signbit(im) == 0)
158 put(w, "+");
159 formatValue(w, im, formatSpec);
160 put(w, "i");
161 }
162
163 @safe pure nothrow @nogc:
164
165 /** Construct a complex number with the specified real and
166 imaginary parts. In the case where a single argument is passed
167 that is not complex, the imaginary part of the result will be
168 zero.
169 */
170 this(R : T)(Complex!R z)
171 {
172 re = z.re;
173 im = z.im;
174 }
175
176 /// ditto
177 this(Rx : T, Ry : T)(Rx x, Ry y)
178 {
179 re = x;
180 im = y;
181 }
182
183 /// ditto
184 this(R : T)(R r)
185 {
186 re = r;
187 im = 0;
188 }
189
190 // ASSIGNMENT OPERATORS
191
192 // this = complex
193 ref Complex opAssign(R : T)(Complex!R z)
194 {
195 re = z.re;
196 im = z.im;
197 return this;
198 }
199
200 // this = numeric
201 ref Complex opAssign(R : T)(R r)
202 {
203 re = r;
204 im = 0;
205 return this;
206 }
207
208 // COMPARISON OPERATORS
209
210 // this == complex
211 bool opEquals(R : T)(Complex!R z) const
212 {
213 return re == z.re && im == z.im;
214 }
215
216 // this == numeric
217 bool opEquals(R : T)(R r) const
218 {
219 return re == r && im == 0;
220 }
221
222 // UNARY OPERATORS
223
224 // +complex
225 Complex opUnary(string op)() const
226 if (op == "+")
227 {
228 return this;
229 }
230
231 // -complex
232 Complex opUnary(string op)() const
233 if (op == "-")
234 {
235 return Complex(-re, -im);
236 }
237
238 // BINARY OPERATORS
239
240 // complex op complex
241 Complex!(CommonType!(T,R)) opBinary(string op, R)(Complex!R z) const
242 {
243 alias C = typeof(return);
244 auto w = C(this.re, this.im);
245 return w.opOpAssign!(op)(z);
246 }
247
248 // complex op numeric
249 Complex!(CommonType!(T,R)) opBinary(string op, R)(R r) const
250 if (isNumeric!R)
251 {
252 alias C = typeof(return);
253 auto w = C(this.re, this.im);
254 return w.opOpAssign!(op)(r);
255 }
256
257 // numeric + complex, numeric * complex
258 Complex!(CommonType!(T, R)) opBinaryRight(string op, R)(R r) const
259 if ((op == "+" || op == "*") && (isNumeric!R))
260 {
261 return opBinary!(op)(r);
262 }
263
264 // numeric - complex
265 Complex!(CommonType!(T, R)) opBinaryRight(string op, R)(R r) const
266 if (op == "-" && isNumeric!R)
267 {
268 return Complex(r - re, -im);
269 }
270
271 // numeric / complex
272 Complex!(CommonType!(T, R)) opBinaryRight(string op, R)(R r) const
273 if (op == "/" && isNumeric!R)
274 {
275 import std.math : fabs;
276 typeof(return) w = void;
277 if (fabs(re) < fabs(im))
278 {
279 immutable ratio = re/im;
280 immutable rdivd = r/(re*ratio + im);
281
282 w.re = rdivd*ratio;
283 w.im = -rdivd;
284 }
285 else
286 {
287 immutable ratio = im/re;
288 immutable rdivd = r/(re + im*ratio);
289
290 w.re = rdivd;
291 w.im = -rdivd*ratio;
292 }
293
294 return w;
295 }
296
297 // numeric ^^ complex
298 Complex!(CommonType!(T, R)) opBinaryRight(string op, R)(R lhs) const
299 if (op == "^^" && isNumeric!R)
300 {
301 import std.math : cos, exp, log, sin, PI;
302 Unqual!(CommonType!(T, R)) ab = void, ar = void;
303
304 if (lhs >= 0)
305 {
306 // r = lhs
307 // theta = 0
308 ab = lhs ^^ this.re;
309 ar = log(lhs) * this.im;
310 }
311 else
312 {
313 // r = -lhs
314 // theta = PI
315 ab = (-lhs) ^^ this.re * exp(-PI * this.im);
316 ar = PI * this.re + log(-lhs) * this.im;
317 }
318
319 return typeof(return)(ab * cos(ar), ab * sin(ar));
320 }
321
322 // OP-ASSIGN OPERATORS
323
324 // complex += complex, complex -= complex
325 ref Complex opOpAssign(string op, C)(C z)
326 if ((op == "+" || op == "-") && is(C R == Complex!R))
327 {
328 mixin ("re "~op~"= z.re;");
329 mixin ("im "~op~"= z.im;");
330 return this;
331 }
332
333 // complex *= complex
334 ref Complex opOpAssign(string op, C)(C z)
335 if (op == "*" && is(C R == Complex!R))
336 {
337 auto temp = re*z.re - im*z.im;
338 im = im*z.re + re*z.im;
339 re = temp;
340 return this;
341 }
342
343 // complex /= complex
344 ref Complex opOpAssign(string op, C)(C z)
345 if (op == "/" && is(C R == Complex!R))
346 {
347 import std.math : fabs;
348 if (fabs(z.re) < fabs(z.im))
349 {
350 immutable ratio = z.re/z.im;
351 immutable denom = z.re*ratio + z.im;
352
353 immutable temp = (re*ratio + im)/denom;
354 im = (im*ratio - re)/denom;
355 re = temp;
356 }
357 else
358 {
359 immutable ratio = z.im/z.re;
360 immutable denom = z.re + z.im*ratio;
361
362 immutable temp = (re + im*ratio)/denom;
363 im = (im - re*ratio)/denom;
364 re = temp;
365 }
366 return this;
367 }
368
369 // complex ^^= complex
370 ref Complex opOpAssign(string op, C)(C z)
371 if (op == "^^" && is(C R == Complex!R))
372 {
373 import std.math : exp, log, cos, sin;
374 immutable r = abs(this);
375 immutable t = arg(this);
376 immutable ab = r^^z.re * exp(-t*z.im);
377 immutable ar = t*z.re + log(r)*z.im;
378
379 re = ab*cos(ar);
380 im = ab*sin(ar);
381 return this;
382 }
383
384 // complex += numeric, complex -= numeric
385 ref Complex opOpAssign(string op, U : T)(U a)
386 if (op == "+" || op == "-")
387 {
388 mixin ("re "~op~"= a;");
389 return this;
390 }
391
392 // complex *= numeric, complex /= numeric
393 ref Complex opOpAssign(string op, U : T)(U a)
394 if (op == "*" || op == "/")
395 {
396 mixin ("re "~op~"= a;");
397 mixin ("im "~op~"= a;");
398 return this;
399 }
400
401 // complex ^^= real
402 ref Complex opOpAssign(string op, R)(R r)
403 if (op == "^^" && isFloatingPoint!R)
404 {
405 import std.math : cos, sin;
406 immutable ab = abs(this)^^r;
407 immutable ar = arg(this)*r;
408 re = ab*cos(ar);
409 im = ab*sin(ar);
410 return this;
411 }
412
413 // complex ^^= int
414 ref Complex opOpAssign(string op, U)(U i)
415 if (op == "^^" && isIntegral!U)
416 {
417 switch (i)
418 {
419 case 0:
420 re = 1.0;
421 im = 0.0;
422 break;
423 case 1:
424 // identity; do nothing
425 break;
426 case 2:
427 this *= this;
428 break;
429 case 3:
430 auto z = this;
431 this *= z;
432 this *= z;
433 break;
434 default:
435 this ^^= cast(real) i;
436 }
437 return this;
438 }
439 }
440
441 @safe pure nothrow unittest
442 {
443 import std.complex;
444 import std.math;
445
446 enum EPS = double.epsilon;
447 auto c1 = complex(1.0, 1.0);
448
449 // Check unary operations.
450 auto c2 = Complex!double(0.5, 2.0);
451
452 assert(c2 == +c2);
453
454 assert((-c2).re == -(c2.re));
455 assert((-c2).im == -(c2.im));
456 assert(c2 == -(-c2));
457
458 // Check complex-complex operations.
459 auto cpc = c1 + c2;
460 assert(cpc.re == c1.re + c2.re);
461 assert(cpc.im == c1.im + c2.im);
462
463 auto cmc = c1 - c2;
464 assert(cmc.re == c1.re - c2.re);
465 assert(cmc.im == c1.im - c2.im);
466
467 auto ctc = c1 * c2;
468 assert(approxEqual(abs(ctc), abs(c1)*abs(c2), EPS));
469 assert(approxEqual(arg(ctc), arg(c1)+arg(c2), EPS));
470
471 auto cdc = c1 / c2;
472 assert(approxEqual(abs(cdc), abs(c1)/abs(c2), EPS));
473 assert(approxEqual(arg(cdc), arg(c1)-arg(c2), EPS));
474
475 auto cec = c1^^c2;
476 assert(approxEqual(cec.re, 0.11524131979943839881, EPS));
477 assert(approxEqual(cec.im, 0.21870790452746026696, EPS));
478
479 // Check complex-real operations.
480 double a = 123.456;
481
482 auto cpr = c1 + a;
483 assert(cpr.re == c1.re + a);
484 assert(cpr.im == c1.im);
485
486 auto cmr = c1 - a;
487 assert(cmr.re == c1.re - a);
488 assert(cmr.im == c1.im);
489
490 auto ctr = c1 * a;
491 assert(ctr.re == c1.re*a);
492 assert(ctr.im == c1.im*a);
493
494 auto cdr = c1 / a;
495 assert(approxEqual(abs(cdr), abs(c1)/a, EPS));
496 assert(approxEqual(arg(cdr), arg(c1), EPS));
497
498 auto cer = c1^^3.0;
499 assert(approxEqual(abs(cer), abs(c1)^^3, EPS));
500 assert(approxEqual(arg(cer), arg(c1)*3, EPS));
501
502 auto rpc = a + c1;
503 assert(rpc == cpr);
504
505 auto rmc = a - c1;
506 assert(rmc.re == a-c1.re);
507 assert(rmc.im == -c1.im);
508
509 auto rtc = a * c1;
510 assert(rtc == ctr);
511
512 auto rdc = a / c1;
513 assert(approxEqual(abs(rdc), a/abs(c1), EPS));
514 assert(approxEqual(arg(rdc), -arg(c1), EPS));
515
516 rdc = a / c2;
517 assert(approxEqual(abs(rdc), a/abs(c2), EPS));
518 assert(approxEqual(arg(rdc), -arg(c2), EPS));
519
520 auto rec1a = 1.0 ^^ c1;
521 assert(rec1a.re == 1.0);
522 assert(rec1a.im == 0.0);
523
524 auto rec2a = 1.0 ^^ c2;
525 assert(rec2a.re == 1.0);
526 assert(rec2a.im == 0.0);
527
528 auto rec1b = (-1.0) ^^ c1;
529 assert(approxEqual(abs(rec1b), std.math.exp(-PI * c1.im), EPS));
530 auto arg1b = arg(rec1b);
531 /* The argument _should_ be PI, but floating-point rounding error
532 * means that in fact the imaginary part is very slightly negative.
533 */
534 assert(approxEqual(arg1b, PI, EPS) || approxEqual(arg1b, -PI, EPS));
535
536 auto rec2b = (-1.0) ^^ c2;
537 assert(approxEqual(abs(rec2b), std.math.exp(-2 * PI), EPS));
538 assert(approxEqual(arg(rec2b), PI_2, EPS));
539
540 auto rec3a = 0.79 ^^ complex(6.8, 5.7);
541 auto rec3b = complex(0.79, 0.0) ^^ complex(6.8, 5.7);
542 assert(approxEqual(rec3a.re, rec3b.re, EPS));
543 assert(approxEqual(rec3a.im, rec3b.im, EPS));
544
545 auto rec4a = (-0.79) ^^ complex(6.8, 5.7);
546 auto rec4b = complex(-0.79, 0.0) ^^ complex(6.8, 5.7);
547 assert(approxEqual(rec4a.re, rec4b.re, EPS));
548 assert(approxEqual(rec4a.im, rec4b.im, EPS));
549
550 auto rer = a ^^ complex(2.0, 0.0);
551 auto rcheck = a ^^ 2.0;
552 static assert(is(typeof(rcheck) == double));
553 assert(feqrel(rer.re, rcheck) == double.mant_dig);
554 assert(isIdentical(rer.re, rcheck));
555 assert(rer.im == 0.0);
556
557 auto rer2 = (-a) ^^ complex(2.0, 0.0);
558 rcheck = (-a) ^^ 2.0;
559 assert(feqrel(rer2.re, rcheck) == double.mant_dig);
560 assert(isIdentical(rer2.re, rcheck));
561 assert(approxEqual(rer2.im, 0.0, EPS));
562
563 auto rer3 = (-a) ^^ complex(-2.0, 0.0);
564 rcheck = (-a) ^^ (-2.0);
565 assert(feqrel(rer3.re, rcheck) == double.mant_dig);
566 assert(isIdentical(rer3.re, rcheck));
567 assert(approxEqual(rer3.im, 0.0, EPS));
568
569 auto rer4 = a ^^ complex(-2.0, 0.0);
570 rcheck = a ^^ (-2.0);
571 assert(feqrel(rer4.re, rcheck) == double.mant_dig);
572 assert(isIdentical(rer4.re, rcheck));
573 assert(rer4.im == 0.0);
574
575 // Check Complex-int operations.
576 foreach (i; 0 .. 6)
577 {
578 auto cei = c1^^i;
579 assert(approxEqual(abs(cei), abs(c1)^^i, EPS));
580 // Use cos() here to deal with arguments that go outside
581 // the (-pi,pi] interval (only an issue for i>3).
582 assert(approxEqual(std.math.cos(arg(cei)), std.math.cos(arg(c1)*i), EPS));
583 }
584
585 // Check operations between different complex types.
586 auto cf = Complex!float(1.0, 1.0);
587 auto cr = Complex!real(1.0, 1.0);
588 auto c1pcf = c1 + cf;
589 auto c1pcr = c1 + cr;
590 static assert(is(typeof(c1pcf) == Complex!double));
591 static assert(is(typeof(c1pcr) == Complex!real));
592 assert(c1pcf.re == c1pcr.re);
593 assert(c1pcf.im == c1pcr.im);
594
595 auto c1c = c1;
596 auto c2c = c2;
597
598 c1c /= c1;
599 assert(approxEqual(c1c.re, 1.0, EPS));
600 assert(approxEqual(c1c.im, 0.0, EPS));
601
602 c1c = c1;
603 c1c /= c2;
604 assert(approxEqual(c1c.re, 0.588235, EPS));
605 assert(approxEqual(c1c.im, -0.352941, EPS));
606
607 c2c /= c1;
608 assert(approxEqual(c2c.re, 1.25, EPS));
609 assert(approxEqual(c2c.im, 0.75, EPS));
610
611 c2c = c2;
612 c2c /= c2;
613 assert(approxEqual(c2c.re, 1.0, EPS));
614 assert(approxEqual(c2c.im, 0.0, EPS));
615 }
616
617 @safe pure nothrow unittest
618 {
619 // Initialization
620 Complex!double a = 1;
621 assert(a.re == 1 && a.im == 0);
622 Complex!double b = 1.0;
623 assert(b.re == 1.0 && b.im == 0);
624 Complex!double c = Complex!real(1.0, 2);
625 assert(c.re == 1.0 && c.im == 2);
626 }
627
628 @safe pure nothrow unittest
629 {
630 // Assignments and comparisons
631 Complex!double z;
632
633 z = 1;
634 assert(z == 1);
635 assert(z.re == 1.0 && z.im == 0.0);
636
637 z = 2.0;
638 assert(z == 2.0);
639 assert(z.re == 2.0 && z.im == 0.0);
640
641 z = 1.0L;
642 assert(z == 1.0L);
643 assert(z.re == 1.0 && z.im == 0.0);
644
645 auto w = Complex!real(1.0, 1.0);
646 z = w;
647 assert(z == w);
648 assert(z.re == 1.0 && z.im == 1.0);
649
650 auto c = Complex!float(2.0, 2.0);
651 z = c;
652 assert(z == c);
653 assert(z.re == 2.0 && z.im == 2.0);
654 }
655
656
657 /* Makes Complex!(Complex!T) fold to Complex!T.
658
659 The rationale for this is that just like the real line is a
660 subspace of the complex plane, the complex plane is a subspace
661 of itself. Example of usage:
662 ---
663 Complex!T addI(T)(T x)
664 {
665 return x + Complex!T(0.0, 1.0);
666 }
667 ---
668 The above will work if T is both real and complex.
669 */
670 template Complex(T)
671 if (is(T R == Complex!R))
672 {
673 alias Complex = T;
674 }
675
676 @safe pure nothrow unittest
677 {
678 static assert(is(Complex!(Complex!real) == Complex!real));
679
680 Complex!T addI(T)(T x)
681 {
682 return x + Complex!T(0.0, 1.0);
683 }
684
685 auto z1 = addI(1.0);
686 assert(z1.re == 1.0 && z1.im == 1.0);
687
688 enum one = Complex!double(1.0, 0.0);
689 auto z2 = addI(one);
690 assert(z1 == z2);
691 }
692
693
694 /**
695 Params: z = A complex number.
696 Returns: The absolute value (or modulus) of `z`.
697 */
abs(T)698 T abs(T)(Complex!T z) @safe pure nothrow @nogc
699 {
700 import std.math : hypot;
701 return hypot(z.re, z.im);
702 }
703
704 ///
705 @safe pure nothrow unittest
706 {
707 static import std.math;
708 assert(abs(complex(1.0)) == 1.0);
709 assert(abs(complex(0.0, 1.0)) == 1.0);
710 assert(abs(complex(1.0L, -2.0L)) == std.math.sqrt(5.0L));
711 }
712
713
714 /++
715 Params:
716 z = A complex number.
717 x = A real number.
718 Returns: The squared modulus of `z`.
719 For genericity, if called on a real number, returns its square.
720 +/
sqAbs(T)721 T sqAbs(T)(Complex!T z) @safe pure nothrow @nogc
722 {
723 return z.re*z.re + z.im*z.im;
724 }
725
726 ///
727 @safe pure nothrow unittest
728 {
729 import std.math;
730 assert(sqAbs(complex(0.0)) == 0.0);
731 assert(sqAbs(complex(1.0)) == 1.0);
732 assert(sqAbs(complex(0.0, 1.0)) == 1.0);
733 assert(approxEqual(sqAbs(complex(1.0L, -2.0L)), 5.0L));
734 assert(approxEqual(sqAbs(complex(-3.0L, 1.0L)), 10.0L));
735 assert(approxEqual(sqAbs(complex(1.0f,-1.0f)), 2.0f));
736 }
737
738 /// ditto
739 T sqAbs(T)(T x) @safe pure nothrow @nogc
740 if (isFloatingPoint!T)
741 {
742 return x*x;
743 }
744
745 @safe pure nothrow unittest
746 {
747 import std.math;
748 assert(sqAbs(0.0) == 0.0);
749 assert(sqAbs(-1.0) == 1.0);
750 assert(approxEqual(sqAbs(-3.0L), 9.0L));
751 assert(approxEqual(sqAbs(-5.0f), 25.0f));
752 }
753
754
755 /**
756 Params: z = A complex number.
757 Returns: The argument (or phase) of `z`.
758 */
arg(T)759 T arg(T)(Complex!T z) @safe pure nothrow @nogc
760 {
761 import std.math : atan2;
762 return atan2(z.im, z.re);
763 }
764
765 ///
766 @safe pure nothrow unittest
767 {
768 import std.math;
769 assert(arg(complex(1.0)) == 0.0);
770 assert(arg(complex(0.0L, 1.0L)) == PI_2);
771 assert(arg(complex(1.0L, 1.0L)) == PI_4);
772 }
773
774
775 /**
776 Params: z = A complex number.
777 Returns: The complex conjugate of `z`.
778 */
779 Complex!T conj(T)(Complex!T z) @safe pure nothrow @nogc
780 {
781 return Complex!T(z.re, -z.im);
782 }
783
784 ///
785 @safe pure nothrow unittest
786 {
787 assert(conj(complex(1.0)) == complex(1.0));
788 assert(conj(complex(1.0, 2.0)) == complex(1.0, -2.0));
789 }
790
791
792 /**
793 Constructs a complex number given its absolute value and argument.
794 Params:
795 modulus = The modulus
796 argument = The argument
797 Returns: The complex number with the given modulus and argument.
798 */
799 Complex!(CommonType!(T, U)) fromPolar(T, U)(T modulus, U argument)
800 @safe pure nothrow @nogc
801 {
802 import std.math : sin, cos;
803 return Complex!(CommonType!(T,U))
804 (modulus*cos(argument), modulus*sin(argument));
805 }
806
807 ///
808 @safe pure nothrow unittest
809 {
810 import std.math;
811 auto z = fromPolar(std.math.sqrt(2.0), PI_4);
812 assert(approxEqual(z.re, 1.0L, real.epsilon));
813 assert(approxEqual(z.im, 1.0L, real.epsilon));
814 }
815
816
817 /**
818 Trigonometric functions on complex numbers.
819
820 Params: z = A complex number.
821 Returns: The sine and cosine of `z`, respectively.
822 */
823 Complex!T sin(T)(Complex!T z) @safe pure nothrow @nogc
824 {
825 import std.math : expi, coshisinh;
826 auto cs = expi(z.re);
827 auto csh = coshisinh(z.im);
828 return typeof(return)(cs.im * csh.re, cs.re * csh.im);
829 }
830
831 ///
832 @safe pure nothrow unittest
833 {
834 static import std.math;
835 import std.math : feqrel;
836 assert(sin(complex(0.0)) == 0.0);
837 assert(sin(complex(2.0, 0)) == std.math.sin(2.0));
838 auto c1 = sin(complex(2.0L, 0));
839 auto c2 = complex(std.math.sin(2.0L), 0);
840 assert(feqrel(c1.re, c2.re) >= real.mant_dig - 1 &&
841 feqrel(c1.im, c2.im) >= real.mant_dig - 1);
842 }
843
844
845 /// ditto
846 Complex!T cos(T)(Complex!T z) @safe pure nothrow @nogc
847 {
848 import std.math : expi, coshisinh;
849 auto cs = expi(z.re);
850 auto csh = coshisinh(z.im);
851 return typeof(return)(cs.re * csh.re, - cs.im * csh.im);
852 }
853
854 ///
855 @safe pure nothrow unittest
856 {
857 static import std.math;
858 import std.math : feqrel;
859 assert(cos(complex(0.0)) == 1.0);
860 assert(cos(complex(1.3)) == std.math.cos(1.3));
861 auto c1 = cos(complex(0, 5.2L));
862 auto c2 = complex(std.math.cosh(5.2L), 0.0L);
863 assert(feqrel(c1.re, c2.re) >= real.mant_dig - 1 &&
864 feqrel(c1.im, c2.im) >= real.mant_dig - 1);
865 auto c3 = cos(complex(1.3L));
866 auto c4 = complex(std.math.cos(1.3L), 0.0L);
867 assert(feqrel(c3.re, c4.re) >= real.mant_dig - 1 &&
868 feqrel(c3.im, c4.im) >= real.mant_dig - 1);
869 }
870
871 /**
872 Params: y = A real number.
873 Returns: The value of cos(y) + i sin(y).
874
875 Note:
876 $(D expi) is included here for convenience and for easy migration of code
877 that uses $(REF _expi, std,math). Unlike $(REF _expi, std,math), which uses the
878 x87 $(I fsincos) instruction when possible, this function is no faster
879 than calculating cos(y) and sin(y) separately.
880 */
881 Complex!real expi(real y) @trusted pure nothrow @nogc
882 {
883 import std.math : cos, sin;
884 return Complex!real(cos(y), sin(y));
885 }
886
887 ///
888 @safe pure nothrow unittest
889 {
890 static import std.math;
891
892 assert(expi(1.3e5L) == complex(std.math.cos(1.3e5L), std.math.sin(1.3e5L)));
893 assert(expi(0.0L) == 1.0L);
894 auto z1 = expi(1.234);
895 auto z2 = std.math.expi(1.234);
896 assert(z1.re == z2.re && z1.im == z2.im);
897 }
898
899
900 /**
901 Params: z = A complex number.
902 Returns: The square root of `z`.
903 */
904 Complex!T sqrt(T)(Complex!T z) @safe pure nothrow @nogc
905 {
906 static import std.math;
907 typeof(return) c;
908 real x,y,w,r;
909
910 if (z == 0)
911 {
912 c = typeof(return)(0, 0);
913 }
914 else
915 {
916 real z_re = z.re;
917 real z_im = z.im;
918
919 x = std.math.fabs(z_re);
920 y = std.math.fabs(z_im);
921 if (x >= y)
922 {
923 r = y / x;
924 w = std.math.sqrt(x)
925 * std.math.sqrt(0.5 * (1 + std.math.sqrt(1 + r * r)));
926 }
927 else
928 {
929 r = x / y;
930 w = std.math.sqrt(y)
931 * std.math.sqrt(0.5 * (r + std.math.sqrt(1 + r * r)));
932 }
933
934 if (z_re >= 0)
935 {
936 c = typeof(return)(w, z_im / (w + w));
937 }
938 else
939 {
940 if (z_im < 0)
941 w = -w;
942 c = typeof(return)(z_im / (w + w), w);
943 }
944 }
945 return c;
946 }
947
948 ///
949 @safe pure nothrow unittest
950 {
951 static import std.math;
952 assert(sqrt(complex(0.0)) == 0.0);
953 assert(sqrt(complex(1.0L, 0)) == std.math.sqrt(1.0L));
954 assert(sqrt(complex(-1.0L, 0)) == complex(0, 1.0L));
955 }
956
957 @safe pure nothrow unittest
958 {
959 import std.math : approxEqual;
960
961 auto c1 = complex(1.0, 1.0);
962 auto c2 = Complex!double(0.5, 2.0);
963
964 auto c1s = sqrt(c1);
965 assert(approxEqual(c1s.re, 1.09868411));
966 assert(approxEqual(c1s.im, 0.45508986));
967
968 auto c2s = sqrt(c2);
969 assert(approxEqual(c2s.re, 1.1317134));
970 assert(approxEqual(c2s.im, 0.8836155));
971 }
972
973 // Issue 10881: support %f formatting of complex numbers
974 @safe unittest
975 {
976 import std.format : format;
977
978 auto x = complex(1.2, 3.4);
979 assert(format("%.2f", x) == "1.20+3.40i");
980
981 auto y = complex(1.2, -3.4);
982 assert(format("%.2f", y) == "1.20-3.40i");
983 }
984
985 @safe unittest
986 {
987 // Test wide string formatting
988 import std.format;
wformat(T)989 wstring wformat(T)(string format, Complex!T c)
990 {
991 import std.array : appender;
992 auto w = appender!wstring();
993 auto n = formattedWrite(w, format, c);
994 return w.data;
995 }
996
997 auto x = complex(1.2, 3.4);
998 assert(wformat("%.2f", x) == "1.20+3.40i"w);
999 }
1000
1001 @safe unittest
1002 {
1003 // Test ease of use (vanilla toString() should be supported)
1004 assert(complex(1.2, 3.4).toString() == "1.2+3.4i");
1005 }
1006