diff options
Diffstat (limited to 'libphobos/src/std/complex.d')
-rw-r--r-- | libphobos/src/std/complex.d | 1235 |
1 files changed, 1067 insertions, 168 deletions
diff --git a/libphobos/src/std/complex.d b/libphobos/src/std/complex.d index 8e488db4162..756d1ca94bb 100644 --- a/libphobos/src/std/complex.d +++ b/libphobos/src/std/complex.d @@ -1,23 +1,32 @@ // Written in the D programming language. /** This module contains the $(LREF Complex) type, which is used to represent - _complex numbers, along with related mathematical operations and functions. + complex numbers, along with related mathematical operations and functions. $(LREF Complex) will eventually $(DDLINK deprecate, Deprecated Features, replace) - the built-in types $(D cfloat), $(D cdouble), $(D creal), $(D ifloat), - $(D idouble), and $(D ireal). + the built-in types `cfloat`, `cdouble`, `creal`, `ifloat`, + `idouble`, and `ireal`. + + Macros: + TABLE_SV = <table border="1" cellpadding="4" cellspacing="0"> + <caption>Special Values</caption> + $0</table> + PLUSMN = ± + NAN = $(RED NAN) + INFIN = ∞ + PI = π Authors: Lars Tandle Kyllingstad, Don Clugston Copyright: Copyright (c) 2010, Lars T. Kyllingstad. License: $(HTTP boost.org/LICENSE_1_0.txt, Boost License 1.0) - Source: $(PHOBOSSRC std/_complex.d) + Source: $(PHOBOSSRC std/complex.d) */ module std.complex; import std.traits; -/** Helper function that returns a _complex number with the specified +/** Helper function that returns a complex number with the specified real and imaginary parts. Params: @@ -28,13 +37,13 @@ import std.traits; im = (optional) imaginary part of complex number, 0 if omitted. Returns: - $(D Complex) instance with real and imaginary parts set - to the values provided as input. If neither $(D re) nor - $(D im) are floating-point numbers, the return type will - be $(D Complex!double). Otherwise, the return type is + `Complex` instance with real and imaginary parts set + to the values provided as input. If neither `re` nor + `im` are floating-point numbers, the return type will + be `Complex!double`. Otherwise, the return type is deduced using $(D std.traits.CommonType!(R, I)). */ -auto complex(R)(R re) @safe pure nothrow @nogc +auto complex(R)(const R re) @safe pure nothrow @nogc if (is(R : double)) { static if (isFloatingPoint!R) @@ -44,7 +53,7 @@ if (is(R : double)) } /// ditto -auto complex(R, I)(R re, I im) @safe pure nothrow @nogc +auto complex(R, I)(const R re, const I im) @safe pure nothrow @nogc if (is(R : double) && is(I : double)) { static if (isFloatingPoint!R || isFloatingPoint!I) @@ -93,13 +102,13 @@ if (is(R : double) && is(I : double)) } -/** A complex number parametrised by a type $(D T), which must be either - $(D float), $(D double) or $(D real). +/** A complex number parametrised by a type `T`, which must be either + `float`, `double` or `real`. */ struct Complex(T) if (isFloatingPoint!T) { - import std.format : FormatSpec; + import std.format.spec : FormatSpec; import std.range.primitives : isOutputRange; /** The real part of the number. */ @@ -146,12 +155,11 @@ if (isFloatingPoint!T) } /// ditto - void toString(Writer, Char)(scope Writer w, - FormatSpec!Char formatSpec) const + void toString(Writer, Char)(scope Writer w, scope const ref FormatSpec!Char formatSpec) const if (isOutputRange!(Writer, const(Char)[])) { - import std.format : formatValue; - import std.math : signbit; + import std.format.write : formatValue; + import std.math.traits : signbit; import std.range.primitives : put; formatValue(w, re, formatSpec); if (signbit(im) == 0) @@ -174,14 +182,14 @@ if (isFloatingPoint!T) } /// ditto - this(Rx : T, Ry : T)(Rx x, Ry y) + this(Rx : T, Ry : T)(const Rx x, const Ry y) { re = x; im = y; } /// ditto - this(R : T)(R r) + this(R : T)(const R r) { re = r; im = 0; @@ -198,7 +206,7 @@ if (isFloatingPoint!T) } // this = numeric - ref Complex opAssign(R : T)(R r) + ref Complex opAssign(R : T)(const R r) { re = r; im = 0; @@ -214,7 +222,7 @@ if (isFloatingPoint!T) } // this == numeric - bool opEquals(R : T)(R r) const + bool opEquals(R : T)(const R r) const { return re == r && im == 0; } @@ -246,7 +254,7 @@ if (isFloatingPoint!T) } // complex op numeric - Complex!(CommonType!(T,R)) opBinary(string op, R)(R r) const + Complex!(CommonType!(T,R)) opBinary(string op, R)(const R r) const if (isNumeric!R) { alias C = typeof(return); @@ -255,50 +263,68 @@ if (isFloatingPoint!T) } // numeric + complex, numeric * complex - Complex!(CommonType!(T, R)) opBinaryRight(string op, R)(R r) const + Complex!(CommonType!(T, R)) opBinaryRight(string op, R)(const R r) const if ((op == "+" || op == "*") && (isNumeric!R)) { return opBinary!(op)(r); } // numeric - complex - Complex!(CommonType!(T, R)) opBinaryRight(string op, R)(R r) const + Complex!(CommonType!(T, R)) opBinaryRight(string op, R)(const R r) const if (op == "-" && isNumeric!R) { return Complex(r - re, -im); } // numeric / complex - Complex!(CommonType!(T, R)) opBinaryRight(string op, R)(R r) const + Complex!(CommonType!(T, R)) opBinaryRight(string op, R)(const R r) const if (op == "/" && isNumeric!R) { - import std.math : fabs; - typeof(return) w = void; - if (fabs(re) < fabs(im)) + version (FastMath) { - immutable ratio = re/im; - immutable rdivd = r/(re*ratio + im); - - w.re = rdivd*ratio; - w.im = -rdivd; + // Compute norm(this) + immutable norm = re * re + im * im; + // Compute r * conj(this) + immutable prod_re = r * re; + immutable prod_im = r * -im; + // Divide the product by the norm + typeof(return) w = void; + w.re = prod_re / norm; + w.im = prod_im / norm; + return w; } else { - immutable ratio = im/re; - immutable rdivd = r/(re + im*ratio); - - w.re = rdivd; - w.im = -rdivd*ratio; + import core.math : fabs; + typeof(return) w = void; + if (fabs(re) < fabs(im)) + { + immutable ratio = re/im; + immutable rdivd = r/(re*ratio + im); + + w.re = rdivd*ratio; + w.im = -rdivd; + } + else + { + immutable ratio = im/re; + immutable rdivd = r/(re + im*ratio); + + w.re = rdivd; + w.im = -rdivd*ratio; + } + + return w; } - - return w; } // numeric ^^ complex - Complex!(CommonType!(T, R)) opBinaryRight(string op, R)(R lhs) const + Complex!(CommonType!(T, R)) opBinaryRight(string op, R)(const R lhs) const if (op == "^^" && isNumeric!R) { - import std.math : cos, exp, log, sin, PI; + import core.math : cos, sin; + import std.math.exponential : exp, log; + import std.math.constants : PI; Unqual!(CommonType!(T, R)) ab = void, ar = void; if (lhs >= 0) @@ -322,7 +348,7 @@ if (isFloatingPoint!T) // OP-ASSIGN OPERATORS // complex += complex, complex -= complex - ref Complex opOpAssign(string op, C)(C z) + ref Complex opOpAssign(string op, C)(const C z) if ((op == "+" || op == "-") && is(C R == Complex!R)) { mixin ("re "~op~"= z.re;"); @@ -331,7 +357,7 @@ if (isFloatingPoint!T) } // complex *= complex - ref Complex opOpAssign(string op, C)(C z) + ref Complex opOpAssign(string op, C)(const C z) if (op == "*" && is(C R == Complex!R)) { auto temp = re*z.re - im*z.im; @@ -341,36 +367,52 @@ if (isFloatingPoint!T) } // complex /= complex - ref Complex opOpAssign(string op, C)(C z) + ref Complex opOpAssign(string op, C)(const C z) if (op == "/" && is(C R == Complex!R)) { - import std.math : fabs; - if (fabs(z.re) < fabs(z.im)) + version (FastMath) { - immutable ratio = z.re/z.im; - immutable denom = z.re*ratio + z.im; - - immutable temp = (re*ratio + im)/denom; - im = (im*ratio - re)/denom; - re = temp; + // Compute norm(z) + immutable norm = z.re * z.re + z.im * z.im; + // Compute this * conj(z) + immutable prod_re = re * z.re - im * -z.im; + immutable prod_im = im * z.re + re * -z.im; + // Divide the product by the norm + re = prod_re / norm; + im = prod_im / norm; + return this; } else { - immutable ratio = z.im/z.re; - immutable denom = z.re + z.im*ratio; - - immutable temp = (re + im*ratio)/denom; - im = (im - re*ratio)/denom; - re = temp; + import core.math : fabs; + if (fabs(z.re) < fabs(z.im)) + { + immutable ratio = z.re/z.im; + immutable denom = z.re*ratio + z.im; + + immutable temp = (re*ratio + im)/denom; + im = (im*ratio - re)/denom; + re = temp; + } + else + { + immutable ratio = z.im/z.re; + immutable denom = z.re + z.im*ratio; + + immutable temp = (re + im*ratio)/denom; + im = (im - re*ratio)/denom; + re = temp; + } + return this; } - return this; } // complex ^^= complex - ref Complex opOpAssign(string op, C)(C z) + ref Complex opOpAssign(string op, C)(const C z) if (op == "^^" && is(C R == Complex!R)) { - import std.math : exp, log, cos, sin; + import core.math : cos, sin; + import std.math.exponential : exp, log; immutable r = abs(this); immutable t = arg(this); immutable ab = r^^z.re * exp(-t*z.im); @@ -382,7 +424,7 @@ if (isFloatingPoint!T) } // complex += numeric, complex -= numeric - ref Complex opOpAssign(string op, U : T)(U a) + ref Complex opOpAssign(string op, U : T)(const U a) if (op == "+" || op == "-") { mixin ("re "~op~"= a;"); @@ -390,7 +432,7 @@ if (isFloatingPoint!T) } // complex *= numeric, complex /= numeric - ref Complex opOpAssign(string op, U : T)(U a) + ref Complex opOpAssign(string op, U : T)(const U a) if (op == "*" || op == "/") { mixin ("re "~op~"= a;"); @@ -399,10 +441,10 @@ if (isFloatingPoint!T) } // complex ^^= real - ref Complex opOpAssign(string op, R)(R r) + ref Complex opOpAssign(string op, R)(const R r) if (op == "^^" && isFloatingPoint!R) { - import std.math : cos, sin; + import core.math : cos, sin; immutable ab = abs(this)^^r; immutable ar = arg(this)*r; re = ab*cos(ar); @@ -411,7 +453,7 @@ if (isFloatingPoint!T) } // complex ^^= int - ref Complex opOpAssign(string op, U)(U i) + ref Complex opOpAssign(string op, U)(const U i) if (op == "^^" && isIntegral!U) { switch (i) @@ -441,6 +483,7 @@ if (isFloatingPoint!T) @safe pure nothrow unittest { import std.complex; + static import core.math; import std.math; enum EPS = double.epsilon; @@ -465,16 +508,16 @@ if (isFloatingPoint!T) assert(cmc.im == c1.im - c2.im); auto ctc = c1 * c2; - assert(approxEqual(abs(ctc), abs(c1)*abs(c2), EPS)); - assert(approxEqual(arg(ctc), arg(c1)+arg(c2), EPS)); + assert(isClose(abs(ctc), abs(c1)*abs(c2), EPS)); + assert(isClose(arg(ctc), arg(c1)+arg(c2), EPS)); auto cdc = c1 / c2; - assert(approxEqual(abs(cdc), abs(c1)/abs(c2), EPS)); - assert(approxEqual(arg(cdc), arg(c1)-arg(c2), EPS)); + assert(isClose(abs(cdc), abs(c1)/abs(c2), EPS)); + assert(isClose(arg(cdc), arg(c1)-arg(c2), EPS)); auto cec = c1^^c2; - assert(approxEqual(cec.re, 0.11524131979943839881, EPS)); - assert(approxEqual(cec.im, 0.21870790452746026696, EPS)); + assert(isClose(cec.re, 0.1152413197994, 1e-12)); + assert(isClose(cec.im, 0.2187079045274, 1e-12)); // Check complex-real operations. double a = 123.456; @@ -492,12 +535,12 @@ if (isFloatingPoint!T) assert(ctr.im == c1.im*a); auto cdr = c1 / a; - assert(approxEqual(abs(cdr), abs(c1)/a, EPS)); - assert(approxEqual(arg(cdr), arg(c1), EPS)); + assert(isClose(abs(cdr), abs(c1)/a, EPS)); + assert(isClose(arg(cdr), arg(c1), EPS)); auto cer = c1^^3.0; - assert(approxEqual(abs(cer), abs(c1)^^3, EPS)); - assert(approxEqual(arg(cer), arg(c1)*3, EPS)); + assert(isClose(abs(cer), abs(c1)^^3, EPS)); + assert(isClose(arg(cer), arg(c1)*3, EPS)); auto rpc = a + c1; assert(rpc == cpr); @@ -510,12 +553,12 @@ if (isFloatingPoint!T) assert(rtc == ctr); auto rdc = a / c1; - assert(approxEqual(abs(rdc), a/abs(c1), EPS)); - assert(approxEqual(arg(rdc), -arg(c1), EPS)); + assert(isClose(abs(rdc), a/abs(c1), EPS)); + assert(isClose(arg(rdc), -arg(c1), EPS)); rdc = a / c2; - assert(approxEqual(abs(rdc), a/abs(c2), EPS)); - assert(approxEqual(arg(rdc), -arg(c2), EPS)); + assert(isClose(abs(rdc), a/abs(c2), EPS)); + assert(isClose(arg(rdc), -arg(c2), EPS)); auto rec1a = 1.0 ^^ c1; assert(rec1a.re == 1.0); @@ -526,26 +569,26 @@ if (isFloatingPoint!T) assert(rec2a.im == 0.0); auto rec1b = (-1.0) ^^ c1; - assert(approxEqual(abs(rec1b), std.math.exp(-PI * c1.im), EPS)); + assert(isClose(abs(rec1b), std.math.exp(-PI * c1.im), EPS)); auto arg1b = arg(rec1b); /* The argument _should_ be PI, but floating-point rounding error * means that in fact the imaginary part is very slightly negative. */ - assert(approxEqual(arg1b, PI, EPS) || approxEqual(arg1b, -PI, EPS)); + assert(isClose(arg1b, PI, EPS) || isClose(arg1b, -PI, EPS)); auto rec2b = (-1.0) ^^ c2; - assert(approxEqual(abs(rec2b), std.math.exp(-2 * PI), EPS)); - assert(approxEqual(arg(rec2b), PI_2, EPS)); + assert(isClose(abs(rec2b), std.math.exp(-2 * PI), EPS)); + assert(isClose(arg(rec2b), PI_2, EPS)); auto rec3a = 0.79 ^^ complex(6.8, 5.7); auto rec3b = complex(0.79, 0.0) ^^ complex(6.8, 5.7); - assert(approxEqual(rec3a.re, rec3b.re, EPS)); - assert(approxEqual(rec3a.im, rec3b.im, EPS)); + assert(isClose(rec3a.re, rec3b.re, 1e-14)); + assert(isClose(rec3a.im, rec3b.im, 1e-14)); auto rec4a = (-0.79) ^^ complex(6.8, 5.7); auto rec4b = complex(-0.79, 0.0) ^^ complex(6.8, 5.7); - assert(approxEqual(rec4a.re, rec4b.re, EPS)); - assert(approxEqual(rec4a.im, rec4b.im, EPS)); + assert(isClose(rec4a.re, rec4b.re, 1e-14)); + assert(isClose(rec4a.im, rec4b.im, 1e-14)); auto rer = a ^^ complex(2.0, 0.0); auto rcheck = a ^^ 2.0; @@ -558,13 +601,13 @@ if (isFloatingPoint!T) rcheck = (-a) ^^ 2.0; assert(feqrel(rer2.re, rcheck) == double.mant_dig); assert(isIdentical(rer2.re, rcheck)); - assert(approxEqual(rer2.im, 0.0, EPS)); + assert(isClose(rer2.im, 0.0, 0.0, 1e-10)); auto rer3 = (-a) ^^ complex(-2.0, 0.0); rcheck = (-a) ^^ (-2.0); assert(feqrel(rer3.re, rcheck) == double.mant_dig); assert(isIdentical(rer3.re, rcheck)); - assert(approxEqual(rer3.im, 0.0, EPS)); + assert(isClose(rer3.im, 0.0, 0.0, EPS)); auto rer4 = a ^^ complex(-2.0, 0.0); rcheck = a ^^ (-2.0); @@ -576,10 +619,10 @@ if (isFloatingPoint!T) foreach (i; 0 .. 6) { auto cei = c1^^i; - assert(approxEqual(abs(cei), abs(c1)^^i, EPS)); + assert(isClose(abs(cei), abs(c1)^^i, 1e-14)); // Use cos() here to deal with arguments that go outside // the (-pi,pi] interval (only an issue for i>3). - assert(approxEqual(std.math.cos(arg(cei)), std.math.cos(arg(c1)*i), EPS)); + assert(isClose(core.math.cos(arg(cei)), core.math.cos(arg(c1)*i), 1e-14)); } // Check operations between different complex types. @@ -596,22 +639,22 @@ if (isFloatingPoint!T) auto c2c = c2; c1c /= c1; - assert(approxEqual(c1c.re, 1.0, EPS)); - assert(approxEqual(c1c.im, 0.0, EPS)); + assert(isClose(c1c.re, 1.0, EPS)); + assert(isClose(c1c.im, 0.0, 0.0, EPS)); c1c = c1; c1c /= c2; - assert(approxEqual(c1c.re, 0.588235, EPS)); - assert(approxEqual(c1c.im, -0.352941, EPS)); + assert(isClose(c1c.re, 0.5882352941177, 1e-12)); + assert(isClose(c1c.im, -0.3529411764706, 1e-12)); c2c /= c1; - assert(approxEqual(c2c.re, 1.25, EPS)); - assert(approxEqual(c2c.im, 0.75, EPS)); + assert(isClose(c2c.re, 1.25, EPS)); + assert(isClose(c2c.im, 0.75, EPS)); c2c = c2; c2c /= c2; - assert(approxEqual(c2c.re, 1.0, EPS)); - assert(approxEqual(c2c.im, 0.0, EPS)); + assert(isClose(c2c.re, 1.0, EPS)); + assert(isClose(c2c.im, 0.0, 0.0, EPS)); } @safe pure nothrow unittest @@ -697,19 +740,39 @@ if (is(T R == Complex!R)) */ T abs(T)(Complex!T z) @safe pure nothrow @nogc { - import std.math : hypot; + import std.math.algebraic : hypot; return hypot(z.re, z.im); } /// @safe pure nothrow unittest { - static import std.math; + static import core.math; assert(abs(complex(1.0)) == 1.0); assert(abs(complex(0.0, 1.0)) == 1.0); - assert(abs(complex(1.0L, -2.0L)) == std.math.sqrt(5.0L)); + assert(abs(complex(1.0L, -2.0L)) == core.math.sqrt(5.0L)); +} + +@safe pure nothrow @nogc unittest +{ + static import core.math; + assert(abs(complex(0.0L, -3.2L)) == 3.2L); + assert(abs(complex(0.0L, 71.6L)) == 71.6L); + assert(abs(complex(-1.0L, 1.0L)) == core.math.sqrt(2.0L)); } +@safe pure nothrow @nogc unittest +{ + import std.meta : AliasSeq; + static foreach (T; AliasSeq!(float, double, real)) + {{ + static import std.math; + Complex!T a = complex(T(-12), T(3)); + T b = std.math.hypot(a.re, a.im); + assert(std.math.isClose(abs(a), b)); + assert(std.math.isClose(abs(-a), b)); + }} +} /++ Params: @@ -726,17 +789,17 @@ T sqAbs(T)(Complex!T z) @safe pure nothrow @nogc /// @safe pure nothrow unittest { - import std.math; + import std.math.operations : isClose; assert(sqAbs(complex(0.0)) == 0.0); assert(sqAbs(complex(1.0)) == 1.0); assert(sqAbs(complex(0.0, 1.0)) == 1.0); - assert(approxEqual(sqAbs(complex(1.0L, -2.0L)), 5.0L)); - assert(approxEqual(sqAbs(complex(-3.0L, 1.0L)), 10.0L)); - assert(approxEqual(sqAbs(complex(1.0f,-1.0f)), 2.0f)); + assert(isClose(sqAbs(complex(1.0L, -2.0L)), 5.0L)); + assert(isClose(sqAbs(complex(-3.0L, 1.0L)), 10.0L)); + assert(isClose(sqAbs(complex(1.0f,-1.0f)), 2.0f)); } /// ditto -T sqAbs(T)(T x) @safe pure nothrow @nogc +T sqAbs(T)(const T x) @safe pure nothrow @nogc if (isFloatingPoint!T) { return x*x; @@ -744,11 +807,11 @@ if (isFloatingPoint!T) @safe pure nothrow unittest { - import std.math; + import std.math.operations : isClose; assert(sqAbs(0.0) == 0.0); assert(sqAbs(-1.0) == 1.0); - assert(approxEqual(sqAbs(-3.0L), 9.0L)); - assert(approxEqual(sqAbs(-5.0f), 25.0f)); + assert(isClose(sqAbs(-3.0L), 9.0L)); + assert(isClose(sqAbs(-5.0f), 25.0f)); } @@ -758,14 +821,14 @@ if (isFloatingPoint!T) */ T arg(T)(Complex!T z) @safe pure nothrow @nogc { - import std.math : atan2; + import std.math.trigonometry : atan2; return atan2(z.im, z.re); } /// @safe pure nothrow unittest { - import std.math; + import std.math.constants : PI_2, PI_4; assert(arg(complex(1.0)) == 0.0); assert(arg(complex(0.0L, 1.0L)) == PI_2); assert(arg(complex(1.0L, 1.0L)) == PI_4); @@ -773,6 +836,30 @@ T arg(T)(Complex!T z) @safe pure nothrow @nogc /** + * Extracts the norm of a complex number. + * Params: + * z = A complex number + * Returns: + * The squared magnitude of `z`. + */ +T norm(T)(Complex!T z) @safe pure nothrow @nogc +{ + return z.re * z.re + z.im * z.im; +} + +/// +@safe pure nothrow @nogc unittest +{ + import std.math.operations : isClose; + import std.math.constants : PI; + assert(norm(complex(3.0, 4.0)) == 25.0); + assert(norm(fromPolar(5.0, 0.0)) == 25.0); + assert(isClose(norm(fromPolar(5.0L, PI / 6)), 25.0L)); + assert(isClose(norm(fromPolar(5.0L, 13 * PI / 6)), 25.0L)); +} + + +/** Params: z = A complex number. Returns: The complex conjugate of `z`. */ @@ -788,6 +875,43 @@ Complex!T conj(T)(Complex!T z) @safe pure nothrow @nogc assert(conj(complex(1.0, 2.0)) == complex(1.0, -2.0)); } +@safe pure nothrow @nogc unittest +{ + import std.meta : AliasSeq; + static foreach (T; AliasSeq!(float, double, real)) + {{ + auto c = Complex!T(7, 3L); + assert(conj(c) == Complex!T(7, -3L)); + auto z = Complex!T(0, -3.2L); + assert(conj(z) == -z); + }} +} + +/** + * Returns the projection of `z` onto the Riemann sphere. + * Params: + * z = A complex number + * Returns: + * The projection of `z` onto the Riemann sphere. + */ +Complex!T proj(T)(Complex!T z) +{ + static import std.math; + + if (std.math.isInfinity(z.re) || std.math.isInfinity(z.im)) + return Complex!T(T.infinity, std.math.copysign(0.0, z.im)); + + return z; +} + +/// +@safe pure nothrow unittest +{ + assert(proj(complex(1.0)) == complex(1.0)); + assert(proj(complex(double.infinity, 5.0)) == complex(double.infinity, 0.0)); + assert(proj(complex(5.0, -double.infinity)) == complex(double.infinity, -0.0)); +} + /** Constructs a complex number given its absolute value and argument. @@ -796,10 +920,10 @@ Complex!T conj(T)(Complex!T z) @safe pure nothrow @nogc argument = The argument Returns: The complex number with the given modulus and argument. */ -Complex!(CommonType!(T, U)) fromPolar(T, U)(T modulus, U argument) +Complex!(CommonType!(T, U)) fromPolar(T, U)(const T modulus, const U argument) @safe pure nothrow @nogc { - import std.math : sin, cos; + import core.math : sin, cos; return Complex!(CommonType!(T,U)) (modulus*cos(argument), modulus*sin(argument)); } @@ -807,22 +931,35 @@ Complex!(CommonType!(T, U)) fromPolar(T, U)(T modulus, U argument) /// @safe pure nothrow unittest { - import std.math; - auto z = fromPolar(std.math.sqrt(2.0), PI_4); - assert(approxEqual(z.re, 1.0L, real.epsilon)); - assert(approxEqual(z.im, 1.0L, real.epsilon)); + import core.math; + import std.math.operations : isClose; + import std.math.algebraic : sqrt; + import std.math.constants : PI_4; + auto z = fromPolar(core.math.sqrt(2.0), PI_4); + assert(isClose(z.re, 1.0L)); + assert(isClose(z.im, 1.0L)); } +version (StdUnittest) +{ + // Helper function for comparing two Complex numbers. + int ceqrel(T)(const Complex!T x, const Complex!T y) @safe pure nothrow @nogc + { + import std.math.operations : feqrel; + const r = feqrel(x.re, y.re); + const i = feqrel(x.im, y.im); + return r < i ? r : i; + } +} /** Trigonometric functions on complex numbers. Params: z = A complex number. - Returns: The sine and cosine of `z`, respectively. + Returns: The sine, cosine and tangent of `z`, respectively. */ Complex!T sin(T)(Complex!T z) @safe pure nothrow @nogc { - import std.math : expi, coshisinh; auto cs = expi(z.re); auto csh = coshisinh(z.im); return typeof(return)(cs.im * csh.re, cs.re * csh.im); @@ -831,21 +968,20 @@ Complex!T sin(T)(Complex!T z) @safe pure nothrow @nogc /// @safe pure nothrow unittest { - static import std.math; - import std.math : feqrel; + static import core.math; assert(sin(complex(0.0)) == 0.0); - assert(sin(complex(2.0, 0)) == std.math.sin(2.0)); - auto c1 = sin(complex(2.0L, 0)); - auto c2 = complex(std.math.sin(2.0L), 0); - assert(feqrel(c1.re, c2.re) >= real.mant_dig - 1 && - feqrel(c1.im, c2.im) >= real.mant_dig - 1); + assert(sin(complex(2.0, 0)) == core.math.sin(2.0)); } +@safe pure nothrow unittest +{ + static import core.math; + assert(ceqrel(sin(complex(2.0L, 0)), complex(core.math.sin(2.0L))) >= real.mant_dig - 1); +} /// ditto Complex!T cos(T)(Complex!T z) @safe pure nothrow @nogc { - import std.math : expi, coshisinh; auto cs = expi(z.re); auto csh = coshisinh(z.im); return typeof(return)(cs.re * csh.re, - cs.im * csh.im); @@ -854,18 +990,235 @@ Complex!T cos(T)(Complex!T z) @safe pure nothrow @nogc /// @safe pure nothrow unittest { + static import core.math; static import std.math; - import std.math : feqrel; assert(cos(complex(0.0)) == 1.0); - assert(cos(complex(1.3)) == std.math.cos(1.3)); - auto c1 = cos(complex(0, 5.2L)); - auto c2 = complex(std.math.cosh(5.2L), 0.0L); - assert(feqrel(c1.re, c2.re) >= real.mant_dig - 1 && - feqrel(c1.im, c2.im) >= real.mant_dig - 1); - auto c3 = cos(complex(1.3L)); - auto c4 = complex(std.math.cos(1.3L), 0.0L); - assert(feqrel(c3.re, c4.re) >= real.mant_dig - 1 && - feqrel(c3.im, c4.im) >= real.mant_dig - 1); + assert(cos(complex(1.3, 0.0)) == core.math.cos(1.3)); + assert(cos(complex(0.0, 5.2)) == std.math.cosh(5.2)); +} + +@safe pure nothrow unittest +{ + static import core.math; + static import std.math; + assert(ceqrel(cos(complex(0, 5.2L)), complex(std.math.cosh(5.2L), 0.0L)) >= real.mant_dig - 1); + assert(ceqrel(cos(complex(1.3L)), complex(core.math.cos(1.3L))) >= real.mant_dig - 1); +} + +/// ditto +Complex!T tan(T)(Complex!T z) @safe pure nothrow @nogc +{ + return sin(z) / cos(z); +} + +/// +@safe pure nothrow @nogc unittest +{ + static import std.math; + assert(ceqrel(tan(complex(1.0, 0.0)), complex(std.math.tan(1.0), 0.0)) >= double.mant_dig - 2); + assert(ceqrel(tan(complex(0.0, 1.0)), complex(0.0, std.math.tanh(1.0))) >= double.mant_dig - 2); +} + +/** + Inverse trigonometric functions on complex numbers. + + Params: z = A complex number. + Returns: The arcsine, arccosine and arctangent of `z`, respectively. +*/ +Complex!T asin(T)(Complex!T z) @safe pure nothrow @nogc +{ + auto ash = asinh(Complex!T(-z.im, z.re)); + return Complex!T(ash.im, -ash.re); +} + +/// +@safe pure nothrow unittest +{ + import std.math.operations : isClose; + import std.math.constants : PI; + assert(asin(complex(0.0)) == 0.0); + assert(isClose(asin(complex(0.5L)), PI / 6)); +} + +@safe pure nothrow unittest +{ + import std.math.operations : isClose; + import std.math.constants : PI; + version (DigitalMars) {} else // Disabled because of issue 21376 + assert(isClose(asin(complex(0.5f)), float(PI) / 6)); +} + +/// ditto +Complex!T acos(T)(Complex!T z) @safe pure nothrow @nogc +{ + static import std.math; + auto as = asin(z); + return Complex!T(T(std.math.PI_2) - as.re, as.im); +} + +/// +@safe pure nothrow unittest +{ + import std.math.operations : isClose; + import std.math.constants : PI; + import std.math.trigonometry : std_math_acos = acos; + assert(acos(complex(0.0)) == std_math_acos(0.0)); + assert(isClose(acos(complex(0.5L)), PI / 3)); +} + +@safe pure nothrow unittest +{ + import std.math.operations : isClose; + import std.math.constants : PI; + version (DigitalMars) {} else // Disabled because of issue 21376 + assert(isClose(acos(complex(0.5f)), float(PI) / 3)); +} + +/// ditto +Complex!T atan(T)(Complex!T z) @safe pure nothrow @nogc +{ + static import std.math; + const T re2 = z.re * z.re; + const T x = 1 - re2 - z.im * z.im; + + T num = z.im + 1; + T den = z.im - 1; + + num = re2 + num * num; + den = re2 + den * den; + + return Complex!T(T(0.5) * std.math.atan2(2 * z.re, x), + T(0.25) * std.math.log(num / den)); +} + +/// +@safe pure nothrow @nogc unittest +{ + import std.math.operations : isClose; + import std.math.constants : PI; + assert(atan(complex(0.0)) == 0.0); + assert(isClose(atan(sqrt(complex(3.0L))), PI / 3)); + assert(isClose(atan(sqrt(complex(3.0f))), float(PI) / 3)); +} + +/** + Hyperbolic trigonometric functions on complex numbers. + + Params: z = A complex number. + Returns: The hyperbolic sine, cosine and tangent of `z`, respectively. +*/ +Complex!T sinh(T)(Complex!T z) @safe pure nothrow @nogc +{ + static import core.math, std.math; + return Complex!T(std.math.sinh(z.re) * core.math.cos(z.im), + std.math.cosh(z.re) * core.math.sin(z.im)); +} + +/// +@safe pure nothrow unittest +{ + static import std.math; + assert(sinh(complex(0.0)) == 0.0); + assert(sinh(complex(1.0L)) == std.math.sinh(1.0L)); + assert(sinh(complex(1.0f)) == std.math.sinh(1.0f)); +} + +/// ditto +Complex!T cosh(T)(Complex!T z) @safe pure nothrow @nogc +{ + static import core.math, std.math; + return Complex!T(std.math.cosh(z.re) * core.math.cos(z.im), + std.math.sinh(z.re) * core.math.sin(z.im)); +} + +/// +@safe pure nothrow unittest +{ + static import std.math; + assert(cosh(complex(0.0)) == 1.0); + assert(cosh(complex(1.0L)) == std.math.cosh(1.0L)); + assert(cosh(complex(1.0f)) == std.math.cosh(1.0f)); +} + +/// ditto +Complex!T tanh(T)(Complex!T z) @safe pure nothrow @nogc +{ + return sinh(z) / cosh(z); +} + +/// +@safe pure nothrow @nogc unittest +{ + import std.math.operations : isClose; + import std.math.trigonometry : std_math_tanh = tanh; + assert(tanh(complex(0.0)) == 0.0); + assert(isClose(tanh(complex(1.0L)), std_math_tanh(1.0L))); + assert(isClose(tanh(complex(1.0f)), std_math_tanh(1.0f))); +} + +/** + Inverse hyperbolic trigonometric functions on complex numbers. + + Params: z = A complex number. + Returns: The hyperbolic arcsine, arccosine and arctangent of `z`, respectively. +*/ +Complex!T asinh(T)(Complex!T z) @safe pure nothrow @nogc +{ + auto t = Complex!T((z.re - z.im) * (z.re + z.im) + 1, 2 * z.re * z.im); + return log(sqrt(t) + z); +} + +/// +@safe pure nothrow unittest +{ + import std.math.operations : isClose; + import std.math.trigonometry : std_math_asinh = asinh; + assert(asinh(complex(0.0)) == 0.0); + assert(isClose(asinh(complex(1.0L)), std_math_asinh(1.0L))); + assert(isClose(asinh(complex(1.0f)), std_math_asinh(1.0f))); +} + +/// ditto +Complex!T acosh(T)(Complex!T z) @safe pure nothrow @nogc +{ + return 2 * log(sqrt(T(0.5) * (z + 1)) + sqrt(T(0.5) * (z - 1))); +} + +/// +@safe pure nothrow unittest +{ + import std.math.operations : isClose; + import std.math.trigonometry : std_math_acosh = acosh; + assert(acosh(complex(1.0)) == 0.0); + assert(isClose(acosh(complex(3.0L)), std_math_acosh(3.0L))); + assert(isClose(acosh(complex(3.0f)), std_math_acosh(3.0f))); +} + +/// ditto +Complex!T atanh(T)(Complex!T z) @safe pure nothrow @nogc +{ + static import std.math; + const T im2 = z.im * z.im; + const T x = 1 - im2 - z.re * z.re; + + T num = 1 + z.re; + T den = 1 - z.re; + + num = im2 + num * num; + den = im2 + den * den; + + return Complex!T(T(0.25) * (std.math.log(num) - std.math.log(den)), + T(0.5) * std.math.atan2(2 * z.im, x)); +} + +/// +@safe pure nothrow @nogc unittest +{ + import std.math.operations : isClose; + import std.math.trigonometry : std_math_atanh = atanh; + assert(atanh(complex(0.0)) == 0.0); + assert(isClose(atanh(complex(0.5L)), std_math_atanh(0.5L))); + assert(isClose(atanh(complex(0.5f)), std_math_atanh(0.5f))); } /** @@ -873,29 +1226,50 @@ Complex!T cos(T)(Complex!T z) @safe pure nothrow @nogc Returns: The value of cos(y) + i sin(y). Note: - $(D expi) is included here for convenience and for easy migration of code - that uses $(REF _expi, std,math). Unlike $(REF _expi, std,math), which uses the - x87 $(I fsincos) instruction when possible, this function is no faster - than calculating cos(y) and sin(y) separately. + `expi` is included here for convenience and for easy migration of code. */ Complex!real expi(real y) @trusted pure nothrow @nogc { - import std.math : cos, sin; + import core.math : cos, sin; return Complex!real(cos(y), sin(y)); } /// @safe pure nothrow unittest { - static import std.math; - - assert(expi(1.3e5L) == complex(std.math.cos(1.3e5L), std.math.sin(1.3e5L))); + import core.math : cos, sin; assert(expi(0.0L) == 1.0L); - auto z1 = expi(1.234); - auto z2 = std.math.expi(1.234); - assert(z1.re == z2.re && z1.im == z2.im); + assert(expi(1.3e5L) == complex(cos(1.3e5L), sin(1.3e5L))); } +/** + Params: y = A real number. + Returns: The value of cosh(y) + i sinh(y) + + Note: + `coshisinh` is included here for convenience and for easy migration of code. +*/ +Complex!real coshisinh(real y) @safe pure nothrow @nogc +{ + static import core.math; + static import std.math; + if (core.math.fabs(y) <= 0.5) + return Complex!real(std.math.cosh(y), std.math.sinh(y)); + else + { + auto z = std.math.exp(y); + auto zi = 0.5 / z; + z = 0.5 * z; + return Complex!real(z + zi, z - zi); + } +} + +/// +@safe pure nothrow @nogc unittest +{ + import std.math.trigonometry : cosh, sinh; + assert(coshisinh(3.0L) == complex(cosh(3.0L), sinh(3.0L))); +} /** Params: z = A complex number. @@ -903,7 +1277,7 @@ Complex!real expi(real y) @trusted pure nothrow @nogc */ Complex!T sqrt(T)(Complex!T z) @safe pure nothrow @nogc { - static import std.math; + static import core.math; typeof(return) c; real x,y,w,r; @@ -916,19 +1290,19 @@ Complex!T sqrt(T)(Complex!T z) @safe pure nothrow @nogc real z_re = z.re; real z_im = z.im; - x = std.math.fabs(z_re); - y = std.math.fabs(z_im); + x = core.math.fabs(z_re); + y = core.math.fabs(z_im); if (x >= y) { r = y / x; - w = std.math.sqrt(x) - * std.math.sqrt(0.5 * (1 + std.math.sqrt(1 + r * r))); + w = core.math.sqrt(x) + * core.math.sqrt(0.5 * (1 + core.math.sqrt(1 + r * r))); } else { r = x / y; - w = std.math.sqrt(y) - * std.math.sqrt(0.5 * (r + std.math.sqrt(1 + r * r))); + w = core.math.sqrt(y) + * core.math.sqrt(0.5 * (r + core.math.sqrt(1 + r * r))); } if (z_re >= 0) @@ -948,29 +1322,31 @@ Complex!T sqrt(T)(Complex!T z) @safe pure nothrow @nogc /// @safe pure nothrow unittest { - static import std.math; + static import core.math; assert(sqrt(complex(0.0)) == 0.0); - assert(sqrt(complex(1.0L, 0)) == std.math.sqrt(1.0L)); + assert(sqrt(complex(1.0L, 0)) == core.math.sqrt(1.0L)); assert(sqrt(complex(-1.0L, 0)) == complex(0, 1.0L)); + assert(sqrt(complex(-8.0, -6.0)) == complex(1.0, -3.0)); } @safe pure nothrow unittest { - import std.math : approxEqual; + import std.math.operations : isClose; auto c1 = complex(1.0, 1.0); auto c2 = Complex!double(0.5, 2.0); auto c1s = sqrt(c1); - assert(approxEqual(c1s.re, 1.09868411)); - assert(approxEqual(c1s.im, 0.45508986)); + assert(isClose(c1s.re, 1.09868411347)); + assert(isClose(c1s.im, 0.455089860562)); auto c2s = sqrt(c2); - assert(approxEqual(c2s.re, 1.1317134)); - assert(approxEqual(c2s.im, 0.8836155)); + assert(isClose(c2s.re, 1.13171392428)); + assert(isClose(c2s.im, 0.883615530876)); } -// Issue 10881: support %f formatting of complex numbers +// support %f formatting of complex numbers +// https://issues.dlang.org/show_bug.cgi?id=10881 @safe unittest { import std.format : format; @@ -985,7 +1361,7 @@ Complex!T sqrt(T)(Complex!T z) @safe pure nothrow @nogc @safe unittest { // Test wide string formatting - import std.format; + import std.format.write : formattedWrite; wstring wformat(T)(string format, Complex!T c) { import std.array : appender; @@ -1003,3 +1379,526 @@ Complex!T sqrt(T)(Complex!T z) @safe pure nothrow @nogc // Test ease of use (vanilla toString() should be supported) assert(complex(1.2, 3.4).toString() == "1.2+3.4i"); } + +@safe pure nothrow @nogc unittest +{ + auto c = complex(3.0L, 4.0L); + c = sqrt(c); + assert(c.re == 2.0L); + assert(c.im == 1.0L); +} + +/** + * Calculates e$(SUPERSCRIPT x). + * Params: + * x = A complex number + * Returns: + * The complex base e exponential of `x` + * + * $(TABLE_SV + * $(TR $(TH x) $(TH exp(x))) + * $(TR $(TD ($(PLUSMN)0, +0)) $(TD (1, +0))) + * $(TR $(TD (any, +$(INFIN))) $(TD ($(NAN), $(NAN)))) + * $(TR $(TD (any, $(NAN)) $(TD ($(NAN), $(NAN))))) + * $(TR $(TD (+$(INFIN), +0)) $(TD (+$(INFIN), +0))) + * $(TR $(TD (-$(INFIN), any)) $(TD ($(PLUSMN)0, cis(x.im)))) + * $(TR $(TD (+$(INFIN), any)) $(TD ($(PLUSMN)$(INFIN), cis(x.im)))) + * $(TR $(TD (-$(INFIN), +$(INFIN))) $(TD ($(PLUSMN)0, $(PLUSMN)0))) + * $(TR $(TD (+$(INFIN), +$(INFIN))) $(TD ($(PLUSMN)$(INFIN), $(NAN)))) + * $(TR $(TD (-$(INFIN), $(NAN))) $(TD ($(PLUSMN)0, $(PLUSMN)0))) + * $(TR $(TD (+$(INFIN), $(NAN))) $(TD ($(PLUSMN)$(INFIN), $(NAN)))) + * $(TR $(TD ($(NAN), +0)) $(TD ($(NAN), +0))) + * $(TR $(TD ($(NAN), any)) $(TD ($(NAN), $(NAN)))) + * $(TR $(TD ($(NAN), $(NAN))) $(TD ($(NAN), $(NAN)))) + * ) + */ +Complex!T exp(T)(Complex!T x) @trusted pure nothrow @nogc // TODO: @safe +{ + static import std.math; + + // Handle special cases explicitly here, as fromPolar will otherwise + // cause them to return Complex!T(NaN, NaN), or with the wrong sign. + if (std.math.isInfinity(x.re)) + { + if (std.math.isNaN(x.im)) + { + if (std.math.signbit(x.re)) + return Complex!T(0, std.math.copysign(0, x.im)); + else + return x; + } + if (std.math.isInfinity(x.im)) + { + if (std.math.signbit(x.re)) + return Complex!T(0, std.math.copysign(0, x.im)); + else + return Complex!T(T.infinity, -T.nan); + } + if (x.im == 0.0) + { + if (std.math.signbit(x.re)) + return Complex!T(0.0); + else + return Complex!T(T.infinity); + } + } + if (std.math.isNaN(x.re)) + { + if (std.math.isNaN(x.im) || std.math.isInfinity(x.im)) + return Complex!T(T.nan, T.nan); + if (x.im == 0.0) + return x; + } + if (x.re == 0.0) + { + if (std.math.isNaN(x.im) || std.math.isInfinity(x.im)) + return Complex!T(T.nan, T.nan); + if (x.im == 0.0) + return Complex!T(1.0, 0.0); + } + + return fromPolar!(T, T)(std.math.exp(x.re), x.im); +} + +/// +@safe pure nothrow @nogc unittest +{ + import std.math.operations : isClose; + import std.math.constants : PI; + + assert(exp(complex(0.0, 0.0)) == complex(1.0, 0.0)); + + auto a = complex(2.0, 1.0); + assert(exp(conj(a)) == conj(exp(a))); + + auto b = exp(complex(0.0L, 1.0L) * PI); + assert(isClose(b, -1.0L, 0.0, 1e-15)); +} + +@safe pure nothrow @nogc unittest +{ + import std.math.traits : isNaN, isInfinity; + + auto a = exp(complex(0.0, double.infinity)); + assert(a.re.isNaN && a.im.isNaN); + auto b = exp(complex(0.0, double.infinity)); + assert(b.re.isNaN && b.im.isNaN); + auto c = exp(complex(0.0, double.nan)); + assert(c.re.isNaN && c.im.isNaN); + + auto d = exp(complex(+double.infinity, 0.0)); + assert(d == complex(double.infinity, 0.0)); + auto e = exp(complex(-double.infinity, 0.0)); + assert(e == complex(0.0)); + auto f = exp(complex(-double.infinity, 1.0)); + assert(f == complex(0.0)); + auto g = exp(complex(+double.infinity, 1.0)); + assert(g == complex(double.infinity, double.infinity)); + auto h = exp(complex(-double.infinity, +double.infinity)); + assert(h == complex(0.0)); + auto i = exp(complex(+double.infinity, +double.infinity)); + assert(i.re.isInfinity && i.im.isNaN); + auto j = exp(complex(-double.infinity, double.nan)); + assert(j == complex(0.0)); + auto k = exp(complex(+double.infinity, double.nan)); + assert(k.re.isInfinity && k.im.isNaN); + + auto l = exp(complex(double.nan, 0)); + assert(l.re.isNaN && l.im == 0.0); + auto m = exp(complex(double.nan, 1)); + assert(m.re.isNaN && m.im.isNaN); + auto n = exp(complex(double.nan, double.nan)); + assert(n.re.isNaN && n.im.isNaN); +} + +@safe pure nothrow @nogc unittest +{ + import std.math.constants : PI; + import std.math.operations : isClose; + + auto a = exp(complex(0.0, -PI)); + assert(isClose(a, -1.0, 0.0, 1e-15)); + + auto b = exp(complex(0.0, -2.0 * PI / 3.0)); + assert(isClose(b, complex(-0.5L, -0.866025403784438646763L))); + + auto c = exp(complex(0.0, PI / 3.0)); + assert(isClose(c, complex(0.5L, 0.866025403784438646763L))); + + auto d = exp(complex(0.0, 2.0 * PI / 3.0)); + assert(isClose(d, complex(-0.5L, 0.866025403784438646763L))); + + auto e = exp(complex(0.0, PI)); + assert(isClose(e, -1.0, 0.0, 1e-15)); +} + +/** + * Calculate the natural logarithm of x. + * The branch cut is along the negative axis. + * Params: + * x = A complex number + * Returns: + * The complex natural logarithm of `x` + * + * $(TABLE_SV + * $(TR $(TH x) $(TH log(x))) + * $(TR $(TD (-0, +0)) $(TD (-$(INFIN), $(PI)))) + * $(TR $(TD (+0, +0)) $(TD (-$(INFIN), +0))) + * $(TR $(TD (any, +$(INFIN))) $(TD (+$(INFIN), $(PI)/2))) + * $(TR $(TD (any, $(NAN))) $(TD ($(NAN), $(NAN)))) + * $(TR $(TD (-$(INFIN), any)) $(TD (+$(INFIN), $(PI)))) + * $(TR $(TD (+$(INFIN), any)) $(TD (+$(INFIN), +0))) + * $(TR $(TD (-$(INFIN), +$(INFIN))) $(TD (+$(INFIN), 3$(PI)/4))) + * $(TR $(TD (+$(INFIN), +$(INFIN))) $(TD (+$(INFIN), $(PI)/4))) + * $(TR $(TD ($(PLUSMN)$(INFIN), $(NAN))) $(TD (+$(INFIN), $(NAN)))) + * $(TR $(TD ($(NAN), any)) $(TD ($(NAN), $(NAN)))) + * $(TR $(TD ($(NAN), +$(INFIN))) $(TD (+$(INFIN), $(NAN)))) + * $(TR $(TD ($(NAN), $(NAN))) $(TD ($(NAN), $(NAN)))) + * ) + */ +Complex!T log(T)(Complex!T x) @safe pure nothrow @nogc +{ + static import std.math; + + // Handle special cases explicitly here for better accuracy. + // The order here is important, so that the correct path is chosen. + if (std.math.isNaN(x.re)) + { + if (std.math.isInfinity(x.im)) + return Complex!T(T.infinity, T.nan); + else + return Complex!T(T.nan, T.nan); + } + if (std.math.isInfinity(x.re)) + { + if (std.math.isNaN(x.im)) + return Complex!T(T.infinity, T.nan); + else if (std.math.isInfinity(x.im)) + { + if (std.math.signbit(x.re)) + return Complex!T(T.infinity, std.math.copysign(3.0 * std.math.PI_4, x.im)); + else + return Complex!T(T.infinity, std.math.copysign(std.math.PI_4, x.im)); + } + else + { + if (std.math.signbit(x.re)) + return Complex!T(T.infinity, std.math.copysign(std.math.PI, x.im)); + else + return Complex!T(T.infinity, std.math.copysign(0.0, x.im)); + } + } + if (std.math.isNaN(x.im)) + return Complex!T(T.nan, T.nan); + if (std.math.isInfinity(x.im)) + return Complex!T(T.infinity, std.math.copysign(std.math.PI_2, x.im)); + if (x.re == 0.0 && x.im == 0.0) + { + if (std.math.signbit(x.re)) + return Complex!T(-T.infinity, std.math.copysign(std.math.PI, x.im)); + else + return Complex!T(-T.infinity, std.math.copysign(0.0, x.im)); + } + + return Complex!T(std.math.log(abs(x)), arg(x)); +} + +/// +@safe pure nothrow @nogc unittest +{ + import core.math : sqrt; + import std.math.constants : PI; + import std.math.operations : isClose; + + auto a = complex(2.0, 1.0); + assert(log(conj(a)) == conj(log(a))); + + auto b = 2.0 * log10(complex(0.0, 1.0)); + auto c = 4.0 * log10(complex(sqrt(2.0) / 2, sqrt(2.0) / 2)); + assert(isClose(b, c, 0.0, 1e-15)); + + assert(log(complex(-1.0L, 0.0L)) == complex(0.0L, PI)); + assert(log(complex(-1.0L, -0.0L)) == complex(0.0L, -PI)); +} + +@safe pure nothrow @nogc unittest +{ + import std.math.traits : isNaN, isInfinity; + import std.math.constants : PI, PI_2, PI_4; + + auto a = log(complex(-0.0L, 0.0L)); + assert(a == complex(-real.infinity, PI)); + auto b = log(complex(0.0L, 0.0L)); + assert(b == complex(-real.infinity, +0.0L)); + auto c = log(complex(1.0L, real.infinity)); + assert(c == complex(real.infinity, PI_2)); + auto d = log(complex(1.0L, real.nan)); + assert(d.re.isNaN && d.im.isNaN); + + auto e = log(complex(-real.infinity, 1.0L)); + assert(e == complex(real.infinity, PI)); + auto f = log(complex(real.infinity, 1.0L)); + assert(f == complex(real.infinity, 0.0L)); + auto g = log(complex(-real.infinity, real.infinity)); + assert(g == complex(real.infinity, 3.0 * PI_4)); + auto h = log(complex(real.infinity, real.infinity)); + assert(h == complex(real.infinity, PI_4)); + auto i = log(complex(real.infinity, real.nan)); + assert(i.re.isInfinity && i.im.isNaN); + + auto j = log(complex(real.nan, 1.0L)); + assert(j.re.isNaN && j.im.isNaN); + auto k = log(complex(real.nan, real.infinity)); + assert(k.re.isInfinity && k.im.isNaN); + auto l = log(complex(real.nan, real.nan)); + assert(l.re.isNaN && l.im.isNaN); +} + +@safe pure nothrow @nogc unittest +{ + import std.math.constants : PI; + import std.math.operations : isClose; + + auto a = log(fromPolar(1.0, PI / 6.0)); + assert(isClose(a, complex(0.0L, 0.523598775598298873077L), 0.0, 1e-15)); + + auto b = log(fromPolar(1.0, PI / 3.0)); + assert(isClose(b, complex(0.0L, 1.04719755119659774615L), 0.0, 1e-15)); + + auto c = log(fromPolar(1.0, PI / 2.0)); + assert(isClose(c, complex(0.0L, 1.57079632679489661923L), 0.0, 1e-15)); + + auto d = log(fromPolar(1.0, 2.0 * PI / 3.0)); + assert(isClose(d, complex(0.0L, 2.09439510239319549230L), 0.0, 1e-15)); + + auto e = log(fromPolar(1.0, 5.0 * PI / 6.0)); + assert(isClose(e, complex(0.0L, 2.61799387799149436538L), 0.0, 1e-15)); + + auto f = log(complex(-1.0L, 0.0L)); + assert(isClose(f, complex(0.0L, PI), 0.0, 1e-15)); +} + +/** + * Calculate the base-10 logarithm of x. + * Params: + * x = A complex number + * Returns: + * The complex base 10 logarithm of `x` + */ +Complex!T log10(T)(Complex!T x) @safe pure nothrow @nogc +{ + static import std.math; + + return log(x) / Complex!T(std.math.log(10.0)); +} + +/// +@safe pure nothrow @nogc unittest +{ + import core.math : sqrt; + import std.math.constants : LN10, PI; + import std.math.operations : isClose; + + auto a = complex(2.0, 1.0); + assert(log10(a) == log(a) / log(complex(10.0))); + + auto b = log10(complex(0.0, 1.0)) * 2.0; + auto c = log10(complex(sqrt(2.0) / 2, sqrt(2.0) / 2)) * 4.0; + assert(isClose(b, c, 0.0, 1e-15)); + + assert(ceqrel(log10(complex(-100.0L, 0.0L)), complex(2.0L, PI / LN10)) >= real.mant_dig - 1); + assert(ceqrel(log10(complex(-100.0L, -0.0L)), complex(2.0L, -PI / LN10)) >= real.mant_dig - 1); +} + +@safe pure nothrow @nogc unittest +{ + import std.math.constants : PI; + import std.math.operations : isClose; + + auto a = log10(fromPolar(1.0, PI / 6.0)); + assert(isClose(a, complex(0.0L, 0.227396058973640224580L), 0.0, 1e-15)); + + auto b = log10(fromPolar(1.0, PI / 3.0)); + assert(isClose(b, complex(0.0L, 0.454792117947280449161L), 0.0, 1e-15)); + + auto c = log10(fromPolar(1.0, PI / 2.0)); + assert(isClose(c, complex(0.0L, 0.682188176920920673742L), 0.0, 1e-15)); + + auto d = log10(fromPolar(1.0, 2.0 * PI / 3.0)); + assert(isClose(d, complex(0.0L, 0.909584235894560898323L), 0.0, 1e-15)); + + auto e = log10(fromPolar(1.0, 5.0 * PI / 6.0)); + assert(isClose(e, complex(0.0L, 1.13698029486820112290L), 0.0, 1e-15)); + + auto f = log10(complex(-1.0L, 0.0L)); + assert(isClose(f, complex(0.0L, 1.36437635384184134748L), 0.0, 1e-15)); +} + +/** + * Calculates x$(SUPERSCRIPT n). + * The branch cut is on the negative axis. + * Params: + * x = base + * n = exponent + * Returns: + * `x` raised to the power of `n` + */ +Complex!T pow(T, Int)(Complex!T x, const Int n) @safe pure nothrow @nogc +if (isIntegral!Int) +{ + alias UInt = Unsigned!(Unqual!Int); + + UInt m = (n < 0) ? -cast(UInt) n : n; + Complex!T y = (m % 2) ? x : Complex!T(1); + + while (m >>= 1) + { + x *= x; + if (m % 2) + y *= x; + } + + return (n < 0) ? Complex!T(1) / y : y; +} + +/// +@safe pure nothrow @nogc unittest +{ + import std.math.operations : isClose; + + auto a = complex(1.0, 2.0); + assert(pow(a, 2) == a * a); + assert(pow(a, 3) == a * a * a); + assert(pow(a, -2) == 1.0 / (a * a)); + assert(isClose(pow(a, -3), 1.0 / (a * a * a))); + + auto b = complex(2.0); + assert(ceqrel(pow(b, 3), exp(3 * log(b))) >= double.mant_dig - 1); +} + +/// ditto +Complex!T pow(T)(Complex!T x, const T n) @trusted pure nothrow @nogc +{ + static import std.math; + + if (x == 0.0) + return Complex!T(0.0); + + if (x.im == 0 && x.re > 0.0) + return Complex!T(std.math.pow(x.re, n)); + + Complex!T t = log(x); + return fromPolar!(T, T)(std.math.exp(n * t.re), n * t.im); +} + +/// +@safe pure nothrow @nogc unittest +{ + import std.math.operations : isClose; + assert(pow(complex(0.0), 2.0) == complex(0.0)); + assert(pow(complex(5.0), 2.0) == complex(25.0)); + + auto a = pow(complex(-1.0, 0.0), 0.5); + assert(isClose(a, complex(0.0, +1.0), 0.0, 1e-16)); + + auto b = pow(complex(-1.0, -0.0), 0.5); + assert(isClose(b, complex(0.0, -1.0), 0.0, 1e-16)); +} + +/// ditto +Complex!T pow(T)(Complex!T x, Complex!T y) @trusted pure nothrow @nogc +{ + return (x == 0) ? Complex!T(0) : exp(y * log(x)); +} + +/// +@safe pure nothrow @nogc unittest +{ + import std.math.operations : isClose; + import std.math.exponential : exp; + import std.math.constants : PI; + auto a = complex(0.0); + auto b = complex(2.0); + assert(pow(a, b) == complex(0.0)); + + auto c = complex(0.0L, 1.0L); + assert(isClose(pow(c, c), exp((-PI) / 2))); +} + +/// ditto +Complex!T pow(T)(const T x, Complex!T n) @trusted pure nothrow @nogc +{ + static import std.math; + + return (x > 0.0) + ? fromPolar!(T, T)(std.math.pow(x, n.re), n.im * std.math.log(x)) + : pow(Complex!T(x), n); +} + +/// +@safe pure nothrow @nogc unittest +{ + import std.math.operations : isClose; + assert(pow(2.0, complex(0.0)) == complex(1.0)); + assert(pow(2.0, complex(5.0)) == complex(32.0)); + + auto a = pow(-2.0, complex(-1.0)); + assert(isClose(a, complex(-0.5), 0.0, 1e-16)); + + auto b = pow(-0.5, complex(-1.0)); + assert(isClose(b, complex(-2.0), 0.0, 1e-15)); +} + +@safe pure nothrow @nogc unittest +{ + import std.math.constants : PI; + import std.math.operations : isClose; + + auto a = pow(complex(3.0, 4.0), 2); + assert(isClose(a, complex(-7.0, 24.0))); + + auto b = pow(complex(3.0, 4.0), PI); + assert(ceqrel(b, complex(-152.91512205297134, 35.547499631917738)) >= double.mant_dig - 3); + + auto c = pow(complex(3.0, 4.0), complex(-2.0, 1.0)); + assert(ceqrel(c, complex(0.015351734187477306, -0.0038407695456661503)) >= double.mant_dig - 3); + + auto d = pow(PI, complex(2.0, -1.0)); + assert(ceqrel(d, complex(4.0790296880118296, -8.9872469554541869)) >= double.mant_dig - 1); +} + +@safe pure nothrow @nogc unittest +{ + import std.meta : AliasSeq; + import std.math : RealFormat, floatTraits; + static foreach (T; AliasSeq!(float, double, real)) + {{ + static if (floatTraits!T.realFormat == RealFormat.ibmExtended) + { + /* For IBM real, epsilon is too small (since 1.0 plus any double is + representable) to be able to expect results within epsilon * 100. */ + } + else + { + T eps = T.epsilon * 100; + + T a = -1.0; + T b = 0.5; + Complex!T ref1 = pow(complex(a), complex(b)); + Complex!T res1 = pow(a, complex(b)); + Complex!T res2 = pow(complex(a), b); + assert(abs(ref1 - res1) < eps); + assert(abs(ref1 - res2) < eps); + assert(abs(res1 - res2) < eps); + + T c = -3.2; + T d = 1.4; + Complex!T ref2 = pow(complex(a), complex(b)); + Complex!T res3 = pow(a, complex(b)); + Complex!T res4 = pow(complex(a), b); + assert(abs(ref2 - res3) < eps); + assert(abs(ref2 - res4) < eps); + assert(abs(res3 - res4) < eps); + } + }} +} |