## digitalmars.D - std.range should support recursion (Was: One-line FFT, nice!)

• Mehrdad (22/22) Sep 25 2012 I thought I'd take the opportunity to point this out:
• Mehrdad (2/3) Sep 25 2012 new types**
• Andrei Alexandrescu (3/6) Sep 25 2012 Ah, better now. Still it would be great to explain it more :o).
• Mehrdad (14/22) Sep 25 2012 Haha ok. :) I mean like, essentially, these need to work:
• monarch_dodra (25/42) Sep 25 2012 I can't comment on the rest of your points, but stride and take
• Mehrdad (44/69) Sep 25 2012 I just wrote down the assert's on the fly, actually.
• jerro (3/41) Sep 25 2012 One possible reason for this error could be that you are
• jerro (22/73) Sep 25 2012 I'd like to add that I don't think you can make this work the way
• jerro (6/17) Sep 25 2012 I was wrong about the complexity. Because each element of the
• Mehrdad (7/36) Sep 25 2012 Great point, I hadn't really thought about the laziness before --
• Andrei Alexandrescu (5/21) Sep 25 2012 I think all of the above are doable and useful. Please file a bug report...
• Mehrdad (5/10) Sep 25 2012 Sure! I'll file them as soon as I know they're bugs --
• Andrei Alexandrescu (5/9) Sep 25 2012 I'm not sure I understand this, and it seems I should. Could you please
I thought I'd take the opportunity to point this out:

The one-line FFT in D is pretty inefficient because it allocates
memory.

If std.range supported recursion (i.e. by providing a different
implementation for ranges that can be implemented without
creating new times, i.e. Stride of Stride == Stride), then it
would make the library a lot more usable and less bloated.

My one-line FFT illustrates it perfectly:

import std.algorithm;
import std.math;
import std.range;
typeof(R.init.stride(0)) dft(R)(R v)
{
return v.length > 1 ?
(p => chain(map!(q => q[0] + q[1])(p), map!(q => q[0] -
q[1])(p)))
(zip(dft(v.stride(2)),
map!(p => p[1] * expi(p[0] * -2 * PI / v.length))
(zip(iota(v.length / 2), dft(v.drop(1).stride(2)))))) : v;
}
void main() { dft([1.0, 2, 3]); }

