Missing math functions in std.complex

complexfunctions.d

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);
    }
}