D issues are now tracked on GitHub. This Bugzilla instance remains as a read-only archive.
Issue 10010 - Some small ideas for std.complex
Summary: Some small ideas for std.complex
Status: RESOLVED WONTFIX
Alias: None
Product: D
Classification: Unclassified
Component: phobos (show other issues)
Version: D2
Hardware: All All
: P2 normal
Assignee: No Owner
URL:
Keywords:
Depends on:
Blocks:
 
Reported: 2013-04-29 15:50 UTC by bearophile_hugs
Modified: 2020-03-21 03:56 UTC (History)
0 users

See Also:


Attachments

Note You need to log in before you can comment on or make changes to this issue.
Description bearophile_hugs 2013-04-29 15:50:44 UTC
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));
}
Comment 1 basile-z 2017-09-15 05:32:02 UTC
So small that nobody cared.