Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Missing math functions in std.complex
import std.complex;
import std.math : log, ceil, floor, round, sinh, cosh, trunc, sin, cos, acos, PI, E, copysign, isInfinity, isNaN,
// sgn needs to be renamed as std.math.sgn greedily accepts any type whatsoever and fails spectacularly when used with non-builtins.
fsgn = sgn;
import std.traits : isFloatingPoint;
// Formulas from
// http://scipp.ucsc.edu/~haber/archives/physics116A10/arc_10.pdf
// http://www.suitcaseofdreams.net/Inverse_Functions.htm
// http://www.suitcaseofdreams.net/inverse__hyperbolic_functions.htm
enum isComplex(T) = is(T == Complex!F, F);
// Useful constant:
static immutable Complex!float I = Complex!float(0, 1);
auto log(T)(Complex!T x) {
return log(abs(x)) + arg(x)*I;
}
unittest {
assert(log(0.8 + 0.19*I).approxEqual(-0.195707 + 0.23318*I));
assert(log(-2 + 0*I).approxEqual(0.69314718 + PI*I));
}
auto log10(T)(Complex!T x) {
return log(x) / log(10);
}
unittest {
assert(log10(5.7 + 2.8*I).approxEqual(0.802814 + 0.198301*I));
}
auto log2(T)(Complex!T x) {
return log(x) / log(2);
}
unittest {
assert(log2(-22.7 + 49.6*I).approxEqual(5.769446 + 2.885395*I));
}
auto exp(T)(Complex!T x) {
return E ^^ x;
}
unittest {
assert(exp(4 + 2*I).approxEqual(-22.7208484 + 49.645957*I));
}
auto tan(T)(Complex!T x) {
return(sin(2 * x.re) + sinh(2 * x.im)*I) /
(cos(2 * x.re) + cosh(2 * x.im));
}
unittest {
assert(tan(1.4 + 0.06*I).approxEqual(5.15475 + 1.85098*I));
}
auto atan(T)(Complex!T x) {
return -I/2 * log((I - x) / (I + x));
}
unittest {
assert(atan(7.19 + 6.4*I).approxEqual(1.49299 + 0.0687654*I));
}
auto asin(T)(Complex!T x) {
return -I*log(I*x + sqrt(1-x*x));
}
unittest {
assert(asin(5.2 + 7.2*I).approxEqual(0.622485 + 2.87812*I));
}
auto acos(T)(Complex!T x) {
return -I*log(x + sqrt(x*x-1));
}
unittest {
assert(acos(5.2 + 7.2*I).approxEqual(0.948311 - 2.87812*I));
assert(acos(0 - 4*I).approxEqual(1.5707963 + 2.094712547*I));
assert(acos(0.5 + 0*I).approxEqual(acos(0.5f)));
}
auto tanh(T)(Complex!T x) {
return (exp(2*x) -1) / (exp(2*x) + 1);
}
unittest {
assert(tanh(1.2 + 13*I).approxEqual(0.8811 + 0.122917*I));
}
auto atanh(T)(Complex!T x) {
return log((1 + x) / (1 - x))/2;
}
unittest {
assert(atanh(16.7 + 3.24*I).approxEqual(0.05776498 + 1.559563*I));
}
auto sinh(T)(Complex!T x) {
return (exp(x) - exp(-x)) / 2;
}
unittest {
assert(sinh(4.3 - 1.5*I).approxEqual(2.60618 - 36.7644*I));
}
auto asinh(T)(Complex!T x) {
return log(x + sqrt(x*x + 1));
}
unittest {
assert(asinh(3 + 17.2*I).approxEqual(3.552268 + 1.397837*I));
}
auto cosh(T)(Complex!T x) {
return (exp(x) + exp(-x)) / 2;
}
unittest {
assert(cosh(5.5 + 3*I).approxEqual(-121.124 + 17.2652*I));
}
auto acosh(T)(Complex!T x) {
return log(x + sqrt(x+1)*sqrt(x-1));
}
unittest {
assert(acosh(-8.4 + 2.7 * I).approxEqual(2.867924 + 2.828709*I));
}
auto cbrt(T)(Complex!T x) {
return x ^^ (1./3);
}
unittest {
assert(cbrt(4.9 + 6.3*I).approxEqual(1.90725 + 0.596781*I));
}
auto pow(T1, T2)(T1 x, T2 y)
if ((isComplex!T1 || isFloatingPoint!T1) && (isComplex!T2 || isFloatingPoint!T2) && (isComplex!T1 || isComplex!T2))
{
return exp(y * log(x));
}
unittest {
assert(pow(2.3+4.7*I, 4 + I).approxEqual(242.99 - 40.4699*I));
}
auto ceil(T)(Complex!T x) {
return ceil(x.re) + ceil(x.im)*I;
}
unittest {
assert(ceil( 1.2+1.6*I) == 2 + 2*I);
assert(ceil(-1.2-1.6*I) == -1 - I);
}
auto floor(T)(Complex!T x) {
return floor(x.re) + floor(x.im)*I;
}
unittest {
assert(floor( 1.2+1.6*I) == 1 + I);
assert(floor(-1.2-1.6*I) == -2 - 2*I);
}
auto round(T)(Complex!T x) {
return round(x.re) + round(x.im)*I;
}
unittest {
assert(round( 1.2+1.6*I) == 1 + 2*I);
assert(round(-1.2-1.6*I) == -1 - 2*I);
}
auto trunc(T)(Complex!T x) {
return trunc(x.re) + trunc(x.im)*I;
}
unittest {
assert(trunc( 1.2+1.6*I) == 1 + I);
assert(trunc(-1.2-1.6*I) == -1 - I);
}
auto atan2(T1, T2)(T1 x, T2 y)
if ((isComplex!T1 || isFloatingPoint!T1) && (isComplex!T2 || isFloatingPoint!T2) && (isComplex!T1 || isComplex!T2))
{
return -I*log((y + x*I) / sqrt(x*x + y*y));
}
unittest {
assert(atan2(1.6 + 2.4*I, .75 + .24*I).approxEqual(1.35469 + .164029*I));
}
auto sgn(T)(Complex!T x) {
return .fsgn(x.re) + .fsgn(x.im)*I;
}
unittest {
assert(sgn( 2 + 3*I) == 1 + I);
assert(sgn( 2 - 3*I) == 1 - I);
assert(sgn(-2 + 3*I) == -1 + I);
assert(sgn(-2 - 3*I) == -1 - I);
}
auto copysign(R, X)(Complex!R to, Complex!X from) {
return .copysign(to.re, from.re) + .copysign(to.im, from.im)*I;
}
unittest {
assert(copysign( 1 + I, 1 + I) == 1 + I);
assert(copysign( 1 + I, 1 - I) == 1 - I);
assert(copysign( 1 + I, -1 + I) == -1 + I);
assert(copysign( 1 + I, -1 - I) == -1 - I);
assert(copysign( 1 - I, 1 + I) == 1 + I);
assert(copysign( 1 - I, 1 - I) == 1 - I);
assert(copysign( 1 - I, -1 + I) == -1 + I);
assert(copysign( 1 - I, -1 - I) == -1 - I);
assert(copysign(-1 + I, 1 + I) == 1 + I);
assert(copysign(-1 + I, 1 - I) == 1 - I);
assert(copysign(-1 + I, -1 + I) == -1 + I);
assert(copysign(-1 + I, -1 - I) == -1 - I);
assert(copysign(-1 - I, 1 + I) == 1 + I);
assert(copysign(-1 - I, 1 - I) == 1 - I);
assert(copysign(-1 - I, -1 + I) == -1 + I);
assert(copysign(-1 - I, -1 - I) == -1 - I);
}
auto poly(T1, T2)(T1 x, in T2[] A)
if ((isComplex!T1 || isFloatingPoint!T1) && (isComplex!T2 || isFloatingPoint!T2) && (isComplex!T1 || isComplex!T2))
{
ptrdiff_t i = A.length - 1;
typeof(x*A[0]) r = A[i];
while (--i >= 0)
{
r *= x;
r += A[i];
}
return r;
}
unittest {
assert(poly(I, [1f,2,3,4,5]).approxEqual(1 + 2*I - 3 - 4*I + 5));
}
auto isNaN(T)(Complex!T x) {
return .isNaN(x.re) || .isNaN(x.im);
}
unittest {
assert(isNaN(I / 0));
assert(isNaN((1 + 0*I) / 0));
assert(!isNaN(I));
}
auto isInfinity(T)(Complex!T x) {
return .isInfinity(x.re) || .isInfinity(x.im) && !isNaN(x);
}
unittest {
assert(isInfinity((1 + 1*I) / 0));
assert(!isInfinity(1 + 1*I));
}
auto approxEqual(T1, T2)(T1 x, T2 y)
if ((isComplex!T1 || isFloatingPoint!T1) && (isComplex!T2 || isFloatingPoint!T2) && (isComplex!T1 || isComplex!T2))
{
alias approxEqual = std.math.approxEqual;
static if (isComplex!T1 && isComplex!T2) {
return approxEqual(x.re, y.re) && approxEqual(x.im, y.im);
} else static if (isComplex!T1) {
return approxEqual(x.re, y) && approxEqual(x.im, 0);
} else {
return approxEqual(x, y.re) && approxEqual(0, y.im);
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.