Wolfgang Thallergcc16 merge934c3885 months agoHistory
1972 lines · 53.3 KB · D
RawDownload
// 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.     $(LREF Complex) will eventually    $(DDLINK deprecate, Deprecated Features, replace)    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 = &plusmn;        NAN = $(RED NAN)        INFIN = &infin;        PI = &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)*/module std.complex; import std.traits; /** Helper function that returns a complex number with the specified    real and imaginary parts.     Params:        R = (template parameter) type of real part of complex number        I = (template parameter) type of imaginary part of complex number         re = real part of complex number to be constructed        im = (optional) imaginary part of complex number, 0 if omitted.     Returns:        `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)(const R re)  @safe pure nothrow @nogcif (is(R : double)){    static if (isFloatingPoint!R)        return Complex!R(re, 0);    else        return Complex!double(re, 0);} /// dittoauto complex(R, I)(const R re, const I im)  @safe pure nothrow @nogcif (is(R : double) && is(I : double)){    static if (isFloatingPoint!R || isFloatingPoint!I)        return Complex!(CommonType!(R, I))(re, im);    else        return Complex!double(re, im);} ///@safe pure nothrow unittest{    auto a = complex(1.0);    static assert(is(typeof(a) == Complex!double));    assert(a.re == 1.0);    assert(a.im == 0.0);     auto b = complex(2.0L);    static assert(is(typeof(b) == Complex!real));    assert(b.re == 2.0L);    assert(b.im == 0.0L);     auto c = complex(1.0, 2.0);    static assert(is(typeof(c) == Complex!double));    assert(c.re == 1.0);    assert(c.im == 2.0);     auto d = complex(3.0, 4.0L);    static assert(is(typeof(d) == Complex!real));    assert(d.re == 3.0);    assert(d.im == 4.0L);     auto e = complex(1);    static assert(is(typeof(e) == Complex!double));    assert(e.re == 1);    assert(e.im == 0);     auto f = complex(1L, 2);    static assert(is(typeof(f) == Complex!double));    assert(f.re == 1L);    assert(f.im == 2);     auto g = complex(3, 4.0L);    static assert(is(typeof(g) == Complex!real));    assert(g.re == 3);    assert(g.im == 4.0L);}  /** 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.spec : FormatSpec;    import std.range.primitives : isOutputRange;     /** The real part of the number. */    T re;     /** The imaginary part of the number. */    T im;     /** Converts the complex number to a string representation.     The second form of this function is usually not called directly;    instead, it is used via $(REF format, std,string), as shown in the examples    below.  Supported format characters are 'e', 'f', 'g', 'a', and 's'.     See the $(MREF std, format) and $(REF format, std,string)    documentation for more information.    */    string toString() const @safe /* TODO: pure nothrow */    {        import std.exception : assumeUnique;        char[] buf;        buf.reserve(100);        auto fmt = FormatSpec!char("%s");        toString((const(char)[] s) { buf ~= s; }, fmt);        static trustedAssumeUnique(T)(T t) @trusted { return assumeUnique(t); }        return trustedAssumeUnique(buf);    }     static if (is(T == double))    ///    @safe unittest    {        auto c = complex(1.2, 3.4);         // Vanilla toString formatting:        assert(c.toString() == "1.2+3.4i");         // Formatting with std.string.format specs: the precision and width        // specifiers apply to both the real and imaginary parts of the        // complex number.        import std.format : format;        assert(format("%.2f", c)  == "1.20+3.40i");        assert(format("%4.1f", c) == " 1.2+ 3.4i");    }     /// ditto    void toString(Writer, Char)(scope Writer w, scope const ref FormatSpec!Char formatSpec) const    if (isOutputRange!(Writer, const(Char)[]))    {        import std.format.write : formatValue;        import std.math.traits : signbit;        import std.range.primitives : put;        formatValue(w, re, formatSpec);        if (signbit(im) == 0)           put(w, "+");        formatValue(w, im, formatSpec);        put(w, "i");    } @safe pure nothrow @nogc:     /** Construct a complex number with the specified real and    imaginary parts. In the case where a single argument is passed    that is not complex, the imaginary part of the result will be    zero.    */    this(R : T)(Complex!R z)    {        re = z.re;        im = z.im;    }     /// ditto    this(Rx : T, Ry : T)(const Rx x, const Ry y)    {        re = x;        im = y;    }     /// ditto    this(R : T)(const R r)    {        re = r;        im = 0;    }     // ASSIGNMENT OPERATORS     // this = complex    ref Complex opAssign(R : T)(Complex!R z)    {        re = z.re;        im = z.im;        return this;    }     // this = numeric    ref Complex opAssign(R : T)(const R r)    {        re = r;        im = 0;        return this;    }     // COMPARISON OPERATORS     // this == complex    bool opEquals(R : T)(Complex!R z) const    {        return re == z.re && im == z.im;    }     // this == numeric    bool opEquals(R : T)(const R r) const    {        return re == r && im == 0;    }     // UNARY OPERATORS     // +complex    Complex opUnary(string op)() const    if (op == "+")    {        return this;    }     // -complex    Complex opUnary(string op)() const    if (op == "-")    {        return Complex(-re, -im);    }     // BINARY OPERATORS     // complex op complex    Complex!(CommonType!(T,R)) opBinary(string op, R)(Complex!R z) const    {        alias C = typeof(return);        auto w = C(this.re, this.im);        return w.opOpAssign!(op)(z);    }     // complex op numeric    Complex!(CommonType!(T,R)) opBinary(string op, R)(const R r) const    if (isNumeric!R)    {        alias C = typeof(return);        auto w = C(this.re, this.im);        return w.opOpAssign!(op)(r);    }     // numeric + complex,  numeric * complex    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)(const R r) const    if (op == "-" && isNumeric!R)    {        return Complex(r - re, -im);    }     // numeric / complex    Complex!(CommonType!(T, R)) opBinaryRight(string op, R)(const R r) const    if (op == "/" && isNumeric!R)    {        version (FastMath)        {            // 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        {            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;        }    }     // numeric ^^ complex    Complex!(CommonType!(T, R)) opBinaryRight(string op, R)(const R lhs) const    if (op == "^^" && isNumeric!R)    {        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)        {            // r = lhs            // theta = 0            ab = lhs ^^ this.re;            ar = log(lhs) * this.im;        }        else        {            // r = -lhs            // theta = PI            ab = (-lhs) ^^ this.re * exp(-PI * this.im);            ar = PI * this.re + log(-lhs) * this.im;        }         return typeof(return)(ab * cos(ar), ab * sin(ar));    }     // OP-ASSIGN OPERATORS     // complex += complex,  complex -= complex    ref Complex opOpAssign(string op, C)(const C z)    if ((op == "+" || op == "-") && is(C R == Complex!R))    {        mixin ("re "~op~"= z.re;");        mixin ("im "~op~"= z.im;");        return this;    }     // complex *= complex    ref Complex opOpAssign(string op, C)(const C z)    if (op == "*" && is(C R == Complex!R))    {        auto temp = re*z.re - im*z.im;        im = im*z.re + re*z.im;        re = temp;        return this;    }     // complex /= complex    ref Complex opOpAssign(string op, C)(const C z)    if (op == "/" && is(C R == Complex!R))    {        version (FastMath)        {            // 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        {            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;        }    }     // complex ^^= complex    ref Complex opOpAssign(string op, C)(const C z)    if (op == "^^" && is(C R == Complex!R))    {        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);        immutable ar = t*z.re + log(r)*z.im;         re = ab*cos(ar);        im = ab*sin(ar);        return this;    }     // complex += numeric,  complex -= numeric    ref Complex opOpAssign(string op, U : T)(const U a)    if (op == "+" || op == "-")    {        mixin ("re "~op~"= a;");        return this;    }     // complex *= numeric,  complex /= numeric    ref Complex opOpAssign(string op, U : T)(const U a)    if (op == "*" || op == "/")    {        mixin ("re "~op~"= a;");        mixin ("im "~op~"= a;");        return this;    }     // complex ^^= real    ref Complex opOpAssign(string op, R)(const R r)    if (op == "^^" && isFloatingPoint!R)    {        import core.math : cos, sin;        immutable ab = abs(this)^^r;        immutable ar = arg(this)*r;        re = ab*cos(ar);        im = ab*sin(ar);        return this;    }     // complex ^^= int    ref Complex opOpAssign(string op, U)(const U i)    if (op == "^^" && isIntegral!U)    {        switch (i)        {        case 0:            re = 1.0;            im = 0.0;            break;        case 1:            // identity; do nothing            break;        case 2:            this *= this;            break;        case 3:            auto z = this;            this *= z;            this *= z;            break;        default:            this ^^= cast(real) i;        }        return this;    }     /** Returns a complex number instance that correponds in size and in ABI        to the associated C compiler's `_Complex` type.     */    auto toNative()    {        import core.stdc.config : c_complex_float, c_complex_double, c_complex_real;        static if (is(T == float))            return c_complex_float(re, im);        else static if (is(T == double))            return c_complex_double(re, im);        else            return c_complex_real(re, im);    }} @safe pure nothrow unittest{    import std.complex;    static import core.math;

Lines 1–500 of 1972.

Keyboard shortcuts

General

Show keyboard shortcuts
?
Focus the current search
/

Issues and pull requests

Open filters
f

Issues, pull requests and Inbox

Select next item
j↓
Select previous item
k↑
Open selected item
Enter
Dismiss filters or clear selection
Esc

Repository code

Go to file
t