Side note: the error messages are also hard to read.
Sep 25 2012
On Tuesday, 25 September 2012 at 08:21:39 UTC, Mehrdad wrote:
without creating new times
new types**
Sep 25 2012
Andrei Alexandrescu <SeeWebsiteForEmail erdani.org> writes:
On 9/25/12 4:23 AM, Mehrdad wrote:
On Tuesday, 25 September 2012 at 08:21:39 UTC, Mehrdad wrote:
without creating new times
new types**
Ah, better now. Still it would be great to explain it more :o). Andrei
Sep 25 2012
On Tuesday, 25 September 2012 at 13:34:28 UTC, Andrei
Alexandrescu wrote:
On 9/25/12 4:23 AM, Mehrdad wrote:
On Tuesday, 25 September 2012 at 08:21:39 UTC, Mehrdad wrote:
without creating new times
new types**
Ah, better now. Still it would be great to explain it more :o). Andrei
Haha ok. :) I mean like, essentially, these need to work: assert(is(typeof(foo.stride(1)) == typeof(foo.stride(2).stride(3)))); assert(is(typeof(foo.drop(1)) == typeof(foo.drop(2).drop(3)))); assert(is(typeof(foo.take(1)) == typeof(foo.take(2).take(3)))); otherwise recursion with these ranges is impossible. The FFT example took the odd- and even-indexed numbers with stride(), but it couldn't recursively do this because the type system prevented it from doing so. So I was forced to copy the array unnecessarily every time. Also, foo should be implicitly convertible to typeof(foo.stride(1)), which also makes recursion easier.
Sep 25 2012
"monarch_dodra" <monarchdodra gmail.com> writes:
On Tuesday, 25 September 2012 at 15:41:42 UTC, Mehrdad wrote:
On Tuesday, 25 September 2012 at 13:34:28 UTC, Andrei
Alexandrescu wrote:
On 9/25/12 4:23 AM, Mehrdad wrote:
On Tuesday, 25 September 2012 at 08:21:39 UTC, Mehrdad wrote:
without creating new times
new types**
Ah, better now. Still it would be great to explain it more :o). Andrei
Haha ok. :) I mean like, essentially, these need to work: assert(is(typeof(foo.stride(1)) == typeof(foo.stride(2).stride(3)))); assert(is(typeof(foo.drop(1)) == typeof(foo.drop(2).drop(3)))); assert(is(typeof(foo.take(1)) == typeof(foo.take(2).take(3))));
I can't comment on the rest of your points, but stride and take DO check for type recursivity, and drop always returns the same type as input anyways. Failure of ANY of these asserts is a bug. What where your inputs? //----------------------- import std.range; struct S { enum empty = false; void popFront(){}; property int front(){return 1;} } void main() { S foo; static assert(is(typeof(foo.stride(1)) == typeof(foo.stride(2).stride(3)))); static assert(is(typeof(foo.drop(1)) == typeof(foo.drop(2).drop(3)))); static assert(is(typeof(foo) == typeof(foo.drop(2)))); //Or this static assert(is(typeof(foo.take(1)) == typeof(foo.take(2).take(3)))); } //-----------------------
Sep 25 2012
On Tuesday, 25 September 2012 at 16:03:22 UTC, monarch_dodra
wrote:
On Tuesday, 25 September 2012 at 15:41:42 UTC, Mehrdad wrote:
On Tuesday, 25 September 2012 at 13:34:28 UTC, Andrei
Alexandrescu wrote:
On 9/25/12 4:23 AM, Mehrdad wrote:
On Tuesday, 25 September 2012 at 08:21:39 UTC, Mehrdad wrote:
without creating new times
new types**
Ah, better now. Still it would be great to explain it more :o). Andrei
Haha ok. :) I mean like, essentially, these need to work: assert(is(typeof(foo.stride(1)) == typeof(foo.stride(2).stride(3)))); assert(is(typeof(foo.drop(1)) == typeof(foo.drop(2).drop(3)))); assert(is(typeof(foo.take(1)) == typeof(foo.take(2).take(3))));
I can't comment on the rest of your points, but stride and take DO check for type recursivity, and drop always returns the same type as input anyways. Failure of ANY of these asserts is a bug. What where your inputs?
I just wrote down the assert's on the fly, actually. Maybe I'm just misinterpreting the error then, and the problem is somewhere else? The code I was trying to compile is this (sorry it's ugly): import std.algorithm, std.math, std.range; typeof(R.init.stride(0)) dft(R)(R v) { return v.length <= 1 ? v.stride(1) : (p => chain(map!(q => q[0] + q[1])(p), map!(q => q[0] - q[1])(p))) (zip(dft(v.stride(2)), map!(p => p[1] * expi(p[0] * -2 * PI / v.length)) (zip(iota(v.length / 2), dft(v.drop(1).stride(2)))))); } void main() { dft([1.0, 2, 3]); } Which gives the following error: Test.d(5): Error: incompatible types for ((stride(v,1u)) ? ((*delegate system Result(Zip!(Result,MapResult!(__lambda8,Zip!(Result,Result))) p) { return chain(map(p),map(p)); } )(zip(dft(stride(v,2u)),map(zip(iota(v.length() / 2u),dft(stride(drop(v,1u),2u)))))))): 'Result' and 'Result' Test.d(10): Error: template instance Test.dft!(Result) error instantiating Test.d:19: instantiated from here: dft!(double[]) Test.d(5): Error: incompatible types for ((stride(v,1u)) ? ((*delegate system Result(Zip!(Result,MapResult!(__lambda8,Zip!(Result,Result))) p) { return chain(map(p),map(p)); } )(zip(dft(stride(v,2u)),map(zip(iota(v.length / 2u),dft(stride(drop(v,1u),2u)))))))): 'Result' and 'Result' Test.d(19): Error: template instance Test.dft!(double[]) error instantiating How should I interpret it? Thanks!
Sep 25 2012
"jerro" <a a.com> writes:
import std.algorithm, std.math, std.range;
typeof(R.init.stride(0)) dft(R)(R v)
{
return v.length <= 1
? v.stride(1)
: (p =>
chain(map!(q => q[0] + q[1])(p),
map!(q => q[0] - q[1])(p)))
(zip(dft(v.stride(2)),
map!(p => p[1] * expi(p[0] * -2 * PI / v.length))
(zip(iota(v.length / 2),
dft(v.drop(1).stride(2))))));
}
void main() { dft([1.0, 2, 3]); }

Which gives the following error:

Test.d(5): Error: incompatible types for ((stride(v,1u)) ?
((*delegate  system
Result(Zip!(Result,MapResult!(__lambda8,Zip!(Result,Result))) p)
{
return chain(map(p),map(p));
}
)(zip(dft(stride(v,2u)),map(zip(iota(v.length() /
2u),dft(stride(drop(v,1u),2u)))))))): 'Result' and 'Result'
Test.d(10): Error: template instance Test.dft!(Result) error
instantiating
Test.d:19: instantiated from here: dft!(double[])
Test.d(5): Error: incompatible types for ((stride(v,1u)) ?
((*delegate  system
Result(Zip!(Result,MapResult!(__lambda8,Zip!(Result,Result))) p)
{
return chain(map(p),map(p));
}
)(zip(dft(stride(v,2u)),map(zip(iota(v.length /
2u),dft(stride(drop(v,1u),2u)))))))): 'Result' and 'Result'
Test.d(19): Error: template instance Test.dft!(double[]) error
instantiating

How should I interpret it?

Thanks!
One possible reason for this error could be that you are returning a result of chain() if length is larger than 1 and a result of stride() otherwise.
Sep 25 2012
"jerro" <a a.com> writes:
On Tuesday, 25 September 2012 at 17:48:49 UTC, jerro wrote:
import std.algorithm, std.math, std.range;
typeof(R.init.stride(0)) dft(R)(R v)
{
return v.length <= 1
? v.stride(1)
: (p =>
chain(map!(q => q[0] + q[1])(p),
map!(q => q[0] - q[1])(p)))
(zip(dft(v.stride(2)),
map!(p => p[1] * expi(p[0] * -2 * PI / v.length))
(zip(iota(v.length / 2),
dft(v.drop(1).stride(2))))));
}
void main() { dft([1.0, 2, 3]); }

Which gives the following error:

Test.d(5): Error: incompatible types for ((stride(v,1u)) ?
((*delegate  system
Result(Zip!(Result,MapResult!(__lambda8,Zip!(Result,Result)))
p)
{
return chain(map(p),map(p));
}
)(zip(dft(stride(v,2u)),map(zip(iota(v.length() /
2u),dft(stride(drop(v,1u),2u)))))))): 'Result' and 'Result'
Test.d(10): Error: template instance Test.dft!(Result) error
instantiating
Test.d:19: instantiated from here: dft!(double[])
Test.d(5): Error: incompatible types for ((stride(v,1u)) ?
((*delegate  system
Result(Zip!(Result,MapResult!(__lambda8,Zip!(Result,Result)))
p)
{
return chain(map(p),map(p));
}
)(zip(dft(stride(v,2u)),map(zip(iota(v.length /
2u),dft(stride(drop(v,1u),2u)))))))): 'Result' and 'Result'
Test.d(19): Error: template instance Test.dft!(double[]) error
instantiating

How should I interpret it?

Thanks!
One possible reason for this error could be that you are returning a result of chain() if length is larger than 1 and a result of stride() otherwise.
I'd like to add that I don't think you can make this work the way you meant it to. The problem is that you return a chain when the length is 2, a chain of chains when the length is 4, and so on. What fft of length n would actually need to return is a binary tree of ranges with n leaves. The size of memory needed for the value returned from fft therefore depends on the length of the range it was given as an argument. So you need to either use heap allocated memory for the return type, or the return type needs to depend on the parameter range's length, which would mean that the parameter range's length needs to be a template parameter. I think there is one thing in this code that will hurt performance much, much, more than allocations. This code will compute elements of the result lazily. So each time you want to read an element from the resulting range, O(log(n)) functions passed to map() will need to be computed. The problem is that each of those functions computes sine and cosine, so sine and cosine need to be computed O(log(n)) times for each element. To get all n elements, you will need to compute them O(n log(n)). Because computing sine and cosine is about two orders of magnitude slower than multiplication and division, this will be very slow.
Sep 25 2012
"jerro" <a a.com> writes:
I think there is one thing in this code that will hurt
performance much, much, more than allocations. This code will
compute elements of the result lazily. So each time you want to
read an element from the resulting range, O(log(n)) functions
passed to map() will need to be computed. The problem is that
each of those functions computes sine and cosine, so sine and
cosine need to be computed O(log(n)) times for each element. To
get all n elements, you will need to compute them O(n log(n)).
Because computing sine and cosine is about two orders of
magnitude slower than multiplication and division, this will be
very slow.
I was wrong about the complexity. Because each element of the result depends on all the elements of the argument range, you actually need O(n) function calls to compute each element of the result and O(n*n) function calls(and sine and cosine computations) to compute all of them. You would need to use memoization to get reasonable complexity.
Sep 25 2012
On Tuesday, 25 September 2012 at 18:33:45 UTC, jerro wrote:
I think there is one thing in this code that will hurt
performance much, much, more than allocations. This code will
compute elements of the result lazily. So each time you want
to read an element from the resulting range, O(log(n))
functions passed to map() will need to be computed. The
problem is that each of those functions computes sine and
cosine, so sine and cosine need to be computed O(log(n)) times
for each element. To get all n elements, you will need to
compute them O(n log(n)). Because computing sine and cosine is
about two orders of magnitude slower than multiplication and
division, this will be very slow.
I was wrong about the complexity. Because each element of the result depends on all the elements of the argument range, you actually need O(n) function calls to compute each element of the result and O(n*n) function calls(and sine and cosine computations) to compute all of them. You would need to use memoization to get reasonable complexity.
Great point, I hadn't really thought about the laziness before -- I was doing this in Python and just thought it might be fun to translate it to D. :P
I'd like to add that I don't think you can make this work the
way you meant it to. The problem is that you return a chain
when the length is 2, a chain of chains when the length is 4,
and so on. What fft of length n would actually need to return
is a binary tree of ranges with n leaves. The size of memory
needed for the value returned from fft therefore depends on the
length of the range it was given as an argument. So you need to
either use heap allocated memory for the return type, or the
return type needs to depend on the parameter range's length,
which would mean that the parameter range's length needs to be
a template parameter.
Ohh huh... sounds like you're right, lemme think about it a bit more though. Thanks! :)
Sep 25 2012
Andrei Alexandrescu <SeeWebsiteForEmail erdani.org> writes:
On 9/25/12 11:42 AM, Mehrdad wrote:
On Tuesday, 25 September 2012 at 13:34:28 UTC, Andrei Alexandrescu wrote:
On 9/25/12 4:23 AM, Mehrdad wrote:
On Tuesday, 25 September 2012 at 08:21:39 UTC, Mehrdad wrote:
without creating new times
new types**
Ah, better now. Still it would be great to explain it more :o). Andrei
Haha ok. :) I mean like, essentially, these need to work: assert(is(typeof(foo.stride(1)) == typeof(foo.stride(2).stride(3)))); assert(is(typeof(foo.drop(1)) == typeof(foo.drop(2).drop(3)))); assert(is(typeof(foo.take(1)) == typeof(foo.take(2).take(3)))); otherwise recursion with these ranges is impossible.
I think all of the above are doable and useful. Please file a bug report containing these and any others you could reasonably think of. Thanks! Andrei
Sep 25 2012
On Tuesday, 25 September 2012 at 16:31:15 UTC, Andrei
Alexandrescu wrote:
I think all of the above are doable and useful. Please file a
bug report containing these and any others you could reasonably
think of.

Thanks!

Andrei
Sure! I'll file them as soon as I know they're bugs -- monarch_dodra's comment makes me think it might just be my misinterpretation of the error.
Sep 25 2012
Andrei Alexandrescu <SeeWebsiteForEmail erdani.org> writes:
On 9/25/12 4:22 AM, Mehrdad wrote:
If std.range supported recursion (i.e. by providing a different
implementation for ranges that can be implemented without creating new
times, i.e. Stride of Stride == Stride), then it would make the library
a lot more usable and less bloated.
I'm not sure I understand this, and it seems I should. Could you please explain (perhaps with a simpler example than FFT)? Thanks, Andrei
Sep 25 2012