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

- "Mehrdad" <wfunction hotmail.com> Sep 25 2012
- "Mehrdad" <wfunction hotmail.com> Sep 25 2012
- Andrei Alexandrescu <SeeWebsiteForEmail erdani.org> Sep 25 2012
- "Mehrdad" <wfunction hotmail.com> Sep 25 2012
- "monarch_dodra" <monarchdodra gmail.com> Sep 25 2012
- "Mehrdad" <wfunction hotmail.com> Sep 25 2012
- "Mehrdad" <wfunction hotmail.com> Sep 25 2012
- "jerro" <a a.com> Sep 25 2012
- "jerro" <a a.com> Sep 25 2012
- "jerro" <a a.com> Sep 25 2012
- "Mehrdad" <wfunction hotmail.com> Sep 25 2012

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

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 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 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

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

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 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 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

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

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

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

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. :PI'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