[Issue 10010] New: Some small ideas for std.complex

[d-bugmail at puremagic.com]

http://d.puremagic.com/issues/show_bug.cgi?id=10010

           Summary: Some small ideas for std.complex
           Product: D
           Version: D2
          Platform: All
        OS/Version: All
            Status: NEW
          Severity: normal
          Priority: P2
         Component: Phobos
        AssignedTo: nobody at puremagic.com
        ReportedBy: bearophile_hugs at eml.cc


--- Comment #0 from bearophile_hugs at eml.cc 2013-04-29 15:50:44 PDT ---
1) I'd like std.complex.expi to have this signature:

pure nothrow @trusted Complex!T expi(T)(T y); 

In computations it's useful to have back the type of complex you give it,
otherwise later you will probably have to cast some complex type to make things
fit together.

As example look at this program, it uses Complex!real everywhere instead of
Complex!double because expi returns a Complex!real:


import std.stdio, std.algorithm, std.range, std.complex;
import std.math: PI;

Complex!real[] fft(Complex!real[] x) /*pure nothrow*/ {
    immutable N = x.length;
    if (N <= 1) return x;
    const ev = x.stride(2).array.fft;
    const od = x[1 .. $].stride(2).array.fft;
    auto l = iota(N / 2).map!(k => ev[k] + expi(-2*PI*k/N) * od[k]);
    auto r = iota(N / 2).map!(k => ev[k] - expi(-2*PI*k/N) * od[k]);
    return l.chain(r).array;
}

void main() {
    [1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0]
    .map!(r => Complex!real(r, 0))
    .array
    .fft
    .writeln;
}


(In this program PI was imported directly because if you just import std.math
and std.complex, the two expi clash!).

One current workaround is to cast:


Complex!double[] fft(Complex!double[] x) /*pure nothrow*/ {
    immutable N = x.length;
    if (N <= 1) return x;
    const ev = x.stride(2).array.fft;
    const od = x[1 .. $].stride(2).array.fft;
    auto l = iota(N / 2).map!(k => ev[k] + cast(Complex!double)expi(-2*PI*k/N)
* od[k]);
    auto r = iota(N / 2).map!(k => ev[k] - cast(Complex!double)expi(-2*PI*k/N)
* od[k]);
    return l.chain(r).array;
}


The modified expi allows to write (I have had to define w because PI is a
real):


import std.stdio, std.algorithm, std.range, std.math, std.complex;

Complex!T expi(T)(T y) @trusted pure nothrow
{
    return Complex!T(std.math.cos(y), std.math.sin(y));
}

auto fft(T)(in T[] x) /*pure nothrow*/ {
    immutable N = x.length;
    if (N <= 1) return x;
    const ev = x.stride(2).array.fft;
    const od = x[1 .. $].stride(2).array.fft;
    immutable double w = -2 * PI / N;
    auto l = iota(N / 2).map!(k => ev[k] + expi(k * w) * od[k]);
    auto r = iota(N / 2).map!(k => ev[k] - expi(k * w) * od[k]);
    return l.chain(r).array;
}

void main() {
    [1.0, 1, 1, 1, 0, 0, 0, 0].map!complex.array.fft.writeln;
}

- - - - - - - - - - - - -

2) Maybe it's handy to have an alias similar to this defined inside
std.complex, to be used as "standard" complex value in code:

alias Complexd = Complex!double;

- - - - - - - - - - - - -

3) The documentation of std.complex.expi says:


Unlike std.math.expi, which uses the x87 fsincos instruction when possible,
this function is no faster than calculating cos(y) and sin(y) separately.


But probably on GDC and maybe LDC the compiler is able to recognize the nearby
calls to sin and cos in code like this, and implement it with just one call to
sincos, so I suggest to remove that comment from the std.complex docs:


Complex!real expi(real y)  @trusted pure nothrow
{
    return Complex!real(std.math.cos(y), std.math.sin(y));
}

-- 
Configure issuemail: http://d.puremagic.com/issues/userprefs.cgi?tab=email
------- You are receiving this mail because: -------