Comment #0 by bearophile_hugs — 2013-04-29T15:50:44Z
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));
}