From c0a3ea1b89bf88e80d5624cd33a8594debe5dd5e Mon Sep 17 00:00:00 2001 From: Georgy Markov Date: Thu, 5 Jun 2025 23:46:06 +0300 Subject: [PATCH 01/31] =?UTF-8?q?Fix=20dlang#10796=20-=20FFT=20doesn?= =?UTF-8?q?=E2=80=99t=20work=20with=20user-defined=20complex=20types?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- std/numeric.d | 48 +++++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 43 insertions(+), 5 deletions(-) diff --git a/std/numeric.d b/std/numeric.d index 918984fd52e..7815595bdb2 100644 --- a/std/numeric.d +++ b/std/numeric.d @@ -3733,7 +3733,14 @@ public: } else if (range.length == 1) { - buf[0] = range[0]; + static if (isComplexLike!(ElementType!R)) { + buf[0].re = range[0].re; + buf[0].im = range[0].im; + } + else { + buf[0].re = range[0]; + buf[0].im = 0; + } return; } else if (range.length == 2) @@ -3931,6 +3938,26 @@ void inverseFft(Ret, R)(R range, Ret buf) assert(isClose(twoInv[1].im, 0, 0.0, 1e-10)); } +// https://github.com/dlang/phobos/issues/10796 +@system unittest +{ + import std.algorithm; + import std.range; + static struct C { float re, im; } // User-defined complex + + float[8] arr = [1,2,3,4,5,6,7,8]; + C[8] fft1; + fft(arr[], fft1[]); + assert(isClose(fft1[].map!"a.re", + [36.0, -4, -4, -4, -4, -4, -4, -4], 1e-4)); + assert(isClose(fft1[].map!"a.im", + [0, 9.6568, 4, 1.6568, 0, -1.6568, -4, -9.6568], 1e-4)); + + auto inv = inverseFft(fft1[]); + assert(isClose(inv[].map!"a.re", arr[], 1e-6)); + assert(inv[].map!"a.im".maxElement < 1e-10); +} + // Swaps the real and imaginary parts of a complex number. This is useful // for inverse FFTs. C swapRealImag(C)(C input) @@ -4112,11 +4139,22 @@ struct Stride(R) // using a generic slow DFT. This seems to be the best base case. (Size 1 // can be coded inline as buf[0] = range[0]). void slowFourier2(Ret, R)(R range, Ret buf) +if (isComplexLike!(ElementType!Ret)) +in (range.length == 2) +in (buf.length == 2) { - assert(range.length == 2); - assert(buf.length == 2); - buf[0] = range[0] + range[1]; - buf[1] = range[0] - range[1]; + static if (isComplexLike!(ElementType!R)) { + buf[0].re = range[0].re + range[1].re; + buf[0].im = range[0].im + range[1].im; + buf[1].re = range[0].re - range[1].re; + buf[1].im = range[0].im - range[1].im; + } + else { + buf[0].re = range[0] + range[1]; + buf[0].im = 0; + buf[1].re = range[0] - range[1]; + buf[1].im = 0; + } } // Hard-coded base case for FFT of size 4. Doesn't work as well as the size From c1789b754b1239ac6be73b1785855d3492923d73 Mon Sep 17 00:00:00 2001 From: Georgy Markov Date: Fri, 6 Jun 2025 21:59:23 +0300 Subject: [PATCH 02/31] Style --- std/numeric.d | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/std/numeric.d b/std/numeric.d index 7815595bdb2..9330fa502f3 100644 --- a/std/numeric.d +++ b/std/numeric.d @@ -3733,11 +3733,13 @@ public: } else if (range.length == 1) { - static if (isComplexLike!(ElementType!R)) { + static if (isComplexLike!(ElementType!R)) + { buf[0].re = range[0].re; buf[0].im = range[0].im; } - else { + else + { buf[0].re = range[0]; buf[0].im = 0; } @@ -4143,13 +4145,15 @@ if (isComplexLike!(ElementType!Ret)) in (range.length == 2) in (buf.length == 2) { - static if (isComplexLike!(ElementType!R)) { + static if (isComplexLike!(ElementType!R)) + { buf[0].re = range[0].re + range[1].re; buf[0].im = range[0].im + range[1].im; buf[1].re = range[0].re - range[1].re; buf[1].im = range[0].im - range[1].im; } - else { + else + { buf[0].re = range[0] + range[1]; buf[0].im = 0; buf[1].re = range[0] - range[1]; From 663dd6104081a1276760ae3dbdfb36bb305bfa4c Mon Sep 17 00:00:00 2001 From: Georgy Markov Date: Fri, 6 Jun 2025 22:17:20 +0300 Subject: [PATCH 03/31] Branch by `isNumeric` instead of `isComplexLike` --- std/numeric.d | 32 ++++++++++++++++---------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/std/numeric.d b/std/numeric.d index 9330fa502f3..345d3b7a87b 100644 --- a/std/numeric.d +++ b/std/numeric.d @@ -3733,16 +3733,16 @@ public: } else if (range.length == 1) { - static if (isComplexLike!(ElementType!R)) - { - buf[0].re = range[0].re; - buf[0].im = range[0].im; - } - else - { + static if (isNumeric!(ElementType!R)) + { buf[0].re = range[0]; buf[0].im = 0; } + else + { + buf[0].re = range[0].re; + buf[0].im = range[0].im; + } return; } else if (range.length == 2) @@ -4145,20 +4145,20 @@ if (isComplexLike!(ElementType!Ret)) in (range.length == 2) in (buf.length == 2) { - static if (isComplexLike!(ElementType!R)) - { - buf[0].re = range[0].re + range[1].re; - buf[0].im = range[0].im + range[1].im; - buf[1].re = range[0].re - range[1].re; - buf[1].im = range[0].im - range[1].im; - } - else - { + static if (isNumeric!(ElementType!R)) + { buf[0].re = range[0] + range[1]; buf[0].im = 0; buf[1].re = range[0] - range[1]; buf[1].im = 0; } + else + { + buf[0].re = range[0].re + range[1].re; + buf[0].im = range[0].im + range[1].im; + buf[1].re = range[0].re - range[1].re; + buf[1].im = range[0].im - range[1].im; + } } // Hard-coded base case for FFT of size 4. Doesn't work as well as the size From cef0b7cf4bdb35d13ae4324ceed18b899afc624c Mon Sep 17 00:00:00 2001 From: Georgy Markov Date: Sat, 7 Jun 2025 20:01:50 +0300 Subject: [PATCH 04/31] Reimplement FFT solver as a struct --- std/numeric.d | 234 +++++++++++++++++++++++++++++++++++++------------- 1 file changed, 176 insertions(+), 58 deletions(-) diff --git a/std/numeric.d b/std/numeric.d index 345d3b7a87b..7c09249c7d0 100644 --- a/std/numeric.d +++ b/std/numeric.d @@ -3372,13 +3372,8 @@ private alias lookup_t = float; */ final class Fft { - import core.bitop : bsf; - import std.algorithm.iteration : map; - import std.array : uninitializedArray; - private: - immutable lookup_t[][] negSinLookup; - + FFTImpl impl; void enforceSize(R)(R range) const { import std.conv : text; @@ -3386,6 +3381,109 @@ private: "FFT size mismatch. Expected ", size, ", got ", range.length)); } + // This constructor is used within this module for allocating the + // buffer space elsewhere besides the GC heap. It's definitely **NOT** + // part of the public API and definitely **IS** subject to change. + // + // Also, this is unsafe because the memSpace buffer will be cast + // to immutable. + // + // Public b/c of https://issues.dlang.org/show_bug.cgi?id=4636. + public this(lookup_t[] memSpace) + { + this.impl = FFTImpl(memSpace); + } + +public: + /**Create an `Fft` object for computing fast Fourier transforms of + * power of two sizes of `size` or smaller. `size` must be a + * power of two. + */ + this(size_t size) + { + this.impl = FFTImpl(size); + } + + @property size_t size() const => impl.size; + + /**Compute the Fourier transform of range using the $(BIGOH N log N) + * Cooley-Tukey Algorithm. `range` must be a random-access range with + * slicing and a length equal to `size` as provided at the construction of + * this object. The contents of range can be either numeric types, + * which will be interpreted as pure real values, or complex types with + * properties or members `.re` and `.im` that can be read. + * + * Note: Pure real FFTs are automatically detected and the relevant + * optimizations are performed. + * + * Returns: An array of complex numbers representing the transformed data in + * the frequency domain. + * + * Conventions: The exponent is negative and the factor is one, + * i.e., output[j] := sum[ exp(-2 PI i j k / N) input[k] ]. + */ + Complex!F[] fft(F = double, R)(R range) const + if (isFloatingPoint!F && isRandomAccessRange!R) + { + enforceSize(range); + return impl.fft(range); + } + + /**Same as the overload, but allows for the results to be stored in a user- + * provided buffer. The buffer must be of the same length as range, must be + * a random-access range, must have slicing, and must contain elements that are + * complex-like. This means that they must have a .re and a .im member or + * property that can be both read and written and are floating point numbers. + */ + void fft(Ret, R)(R range, Ret buf) const + if (isRandomAccessRange!Ret && isComplexLike!(ElementType!Ret) && hasSlicing!Ret) + in (buf.length == range.length) + { + enforceSize(range); + impl.fft(range, buf); + } + + /** + * Computes the inverse Fourier transform of a range. The range must be a + * random access range with slicing, have a length equal to the size + * provided at construction of this object, and contain elements that are + * either of type std.complex.Complex or have essentially + * the same compile-time interface. + * + * Returns: The time-domain signal. + * + * Conventions: The exponent is positive and the factor is 1/N, i.e., + * output[j] := (1 / N) sum[ exp(+2 PI i j k / N) input[k] ]. + */ + Complex!F[] inverseFft(F = double, R)(R range) const + if (isRandomAccessRange!R && isComplexLike!(ElementType!R) && isFloatingPoint!F) + { + enforceSize(range); + return impl.inverseFft(range); + } + + /** + * Inverse FFT that allows a user-supplied buffer to be provided. The buffer + * must be a random access range with slicing, and its elements + * must be some complex-like type. + */ + void inverseFft(Ret, R)(R range, Ret buf) const + if (isRandomAccessRange!Ret && isComplexLike!(ElementType!Ret) && hasSlicing!Ret) + { + enforceSize(range); + impl.inverseFft(range, buf); + } +} + +private struct FFTImpl +{ + import core.bitop : bsf; + import std.algorithm.iteration : map; + +private: + lookup_t* pStorage; + immutable lookup_t[] table; + void fftImpl(Ret, R)(Stride!R range, Ret buf) const in { @@ -3537,7 +3635,7 @@ private: do { immutable n = buf.length; - immutable localLookup = negSinLookup[bsf(n)]; + immutable localLookup = table[n .. 2*n]; assert(localLookup.length == n); immutable cosMask = n - 1; @@ -3599,37 +3697,62 @@ private: } } - // This constructor is used within this module for allocating the - // buffer space elsewhere besides the GC heap. It's definitely **NOT** - // part of the public API and definitely **IS** subject to change. - // - // Also, this is unsafe because the memSpace buffer will be cast - // to immutable. - // - // Public b/c of https://issues.dlang.org/show_bug.cgi?id=4636. - public this(lookup_t[] memSpace) - { - immutable size = memSpace.length / 2; +public: - /* Create a lookup table of all negative sine values at a resolution of - * size and all smaller power of two resolutions. This may seem - * inefficient, but having all the lookups be next to each other in - * memory at every level of iteration is a huge win performance-wise. - */ + /**Create an `FFT` object for computing fast Fourier transforms of + * power of two sizes of `size` or smaller. `size` must be a + * power of two. + */ + this(size_t size) + { + import core.stdc.stdlib; if (size == 0) + { + this([]); + } + else { + immutable bufferSize = 2*size; + this.pStorage = cast(lookup_t*)malloc(lookup_t.sizeof*bufferSize); + this(pStorage[0..bufferSize]); + } + } + + /**This constructor takes a preallocated buffer for a lookup table. + * The tablemust be twice as big as the desired FFT size. + * + * Unsafe because the `memSpace` buffer will be cast to `immutable`. + */ + this(return scope lookup_t[] memSpace) + in (memSpace.length == 0 || isPowerOf2(memSpace.length/2), + "Can only do FFTs on ranges with a size that is a power of two.") + { + if (memSpace.length == 0) { return; } - assert(isPowerOf2(size), - "Can only do FFTs on ranges with a size that is a power of two."); + immutable size = memSpace.length / 2; - auto table = new lookup_t[][bsf(size) + 1]; + /* Create a lookup table of all negative sine values at a resolution of + * size and all smaller power of two resolutions. This may seem + * inefficient, but having all the lookups be next to each other in + * memory at every level of iteration is a huge win performance-wise. + * + * Lookups are stored in a dense triangular array: + * [ + * nan, // never used + * nan, + * 0, 0, + * 0, -1, 0, 1, + * 0, -0.7, -1, -0.7, 0, 0.7, 1, 0.7, + * ... + * ] + * The index of the `i`th lookup is `2^^i` and the length is also `2^^i`. + */ - table[$ - 1] = memSpace[$ - size..$]; - memSpace = memSpace[0 .. size]; + memSpace[0..2] = float.nan; - auto lastRow = table[$ - 1]; + auto lastRow = memSpace[$ - size .. $]; lastRow[0] = 0; // -sin(0) == 0. foreach (ptrdiff_t i; 1 .. size) { @@ -3647,41 +3770,36 @@ private: } // Fill in all the other rows with strided versions. - foreach (i; 1 .. table.length - 1) + immutable tableSize = bsf(size) + 1; + foreach (i; 1 .. tableSize - 1) { - immutable strideLength = size / (2 ^^ i); - auto strided = Stride!(lookup_t[])(lastRow, strideLength); - table[i] = memSpace[$ - strided.length..$]; - memSpace = memSpace[0..$ - strided.length]; + immutable idx = 1 << i, len = 1 << i; + auto lookup = memSpace[idx .. idx+len]; + immutable strideLength = size / (1 << i); + auto strided = Stride!(lookup_t[])(lastRow, strideLength); size_t copyIndex; foreach (elem; strided) { - table[i][copyIndex++] = elem; + lookup[copyIndex++] = elem; } } - negSinLookup = cast(immutable) table; + this.table = cast(immutable)memSpace; } -public: - /**Create an `Fft` object for computing fast Fourier transforms of - * power of two sizes of `size` or smaller. `size` must be a - * power of two. - */ - this(size_t size) - { - // Allocate all twiddle factor buffers in one contiguous block so that, - // when one is done being used, the next one is next in cache. - auto memSpace = uninitializedArray!(lookup_t[])(2 * size); - this(memSpace); - } + /// + @disable this(ref FFTImpl); - @property size_t size() const + /// + ~this() @nogc { - return (negSinLookup is null) ? 0 : negSinLookup[$ - 1].length; + import core.stdc.stdlib; + free(pStorage); } + @property size_t size() const => table.length/2; + /**Compute the Fourier transform of range using the $(BIGOH N log N) * Cooley-Tukey Algorithm. `range` must be a random-access range with * slicing and a length equal to `size` as provided at the construction of @@ -3700,8 +3818,9 @@ public: */ Complex!F[] fft(F = double, R)(R range) const if (isFloatingPoint!F && isRandomAccessRange!R) + in (range.length <= size, "FFT size mismatch.") { - enforceSize(range); + import std.array : uninitializedArray; Complex!F[] ret; if (range.length == 0) { @@ -3723,10 +3842,9 @@ public: */ void fft(Ret, R)(R range, Ret buf) const if (isRandomAccessRange!Ret && isComplexLike!(ElementType!Ret) && hasSlicing!Ret) + in (range.length <= size, "FFT size mismatch.") + in (buf.length == range.length) { - assert(buf.length == range.length); - enforceSize(range); - if (range.length == 0) { return; @@ -3781,8 +3899,9 @@ public: */ Complex!F[] inverseFft(F = double, R)(R range) const if (isRandomAccessRange!R && isComplexLike!(ElementType!R) && isFloatingPoint!F) + in (range.length <= size, "FFT size mismatch.") { - enforceSize(range); + import std.array : uninitializedArray; Complex!F[] ret; if (range.length == 0) { @@ -3803,9 +3922,8 @@ public: */ void inverseFft(Ret, R)(R range, Ret buf) const if (isRandomAccessRange!Ret && isComplexLike!(ElementType!Ret) && hasSlicing!Ret) + in (range.length <= size, "FFT size mismatch.") { - enforceSize(range); - auto swapped = map!swapRealImag(range); fft(swapped, buf); @@ -4158,7 +4276,7 @@ in (buf.length == 2) buf[0].im = range[0].im + range[1].im; buf[1].re = range[0].re - range[1].re; buf[1].im = range[0].im - range[1].im; - } +} } // Hard-coded base case for FFT of size 4. Doesn't work as well as the size From 9bf6dbb040be4331cfcef49ae2e0cef1125c2d32 Mon Sep 17 00:00:00 2001 From: Georgy Markov Date: Sat, 7 Jun 2025 20:11:05 +0300 Subject: [PATCH 05/31] FFTImpl(size_t size) : throw error on Out of Memory --- std/numeric.d | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/std/numeric.d b/std/numeric.d index 7c09249c7d0..a43833d6b51 100644 --- a/std/numeric.d +++ b/std/numeric.d @@ -3706,6 +3706,7 @@ public: this(size_t size) { import core.stdc.stdlib; + import core.exception : onOutOfMemoryError; if (size == 0) { this([]); @@ -3713,6 +3714,9 @@ public: else { immutable bufferSize = 2*size; this.pStorage = cast(lookup_t*)malloc(lookup_t.sizeof*bufferSize); + if (!this.pStorage) { + onOutOfMemoryError(); + } this(pStorage[0..bufferSize]); } } From bb489547d0665b0687049aefe619608778ef884b Mon Sep 17 00:00:00 2001 From: Georgy Markov Date: Sat, 7 Jun 2025 20:17:24 +0300 Subject: [PATCH 06/31] Convenience functions use `FFTImpl` instead of `Fft` --- std/numeric.d | 37 ++++--------------------------------- 1 file changed, 4 insertions(+), 33 deletions(-) diff --git a/std/numeric.d b/std/numeric.d index a43833d6b51..263ffe89284 100644 --- a/std/numeric.d +++ b/std/numeric.d @@ -3381,19 +3381,6 @@ private: "FFT size mismatch. Expected ", size, ", got ", range.length)); } - // This constructor is used within this module for allocating the - // buffer space elsewhere besides the GC heap. It's definitely **NOT** - // part of the public API and definitely **IS** subject to change. - // - // Also, this is unsafe because the memSpace buffer will be cast - // to immutable. - // - // Public b/c of https://issues.dlang.org/show_bug.cgi?id=4636. - public this(lookup_t[] memSpace) - { - this.impl = FFTImpl(memSpace); - } - public: /**Create an `Fft` object for computing fast Fourier transforms of * power of two sizes of `size` or smaller. `size` must be a @@ -3941,22 +3928,6 @@ public: } } -// This mixin creates an Fft object in the scope it's mixed into such that all -// memory owned by the object is deterministically destroyed at the end of that -// scope. -private enum string MakeLocalFft = q{ - import core.stdc.stdlib; - import core.exception : onOutOfMemoryError; - - auto lookupBuf = (cast(lookup_t*) malloc(range.length * 2 * lookup_t.sizeof)) - [0 .. 2 * range.length]; - if (!lookupBuf.ptr) - onOutOfMemoryError(); - - scope(exit) free(cast(void*) lookupBuf.ptr); - auto fftObj = scoped!Fft(lookupBuf); -}; - /**Convenience functions that create an `Fft` object, run the FFT or inverse * FFT and return the result. Useful for one-off FFTs. * @@ -3967,28 +3938,28 @@ private enum string MakeLocalFft = q{ */ Complex!F[] fft(F = double, R)(R range) { - mixin(MakeLocalFft); + auto fftObj = FFTImpl(range.length); return fftObj.fft!(F, R)(range); } /// ditto void fft(Ret, R)(R range, Ret buf) { - mixin(MakeLocalFft); + auto fftObj = FFTImpl(range.length); return fftObj.fft!(Ret, R)(range, buf); } /// ditto Complex!F[] inverseFft(F = double, R)(R range) { - mixin(MakeLocalFft); + auto fftObj = FFTImpl(range.length); return fftObj.inverseFft!(F, R)(range); } /// ditto void inverseFft(Ret, R)(R range, Ret buf) { - mixin(MakeLocalFft); + auto fftObj = FFTImpl(range.length); return fftObj.inverseFft!(Ret, R)(range, buf); } From e0e18b37b05bce783cd339cb7c7919744e2d3a10 Mon Sep 17 00:00:00 2001 From: Georgy Markov Date: Sat, 7 Jun 2025 21:06:04 +0300 Subject: [PATCH 07/31] Make `FFTImpl` `@nogc nothrow` --- std/numeric.d | 30 ++++++++++++++++++++++++++---- 1 file changed, 26 insertions(+), 4 deletions(-) diff --git a/std/numeric.d b/std/numeric.d index 263ffe89284..5f08c68df25 100644 --- a/std/numeric.d +++ b/std/numeric.d @@ -3690,7 +3690,7 @@ public: * power of two sizes of `size` or smaller. `size` must be a * power of two. */ - this(size_t size) + this(size_t size) @nogc nothrow { import core.stdc.stdlib; import core.exception : onOutOfMemoryError; @@ -3713,7 +3713,7 @@ public: * * Unsafe because the `memSpace` buffer will be cast to `immutable`. */ - this(return scope lookup_t[] memSpace) + this(return scope lookup_t[] memSpace) @nogc nothrow in (memSpace.length == 0 || isPowerOf2(memSpace.length/2), "Can only do FFTs on ranges with a size that is a power of two.") { @@ -3783,13 +3783,13 @@ public: @disable this(ref FFTImpl); /// - ~this() @nogc + ~this() @nogc nothrow { import core.stdc.stdlib; free(pStorage); } - @property size_t size() const => table.length/2; + @property size_t size() const @nogc nothrow => table.length/2; /**Compute the Fourier transform of range using the $(BIGOH N log N) * Cooley-Tukey Algorithm. `range` must be a random-access range with @@ -4053,6 +4053,28 @@ void inverseFft(Ret, R)(R range, Ret buf) assert(inv[].map!"a.im".maxElement < 1e-10); } +// https://github.com/dlang/phobos/issues/10798 +@nogc nothrow @system unittest +{ + import std.algorithm; + import std.range; + static struct C { float re, im; } // User-defined complex + + float[8] arr = [1,2,3,4,5,6,7,8]; + C[8] fft1; + fft(arr[], fft1[]); + + float[8] re = [36.0, -4, -4, -4, -4, -4, -4, -4]; + float[8] im = [0, 9.6568, 4, 1.6568, 0, -1.6568, -4, -9.6568]; + assert(isClose(fft1[].map!"a.re", re[], 1e-4)); + assert(isClose(fft1[].map!"a.im", im[], 1e-4)); + + C[8] inv; + inverseFft(fft1[], inv[]); + assert(isClose(inv[].map!"a.re", arr[], 1e-6)); + assert(inv[].map!"a.im".maxElement < 1e-10); +} + // Swaps the real and imaginary parts of a complex number. This is useful // for inverse FFTs. C swapRealImag(C)(C input) From 7d85642a983dc42df375ef331feb5fafebb6b4ff Mon Sep 17 00:00:00 2001 From: Georgy Markov Date: Sat, 7 Jun 2025 21:22:50 +0300 Subject: [PATCH 08/31] `FFTImpl`: rename `fft`/`inverseFft` to `compute`/`computeInverse` --- std/numeric.d | 32 ++++++++++++++++---------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/std/numeric.d b/std/numeric.d index 5f08c68df25..44c42482253 100644 --- a/std/numeric.d +++ b/std/numeric.d @@ -3413,7 +3413,7 @@ public: if (isFloatingPoint!F && isRandomAccessRange!R) { enforceSize(range); - return impl.fft(range); + return impl.compute(range); } /**Same as the overload, but allows for the results to be stored in a user- @@ -3427,7 +3427,7 @@ public: in (buf.length == range.length) { enforceSize(range); - impl.fft(range, buf); + impl.compute(range, buf); } /** @@ -3446,7 +3446,7 @@ public: if (isRandomAccessRange!R && isComplexLike!(ElementType!R) && isFloatingPoint!F) { enforceSize(range); - return impl.inverseFft(range); + return impl.computeInverse(range); } /** @@ -3458,7 +3458,7 @@ public: if (isRandomAccessRange!Ret && isComplexLike!(ElementType!Ret) && hasSlicing!Ret) { enforceSize(range); - impl.inverseFft(range, buf); + impl.computeInverse(range, buf); } } @@ -3587,7 +3587,7 @@ private: auto oddsImag = OddToImaginary(range); } - fft(oddsImag, buf[0..$ / 2]); + compute(oddsImag, buf[0..$ / 2]); auto evenFft = buf[0..$ / 2]; auto oddFft = buf[$ / 2..$]; immutable halfN = evenFft.length; @@ -3807,7 +3807,7 @@ public: * Conventions: The exponent is negative and the factor is one, * i.e., output[j] := sum[ exp(-2 PI i j k / N) input[k] ]. */ - Complex!F[] fft(F = double, R)(R range) const + Complex!F[] compute(F = double, R)(R range) const if (isFloatingPoint!F && isRandomAccessRange!R) in (range.length <= size, "FFT size mismatch.") { @@ -3821,7 +3821,7 @@ public: // Don't waste time initializing the memory for ret. ret = uninitializedArray!(Complex!F[])(range.length); - fft(range, ret); + compute(range, ret); return ret; } @@ -3831,7 +3831,7 @@ public: * complex-like. This means that they must have a .re and a .im member or * property that can be both read and written and are floating point numbers. */ - void fft(Ret, R)(R range, Ret buf) const + void compute(Ret, R)(R range, Ret buf) const if (isRandomAccessRange!Ret && isComplexLike!(ElementType!Ret) && hasSlicing!Ret) in (range.length <= size, "FFT size mismatch.") in (buf.length == range.length) @@ -3888,7 +3888,7 @@ public: * Conventions: The exponent is positive and the factor is 1/N, i.e., * output[j] := (1 / N) sum[ exp(+2 PI i j k / N) input[k] ]. */ - Complex!F[] inverseFft(F = double, R)(R range) const + Complex!F[] computeInverse(F = double, R)(R range) const if (isRandomAccessRange!R && isComplexLike!(ElementType!R) && isFloatingPoint!F) in (range.length <= size, "FFT size mismatch.") { @@ -3902,7 +3902,7 @@ public: // Don't waste time initializing the memory for ret. ret = uninitializedArray!(Complex!F[])(range.length); - inverseFft(range, ret); + computeInverse(range, ret); return ret; } @@ -3911,12 +3911,12 @@ public: * must be a random access range with slicing, and its elements * must be some complex-like type. */ - void inverseFft(Ret, R)(R range, Ret buf) const + void computeInverse(Ret, R)(R range, Ret buf) const if (isRandomAccessRange!Ret && isComplexLike!(ElementType!Ret) && hasSlicing!Ret) in (range.length <= size, "FFT size mismatch.") { auto swapped = map!swapRealImag(range); - fft(swapped, buf); + compute(swapped, buf); immutable lenNeg1 = 1.0 / buf.length; foreach (ref elem; buf) @@ -3939,28 +3939,28 @@ public: Complex!F[] fft(F = double, R)(R range) { auto fftObj = FFTImpl(range.length); - return fftObj.fft!(F, R)(range); + return fftObj.compute!(F, R)(range); } /// ditto void fft(Ret, R)(R range, Ret buf) { auto fftObj = FFTImpl(range.length); - return fftObj.fft!(Ret, R)(range, buf); + return fftObj.compute!(Ret, R)(range, buf); } /// ditto Complex!F[] inverseFft(F = double, R)(R range) { auto fftObj = FFTImpl(range.length); - return fftObj.inverseFft!(F, R)(range); + return fftObj.computeInverse!(F, R)(range); } /// ditto void inverseFft(Ret, R)(R range, Ret buf) { auto fftObj = FFTImpl(range.length); - return fftObj.inverseFft!(Ret, R)(range, buf); + return fftObj.computeInverse!(Ret, R)(range, buf); } @system unittest From 434c2e0b71689233c8b4fc8af3e35f5881eaa69c Mon Sep 17 00:00:00 2001 From: Georgy Markov Date: Sat, 7 Jun 2025 23:02:47 +0300 Subject: [PATCH 09/31] Add `std.numeric.FFT` Unlike `Fft`, this can be used in `@nogc` code. --- std/numeric.d | 172 ++++++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 160 insertions(+), 12 deletions(-) diff --git a/std/numeric.d b/std/numeric.d index 44c42482253..cb574c4bcb7 100644 --- a/std/numeric.d +++ b/std/numeric.d @@ -3462,6 +3462,157 @@ public: } } +/**A struct for performing fast Fourier transforms of power of two sizes. + * This class encapsulates a large amount of state that is reusable when + * performing multiple FFTs of sizes smaller than or equal to that specified + * in the constructor. This results in substantial speedups when performing + * multiple FFTs with a known maximum size. However, + * a free function API is provided for convenience if you need to perform a + * one-off FFT. + * + * References: + * $(HTTP en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm) + */ +struct FFT +{ +private: + + alias Impl = SafeRefCounted!(FFTImpl, RefCountedAutoInitialize.no); + Impl impl; + +public: + + /**Create an `FFT` object for computing fast Fourier transforms of + * power of two sizes of `size` or smaller. `size` must be a + * power of two. + */ + this(size_t size) @nogc nothrow + { + this.impl = Impl(size); + } + + /**This constructor takes a preallocated buffer for a lookup table. + * The table must be twice as big as the desired FFT size. + * + * Unsafe because the `memSpace` buffer will be cast to `immutable`. + */ + this(return scope lookup_t[] memSpace) @nogc nothrow + in (memSpace.length == 0 || isPowerOf2(memSpace.length/2), + "Can only do FFTs on ranges with a size that is a power of two.") + { + this.impl = Impl(memSpace); + } + + @property size_t size() const @nogc nothrow + { + if (impl.refCountedStore.isInitialized) + { + return impl.borrow!((const ref fft) => fft.size); + } + else + { + return 0; + } + } + + /**Compute the Fourier transform of range using the $(BIGOH N log N) + * Cooley-Tukey Algorithm. `range` must be a random-access range with + * slicing and a length equal to `size` as provided at the construction of + * this object. The contents of range can be either numeric types, + * which will be interpreted as pure real values, or complex types with + * properties or members `.re` and `.im` that can be read. + * + * Note: Pure real FFTs are automatically detected and the relevant + * optimizations are performed. + * + * Returns: An array of complex numbers representing the transformed data in + * the frequency domain. + * + * Conventions: The exponent is negative and the factor is one, + * i.e., output[j] := sum[ exp(-2 PI i j k / N) input[k] ]. + */ + Complex!F[] compute(F = double, R)(R range) const + if (isFloatingPoint!F && isRandomAccessRange!R) + in (range.length <= size, "FFT size mismatch.") + { + if (impl.refCountedStore.isInitialized) + { + return impl.borrow!((const ref fft) => fft.compute(range)); + } + else + { + return FFTImpl(0).compute(range); + } + } + + /**Same as the overload, but allows for the results to be stored in a user- + * provided buffer. The buffer must be of the same length as range, must be + * a random-access range, must have slicing, and must contain elements that are + * complex-like. This means that they must have a .re and a .im member or + * property that can be both read and written and are floating point numbers. + */ + void compute(Ret, R)(R range, Ret buf) const + if (isRandomAccessRange!Ret && isComplexLike!(ElementType!Ret) && hasSlicing!Ret) + in (range.length <= size, "FFT size mismatch.") + in (buf.length == range.length) + { + if (impl.refCountedStore.isInitialized) + { + impl.borrow!((const ref fft) => fft.compute(range, buf)); + } + else + { + return FFTImpl(0).compute(range, buf); + } + } + + /** + * Computes the inverse Fourier transform of a range. The range must be a + * random access range with slicing, have a length equal to the size + * provided at construction of this object, and contain elements that are + * either of type std.complex.Complex or have essentially + * the same compile-time interface. + * + * Returns: The time-domain signal. + * + * Conventions: The exponent is positive and the factor is 1/N, i.e., + * output[j] := (1 / N) sum[ exp(+2 PI i j k / N) input[k] ]. + */ + Complex!F[] computeInverse(F = double, R)(R range) const + if (isRandomAccessRange!R && isComplexLike!(ElementType!R) && isFloatingPoint!F) + in (range.length <= size, "FFT size mismatch.") + { + if (impl.refCountedStore.isInitialized) + { + impl.borrow!((const ref fft) => fft.computeInverse(range)); + } + else + { + return FFTImpl(0).computeInverse(range); + } + } + + /** + * Inverse FFT that allows a user-supplied buffer to be provided. The buffer + * must be a random access range with slicing, and its elements + * must be some complex-like type. + */ + void computeInverse(Ret, R)(R range, Ret buf) const + if (isRandomAccessRange!Ret && isComplexLike!(ElementType!Ret) && hasSlicing!Ret) + in (range.length <= size, "FFT size mismatch.") + { + if (impl.refCountedStore.isInitialized) + { + impl.borrow!((const ref fft) => fft.computeInverse(range, buf)); + } + else + { + return FFTImpl(0).computeInverse(range, buf); + } + } +} + + private struct FFTImpl { import core.bitop : bsf; @@ -3469,7 +3620,7 @@ private struct FFTImpl private: lookup_t* pStorage; - immutable lookup_t[] table; + immutable(lookup_t)[] table; void fftImpl(Ret, R)(Stride!R range, Ret buf) const in @@ -4056,23 +4207,20 @@ void inverseFft(Ret, R)(R range, Ret buf) // https://github.com/dlang/phobos/issues/10798 @nogc nothrow @system unittest { - import std.algorithm; - import std.range; - static struct C { float re, im; } // User-defined complex + static struct C { float re, im; } float[8] arr = [1,2,3,4,5,6,7,8]; C[8] fft1; + C[8] inv; + fft(arr[], fft1[]); + inverseFft(fft1[], inv[]); - float[8] re = [36.0, -4, -4, -4, -4, -4, -4, -4]; - float[8] im = [0, 9.6568, 4, 1.6568, 0, -1.6568, -4, -9.6568]; - assert(isClose(fft1[].map!"a.re", re[], 1e-4)); - assert(isClose(fft1[].map!"a.im", im[], 1e-4)); + lookup_t[16] buff; + auto fft = FFT(buff); - C[8] inv; - inverseFft(fft1[], inv[]); - assert(isClose(inv[].map!"a.re", arr[], 1e-6)); - assert(inv[].map!"a.im".maxElement < 1e-10); + fft.compute(arr[], fft1[]); + fft.computeInverse(fft1[], inv[]); } // Swaps the real and imaginary parts of a complex number. This is useful From 61ba6c5b02a3668423f2cf36769472773ccfb9f3 Mon Sep 17 00:00:00 2001 From: Georgy Markov Date: Sat, 7 Jun 2025 23:12:13 +0300 Subject: [PATCH 10/31] Make `lookup_t` public --- std/numeric.d | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/std/numeric.d b/std/numeric.d index cb574c4bcb7..a3b5f9b05b0 100644 --- a/std/numeric.d +++ b/std/numeric.d @@ -3357,7 +3357,8 @@ if (!isIntegral!T && // though floats seem accurate enough for all practical purposes, since // they pass the "isClose(inverseFft(fft(arr)), arr)" test even for // size 2 ^^ 22. -private alias lookup_t = float; +/// Value type for `FFT` lookup table +alias lookup_t = float; /**A class for performing fast Fourier transforms of power of two sizes. * This class encapsulates a large amount of state that is reusable when From 134e998c8c1be8b2566522f3db2bedd37be48646 Mon Sep 17 00:00:00 2001 From: Georgy Markov Date: Sat, 7 Jun 2025 23:23:21 +0300 Subject: [PATCH 11/31] FFT.lookupTableSize --- std/numeric.d | 26 +++++++++++++++++++++----- 1 file changed, 21 insertions(+), 5 deletions(-) diff --git a/std/numeric.d b/std/numeric.d index a3b5f9b05b0..7238dde0b17 100644 --- a/std/numeric.d +++ b/std/numeric.d @@ -3483,6 +3483,21 @@ private: public: + /// Returns the size of a lookup table required for the given FFT size. + static size_t lookupTableSize(size_t fftSize) @nogc nothrow => fftSize*2; + + /// + unittest { + enum size = 8; + float[size] arr = [1,2,3,4,5,6,7,8]; + Complex!float[size] fft1; + + lookup_t[FFT.lookupTableSize(size)] buff; + auto fft = FFT(buff); + + fft.compute(arr[], fft1[]); + } + /**Create an `FFT` object for computing fast Fourier transforms of * power of two sizes of `size` or smaller. `size` must be a * power of two. @@ -3493,7 +3508,7 @@ public: } /**This constructor takes a preallocated buffer for a lookup table. - * The table must be twice as big as the desired FFT size. + * Use `lookupTableSize` to calculate the required buffer size. * * Unsafe because the `memSpace` buffer will be cast to `immutable`. */ @@ -4210,14 +4225,15 @@ void inverseFft(Ret, R)(R range, Ret buf) { static struct C { float re, im; } - float[8] arr = [1,2,3,4,5,6,7,8]; - C[8] fft1; - C[8] inv; + enum size = 8; + float[size] arr = [1,2,3,4,5,6,7,8]; + C[size] fft1; + C[size] inv; fft(arr[], fft1[]); inverseFft(fft1[], inv[]); - lookup_t[16] buff; + lookup_t[FFT.lookupTableSize(size)] buff; auto fft = FFT(buff); fft.compute(arr[], fft1[]); From c5265c6b8e4f5d7f042aaeb2c9fa8c87109b11b6 Mon Sep 17 00:00:00 2001 From: Georgy Markov Date: Sat, 7 Jun 2025 23:26:25 +0300 Subject: [PATCH 12/31] =?UTF-8?q?Fix=20dlang#10798=20=E2=80=94=20make=20FF?= =?UTF-8?q?T=20fully=20compatible=20with=20`@nogc`?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit From e5546ac2b001d79c87aba25fb1a1dd269d5a1b1f Mon Sep 17 00:00:00 2001 From: Georgy Markov Date: Sat, 7 Jun 2025 23:37:04 +0300 Subject: [PATCH 13/31] Style --- std/numeric.d | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/std/numeric.d b/std/numeric.d index 7238dde0b17..396307b3de3 100644 --- a/std/numeric.d +++ b/std/numeric.d @@ -3487,7 +3487,8 @@ public: static size_t lookupTableSize(size_t fftSize) @nogc nothrow => fftSize*2; /// - unittest { + unittest + { enum size = 8; float[size] arr = [1,2,3,4,5,6,7,8]; Complex!float[size] fft1; @@ -3865,10 +3866,12 @@ public: { this([]); } - else { + else + { immutable bufferSize = 2*size; this.pStorage = cast(lookup_t*)malloc(lookup_t.sizeof*bufferSize); - if (!this.pStorage) { + if (!this.pStorage) + { onOutOfMemoryError(); } this(pStorage[0..bufferSize]); From 2e17de60f413dc2878f46fc40d4ca31d246c50d5 Mon Sep 17 00:00:00 2001 From: Georgy Markov Date: Sat, 7 Jun 2025 23:46:06 +0300 Subject: [PATCH 14/31] Style --- std/numeric.d | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/std/numeric.d b/std/numeric.d index 396307b3de3..e84c8dc1349 100644 --- a/std/numeric.d +++ b/std/numeric.d @@ -3869,7 +3869,7 @@ public: else { immutable bufferSize = 2*size; - this.pStorage = cast(lookup_t*)malloc(lookup_t.sizeof*bufferSize); + this.pStorage = cast(lookup_t*) malloc(lookup_t.sizeof*bufferSize); if (!this.pStorage) { onOutOfMemoryError(); @@ -3946,7 +3946,7 @@ public: } } - this.table = cast(immutable)memSpace; + this.table = cast(immutable) memSpace; } /// From 86a7c59c5047616a88f588fbfa6ea93a0af130d2 Mon Sep 17 00:00:00 2001 From: Georgy Markov Date: Sat, 7 Jun 2025 23:50:16 +0300 Subject: [PATCH 15/31] Works on my machine but CI reports missing import --- std/numeric.d | 1 + 1 file changed, 1 insertion(+) diff --git a/std/numeric.d b/std/numeric.d index e84c8dc1349..3f1b5dc6d9c 100644 --- a/std/numeric.d +++ b/std/numeric.d @@ -3489,6 +3489,7 @@ public: /// unittest { + import std.complex : Complex; enum size = 8; float[size] arr = [1,2,3,4,5,6,7,8]; Complex!float[size] fft1; From 56574a35e86072d9979393c8dd05a5c096c4c10e Mon Sep 17 00:00:00 2001 From: Georgy Markov Date: Sun, 8 Jun 2025 00:06:56 +0300 Subject: [PATCH 16/31] Style --- std/numeric.d | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/std/numeric.d b/std/numeric.d index 3f1b5dc6d9c..da180233082 100644 --- a/std/numeric.d +++ b/std/numeric.d @@ -3487,7 +3487,7 @@ public: static size_t lookupTableSize(size_t fftSize) @nogc nothrow => fftSize*2; /// - unittest + @system unittest { import std.complex : Complex; enum size = 8; @@ -3861,7 +3861,7 @@ public: */ this(size_t size) @nogc nothrow { - import core.stdc.stdlib; + import core.stdc.stdlib : malloc; import core.exception : onOutOfMemoryError; if (size == 0) { @@ -3875,7 +3875,7 @@ public: { onOutOfMemoryError(); } - this(pStorage[0..bufferSize]); + this(pStorage[0 .. bufferSize]); } } @@ -3912,7 +3912,7 @@ public: * The index of the `i`th lookup is `2^^i` and the length is also `2^^i`. */ - memSpace[0..2] = float.nan; + memSpace[0 .. 2] = float.nan; auto lastRow = memSpace[$ - size .. $]; lastRow[0] = 0; // -sin(0) == 0. @@ -3956,7 +3956,7 @@ public: /// ~this() @nogc nothrow { - import core.stdc.stdlib; + import core.stdc.stdlib : free; free(pStorage); } From 53f11700a22bde910b248a10fec6eaacb3aa4d07 Mon Sep 17 00:00:00 2001 From: Georgy Markov Date: Sun, 8 Jun 2025 01:23:11 +0300 Subject: [PATCH 17/31] Removed `FFT`. `FFT` is now an alias to `FFTImpl` --- std/numeric.d | 180 ++++++-------------------------------------------- 1 file changed, 21 insertions(+), 159 deletions(-) diff --git a/std/numeric.d b/std/numeric.d index da180233082..5951cf3b4df 100644 --- a/std/numeric.d +++ b/std/numeric.d @@ -3463,6 +3463,8 @@ public: } } +alias FFT = FFTImpl; + /**A struct for performing fast Fourier transforms of power of two sizes. * This class encapsulates a large amount of state that is reusable when * performing multiple FFTs of sizes smaller than or equal to that specified @@ -3474,164 +3476,7 @@ public: * References: * $(HTTP en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm) */ -struct FFT -{ -private: - - alias Impl = SafeRefCounted!(FFTImpl, RefCountedAutoInitialize.no); - Impl impl; - -public: - - /// Returns the size of a lookup table required for the given FFT size. - static size_t lookupTableSize(size_t fftSize) @nogc nothrow => fftSize*2; - - /// - @system unittest - { - import std.complex : Complex; - enum size = 8; - float[size] arr = [1,2,3,4,5,6,7,8]; - Complex!float[size] fft1; - - lookup_t[FFT.lookupTableSize(size)] buff; - auto fft = FFT(buff); - - fft.compute(arr[], fft1[]); - } - - /**Create an `FFT` object for computing fast Fourier transforms of - * power of two sizes of `size` or smaller. `size` must be a - * power of two. - */ - this(size_t size) @nogc nothrow - { - this.impl = Impl(size); - } - - /**This constructor takes a preallocated buffer for a lookup table. - * Use `lookupTableSize` to calculate the required buffer size. - * - * Unsafe because the `memSpace` buffer will be cast to `immutable`. - */ - this(return scope lookup_t[] memSpace) @nogc nothrow - in (memSpace.length == 0 || isPowerOf2(memSpace.length/2), - "Can only do FFTs on ranges with a size that is a power of two.") - { - this.impl = Impl(memSpace); - } - - @property size_t size() const @nogc nothrow - { - if (impl.refCountedStore.isInitialized) - { - return impl.borrow!((const ref fft) => fft.size); - } - else - { - return 0; - } - } - - /**Compute the Fourier transform of range using the $(BIGOH N log N) - * Cooley-Tukey Algorithm. `range` must be a random-access range with - * slicing and a length equal to `size` as provided at the construction of - * this object. The contents of range can be either numeric types, - * which will be interpreted as pure real values, or complex types with - * properties or members `.re` and `.im` that can be read. - * - * Note: Pure real FFTs are automatically detected and the relevant - * optimizations are performed. - * - * Returns: An array of complex numbers representing the transformed data in - * the frequency domain. - * - * Conventions: The exponent is negative and the factor is one, - * i.e., output[j] := sum[ exp(-2 PI i j k / N) input[k] ]. - */ - Complex!F[] compute(F = double, R)(R range) const - if (isFloatingPoint!F && isRandomAccessRange!R) - in (range.length <= size, "FFT size mismatch.") - { - if (impl.refCountedStore.isInitialized) - { - return impl.borrow!((const ref fft) => fft.compute(range)); - } - else - { - return FFTImpl(0).compute(range); - } - } - - /**Same as the overload, but allows for the results to be stored in a user- - * provided buffer. The buffer must be of the same length as range, must be - * a random-access range, must have slicing, and must contain elements that are - * complex-like. This means that they must have a .re and a .im member or - * property that can be both read and written and are floating point numbers. - */ - void compute(Ret, R)(R range, Ret buf) const - if (isRandomAccessRange!Ret && isComplexLike!(ElementType!Ret) && hasSlicing!Ret) - in (range.length <= size, "FFT size mismatch.") - in (buf.length == range.length) - { - if (impl.refCountedStore.isInitialized) - { - impl.borrow!((const ref fft) => fft.compute(range, buf)); - } - else - { - return FFTImpl(0).compute(range, buf); - } - } - - /** - * Computes the inverse Fourier transform of a range. The range must be a - * random access range with slicing, have a length equal to the size - * provided at construction of this object, and contain elements that are - * either of type std.complex.Complex or have essentially - * the same compile-time interface. - * - * Returns: The time-domain signal. - * - * Conventions: The exponent is positive and the factor is 1/N, i.e., - * output[j] := (1 / N) sum[ exp(+2 PI i j k / N) input[k] ]. - */ - Complex!F[] computeInverse(F = double, R)(R range) const - if (isRandomAccessRange!R && isComplexLike!(ElementType!R) && isFloatingPoint!F) - in (range.length <= size, "FFT size mismatch.") - { - if (impl.refCountedStore.isInitialized) - { - impl.borrow!((const ref fft) => fft.computeInverse(range)); - } - else - { - return FFTImpl(0).computeInverse(range); - } - } - - /** - * Inverse FFT that allows a user-supplied buffer to be provided. The buffer - * must be a random access range with slicing, and its elements - * must be some complex-like type. - */ - void computeInverse(Ret, R)(R range, Ret buf) const - if (isRandomAccessRange!Ret && isComplexLike!(ElementType!Ret) && hasSlicing!Ret) - in (range.length <= size, "FFT size mismatch.") - { - if (impl.refCountedStore.isInitialized) - { - impl.borrow!((const ref fft) => fft.computeInverse(range, buf)); - } - else - { - return FFTImpl(0).computeInverse(range, buf); - } - } -} - - -private struct FFTImpl +struct FFTImpl { import core.bitop : bsf; import std.algorithm.iteration : map; @@ -3855,6 +3700,23 @@ private: public: + /// Returns the size of a lookup table required for the given FFT size. + static size_t lookupTableSize(size_t fftSize) @nogc nothrow => fftSize*2; + + /// + @system unittest + { + import std.complex : Complex; + enum size = 8; + float[size] arr = [1,2,3,4,5,6,7,8]; + Complex!float[size] fft1; + + lookup_t[FFT.lookupTableSize(size)] buff; + auto fft = FFT(buff); + + fft.compute(arr[], fft1[]); + } + /**Create an `FFT` object for computing fast Fourier transforms of * power of two sizes of `size` or smaller. `size` must be a * power of two. @@ -3880,7 +3742,7 @@ public: } /**This constructor takes a preallocated buffer for a lookup table. - * The tablemust be twice as big as the desired FFT size. + * Use `lookupTableSize` to calculate the required buffer size. * * Unsafe because the `memSpace` buffer will be cast to `immutable`. */ From 3f13129e96f7d8fe2eae38ecef5ec83001a66431 Mon Sep 17 00:00:00 2001 From: Georgy Markov Date: Sun, 8 Jun 2025 01:24:10 +0300 Subject: [PATCH 18/31] Rename `FFTImpl` to `FFT` --- std/numeric.d | 18 ++++++++---------- 1 file changed, 8 insertions(+), 10 deletions(-) diff --git a/std/numeric.d b/std/numeric.d index 5951cf3b4df..21f76d1d68b 100644 --- a/std/numeric.d +++ b/std/numeric.d @@ -3374,7 +3374,7 @@ alias lookup_t = float; final class Fft { private: - FFTImpl impl; + FFT impl; void enforceSize(R)(R range) const { import std.conv : text; @@ -3389,7 +3389,7 @@ public: */ this(size_t size) { - this.impl = FFTImpl(size); + this.impl = FFT(size); } @property size_t size() const => impl.size; @@ -3463,8 +3463,6 @@ public: } } -alias FFT = FFTImpl; - /**A struct for performing fast Fourier transforms of power of two sizes. * This class encapsulates a large amount of state that is reusable when * performing multiple FFTs of sizes smaller than or equal to that specified @@ -3476,7 +3474,7 @@ alias FFT = FFTImpl; * References: * $(HTTP en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm) */ -struct FFTImpl +struct FFT { import core.bitop : bsf; import std.algorithm.iteration : map; @@ -3813,7 +3811,7 @@ public: } /// - @disable this(ref FFTImpl); + @disable this(ref FFT); /// ~this() @nogc nothrow @@ -3971,28 +3969,28 @@ public: */ Complex!F[] fft(F = double, R)(R range) { - auto fftObj = FFTImpl(range.length); + auto fftObj = FFT(range.length); return fftObj.compute!(F, R)(range); } /// ditto void fft(Ret, R)(R range, Ret buf) { - auto fftObj = FFTImpl(range.length); + auto fftObj = FFT(range.length); return fftObj.compute!(Ret, R)(range, buf); } /// ditto Complex!F[] inverseFft(F = double, R)(R range) { - auto fftObj = FFTImpl(range.length); + auto fftObj = FFT(range.length); return fftObj.computeInverse!(F, R)(range); } /// ditto void inverseFft(Ret, R)(R range, Ret buf) { - auto fftObj = FFTImpl(range.length); + auto fftObj = FFT(range.length); return fftObj.computeInverse!(Ret, R)(range, buf); } From e267223f64d834614c07fab420ee8326ba740d2a Mon Sep 17 00:00:00 2001 From: Georgy Markov Date: Mon, 9 Jun 2025 21:57:39 +0300 Subject: [PATCH 19/31] Make `FFT` `pure` --- std/numeric.d | 31 +++++++++++++++++++++++-------- 1 file changed, 23 insertions(+), 8 deletions(-) diff --git a/std/numeric.d b/std/numeric.d index 21f76d1d68b..cf0765740f8 100644 --- a/std/numeric.d +++ b/std/numeric.d @@ -3719,9 +3719,9 @@ public: * power of two sizes of `size` or smaller. `size` must be a * power of two. */ - this(size_t size) @nogc nothrow + this(size_t size) @nogc nothrow pure { - import core.stdc.stdlib : malloc; + import core.memory : pureMalloc; import core.exception : onOutOfMemoryError; if (size == 0) { @@ -3730,7 +3730,7 @@ public: else { immutable bufferSize = 2*size; - this.pStorage = cast(lookup_t*) malloc(lookup_t.sizeof*bufferSize); + this.pStorage = cast(lookup_t*) pureMalloc(lookup_t.sizeof*bufferSize); if (!this.pStorage) { onOutOfMemoryError(); @@ -3744,7 +3744,7 @@ public: * * Unsafe because the `memSpace` buffer will be cast to `immutable`. */ - this(return scope lookup_t[] memSpace) @nogc nothrow + this(return scope lookup_t[] memSpace) @nogc nothrow pure in (memSpace.length == 0 || isPowerOf2(memSpace.length/2), "Can only do FFTs on ranges with a size that is a power of two.") { @@ -3814,13 +3814,13 @@ public: @disable this(ref FFT); /// - ~this() @nogc nothrow + ~this() @nogc nothrow pure { - import core.stdc.stdlib : free; - free(pStorage); + import core.memory : pureFree; + pureFree(pStorage); } - @property size_t size() const @nogc nothrow => table.length/2; + @property size_t size() const @nogc nothrow pure => table.length/2; /**Compute the Fourier transform of range using the $(BIGOH N log N) * Cooley-Tukey Algorithm. `range` must be a random-access range with @@ -4104,6 +4104,21 @@ void inverseFft(Ret, R)(R range, Ret buf) fft.computeInverse(fft1[], inv[]); } +@nogc nothrow pure @system unittest +{ + static struct C { float re, im; } + + enum size = 8; + float[size] arr = [1,2,3,4,5,6,7,8]; + C[size] fft1; + C[size] inv; + + immutable fft = FFT(size); + + fft.compute(arr[], fft1[]); + fft.computeInverse(fft1[], inv[]); +} + // Swaps the real and imaginary parts of a complex number. This is useful // for inverse FFTs. C swapRealImag(C)(C input) From 13441e64446d00ef59b385a7f164083b96db69e6 Mon Sep 17 00:00:00 2001 From: Georgy Markov Date: Mon, 9 Jun 2025 22:21:32 +0300 Subject: [PATCH 20/31] Mark unittests `pure` --- std/numeric.d | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/std/numeric.d b/std/numeric.d index cf0765740f8..96cf1a59f6c 100644 --- a/std/numeric.d +++ b/std/numeric.d @@ -3994,7 +3994,7 @@ void inverseFft(Ret, R)(R range, Ret buf) return fftObj.computeInverse!(Ret, R)(range, buf); } -@system unittest +pure @system unittest { import std.algorithm; import std.conv; @@ -4065,7 +4065,7 @@ void inverseFft(Ret, R)(R range, Ret buf) } // https://github.com/dlang/phobos/issues/10796 -@system unittest +pure @system unittest { import std.algorithm; import std.range; @@ -4085,7 +4085,7 @@ void inverseFft(Ret, R)(R range, Ret buf) } // https://github.com/dlang/phobos/issues/10798 -@nogc nothrow @system unittest +@nogc nothrow pure @system unittest { static struct C { float re, im; } From e3820928df5e41861548c1570231c120e550c8b4 Mon Sep 17 00:00:00 2001 From: Georgy Markov Date: Mon, 9 Jun 2025 22:41:28 +0300 Subject: [PATCH 21/31] Move lookup generation to a function --- std/numeric.d | 121 ++++++++++++++++++++++++++------------------------ 1 file changed, 62 insertions(+), 59 deletions(-) diff --git a/std/numeric.d b/std/numeric.d index 96cf1a59f6c..c5445a33747 100644 --- a/std/numeric.d +++ b/std/numeric.d @@ -3483,6 +3483,67 @@ private: lookup_t* pStorage; immutable(lookup_t)[] table; + /* Create a lookup table of all negative sine values at a resolution of + * size and all smaller power of two resolutions. This may seem + * inefficient, but having all the lookups be next to each other in + * memory at every level of iteration is a huge win performance-wise. + * + * Lookups are stored in a dense triangular array: + * [ + * nan, // never used + * nan, + * 0, 0, + * 0, -1, 0, 1, + * 0, -0.7, -1, -0.7, 0, 0.7, 1, 0.7, + * ... + * ] + * The index of the `i`th lookup is `2^^i` and the length is also `2^^i`. + */ + static void fillLookupTable(scope lookup_t[] memSpace) @nogc nothrow pure + { + if (memSpace.length == 0) + { + return; + } + + immutable size = memSpace.length / 2; + + memSpace[0 .. 2] = float.nan; + + auto lastRow = memSpace[$ - size .. $]; + lastRow[0] = 0; // -sin(0) == 0. + foreach (ptrdiff_t i; 1 .. size) + { + // The hard coded cases are for improved accuracy and to prevent + // annoying non-zeroness when stuff should be zero. + + if (i == size / 4) + lastRow[i] = -1; // -sin(pi / 2) == -1. + else if (i == size / 2) + lastRow[i] = 0; // -sin(pi) == 0. + else if (i == size * 3 / 4) + lastRow[i] = 1; // -sin(pi * 3 / 2) == 1 + else + lastRow[i] = -sin(i * 2.0L * PI / size); + } + + // Fill in all the other rows with strided versions. + immutable tableSize = bsf(size) + 1; + foreach (i; 1 .. tableSize - 1) + { + immutable idx = 1 << i, len = 1 << i; + auto lookup = memSpace[idx .. idx+len]; + + immutable strideLength = size / (1 << i); + auto strided = Stride!(lookup_t[])(lastRow, strideLength); + size_t copyIndex; + foreach (elem; strided) + { + lookup[copyIndex++] = elem; + } + } + } + void fftImpl(Ret, R)(Stride!R range, Ret buf) const in { @@ -3748,65 +3809,7 @@ public: in (memSpace.length == 0 || isPowerOf2(memSpace.length/2), "Can only do FFTs on ranges with a size that is a power of two.") { - if (memSpace.length == 0) - { - return; - } - - immutable size = memSpace.length / 2; - - /* Create a lookup table of all negative sine values at a resolution of - * size and all smaller power of two resolutions. This may seem - * inefficient, but having all the lookups be next to each other in - * memory at every level of iteration is a huge win performance-wise. - * - * Lookups are stored in a dense triangular array: - * [ - * nan, // never used - * nan, - * 0, 0, - * 0, -1, 0, 1, - * 0, -0.7, -1, -0.7, 0, 0.7, 1, 0.7, - * ... - * ] - * The index of the `i`th lookup is `2^^i` and the length is also `2^^i`. - */ - - memSpace[0 .. 2] = float.nan; - - auto lastRow = memSpace[$ - size .. $]; - lastRow[0] = 0; // -sin(0) == 0. - foreach (ptrdiff_t i; 1 .. size) - { - // The hard coded cases are for improved accuracy and to prevent - // annoying non-zeroness when stuff should be zero. - - if (i == size / 4) - lastRow[i] = -1; // -sin(pi / 2) == -1. - else if (i == size / 2) - lastRow[i] = 0; // -sin(pi) == 0. - else if (i == size * 3 / 4) - lastRow[i] = 1; // -sin(pi * 3 / 2) == 1 - else - lastRow[i] = -sin(i * 2.0L * PI / size); - } - - // Fill in all the other rows with strided versions. - immutable tableSize = bsf(size) + 1; - foreach (i; 1 .. tableSize - 1) - { - immutable idx = 1 << i, len = 1 << i; - auto lookup = memSpace[idx .. idx+len]; - - immutable strideLength = size / (1 << i); - auto strided = Stride!(lookup_t[])(lastRow, strideLength); - size_t copyIndex; - foreach (elem; strided) - { - lookup[copyIndex++] = elem; - } - } - + fillLookupTable(memSpace); this.table = cast(immutable) memSpace; } From 107e115d9009c3fac4cc1ff2f378b40385cec0c4 Mon Sep 17 00:00:00 2001 From: Georgy Markov Date: Mon, 9 Jun 2025 23:04:33 +0300 Subject: [PATCH 22/31] Unsafe constructor should take `void[]`, not `lookup_t[]` The type and the size of a lookup table is an implementation detail. --- std/numeric.d | 31 +++++++++++++++++++++---------- 1 file changed, 21 insertions(+), 10 deletions(-) diff --git a/std/numeric.d b/std/numeric.d index c5445a33747..fbf83f31a53 100644 --- a/std/numeric.d +++ b/std/numeric.d @@ -3544,6 +3544,14 @@ private: } } + this(return scope lookup_t[] memSpace) @nogc nothrow pure + in (memSpace.length == 0 || isPowerOf2(memSpace.length/2), + "Can only do FFTs on ranges with a size that is a power of two.") + { + fillLookupTable(memSpace); + this.table = cast(immutable) memSpace; + } + void fftImpl(Ret, R)(Stride!R range, Ret buf) const in { @@ -3759,8 +3767,10 @@ private: public: - /// Returns the size of a lookup table required for the given FFT size. - static size_t lookupTableSize(size_t fftSize) @nogc nothrow => fftSize*2; + /// Returns the size of a buffer required for the given FFT size. + static size_t requiredBufferSize(size_t fftSize) @nogc nothrow pure { + return lookup_t.sizeof*fftSize*2; + } /// @system unittest @@ -3770,8 +3780,8 @@ public: float[size] arr = [1,2,3,4,5,6,7,8]; Complex!float[size] fft1; - lookup_t[FFT.lookupTableSize(size)] buff; - auto fft = FFT(buff); + ubyte[FFT.requiredBufferSize(size)] buff; + auto fft = FFT(size, buff); fft.compute(arr[], fft1[]); } @@ -3805,12 +3815,13 @@ public: * * Unsafe because the `memSpace` buffer will be cast to `immutable`. */ - this(return scope lookup_t[] memSpace) @nogc nothrow pure - in (memSpace.length == 0 || isPowerOf2(memSpace.length/2), + this(size_t size, return scope void[] buffer) @nogc nothrow pure @system + in (size == 0 || isPowerOf2(size), "Can only do FFTs on ranges with a size that is a power of two.") + in (buffer.length == requiredBufferSize(size)) { - fillLookupTable(memSpace); - this.table = cast(immutable) memSpace; + auto memSpace = cast(lookup_t[]) buffer; + this(memSpace); } /// @@ -4100,8 +4111,8 @@ pure @system unittest fft(arr[], fft1[]); inverseFft(fft1[], inv[]); - lookup_t[FFT.lookupTableSize(size)] buff; - auto fft = FFT(buff); + ubyte[FFT.requiredBufferSize(size)] buff; + auto fft = FFT(size, buff); fft.compute(arr[], fft1[]); fft.computeInverse(fft1[], inv[]); From 3bddd34b2fffc56f0ab9f5a24d308ed15c5b4fe1 Mon Sep 17 00:00:00 2001 From: Georgy Markov Date: Mon, 9 Jun 2025 23:05:12 +0300 Subject: [PATCH 23/31] Make `lookup_t` private, as it used to be --- std/numeric.d | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/std/numeric.d b/std/numeric.d index fbf83f31a53..6f69f7c8100 100644 --- a/std/numeric.d +++ b/std/numeric.d @@ -3357,8 +3357,7 @@ if (!isIntegral!T && // though floats seem accurate enough for all practical purposes, since // they pass the "isClose(inverseFft(fft(arr)), arr)" test even for // size 2 ^^ 22. -/// Value type for `FFT` lookup table -alias lookup_t = float; +private alias lookup_t = float; /**A class for performing fast Fourier transforms of power of two sizes. * This class encapsulates a large amount of state that is reusable when From 7dc1d0bdd691fd3d384e01c037272ee12d02514d Mon Sep 17 00:00:00 2001 From: Georgy Markov Date: Mon, 9 Jun 2025 23:19:25 +0300 Subject: [PATCH 24/31] Make `FFT.this(size_t)` `@safe` --- std/numeric.d | 34 +++++++++++++++++++--------------- 1 file changed, 19 insertions(+), 15 deletions(-) diff --git a/std/numeric.d b/std/numeric.d index 6f69f7c8100..03ad951b3fe 100644 --- a/std/numeric.d +++ b/std/numeric.d @@ -3498,7 +3498,7 @@ private: * ] * The index of the `i`th lookup is `2^^i` and the length is also `2^^i`. */ - static void fillLookupTable(scope lookup_t[] memSpace) @nogc nothrow pure + static void fillLookupTable(scope lookup_t[] memSpace) @nogc nothrow pure @safe { if (memSpace.length == 0) { @@ -3543,7 +3543,7 @@ private: } } - this(return scope lookup_t[] memSpace) @nogc nothrow pure + this(return scope lookup_t[] memSpace) @nogc nothrow pure @trusted in (memSpace.length == 0 || isPowerOf2(memSpace.length/2), "Can only do FFTs on ranges with a size that is a power of two.") { @@ -3767,7 +3767,7 @@ private: public: /// Returns the size of a buffer required for the given FFT size. - static size_t requiredBufferSize(size_t fftSize) @nogc nothrow pure { + static size_t requiredBufferSize(size_t fftSize) @nogc nothrow pure @safe { return lookup_t.sizeof*fftSize*2; } @@ -3789,7 +3789,9 @@ public: * power of two sizes of `size` or smaller. `size` must be a * power of two. */ - this(size_t size) @nogc nothrow pure + this(size_t size) @nogc nothrow pure @safe + in (size == 0 || isPowerOf2(size), + "Can only do FFTs on ranges with a size that is a power of two.") { import core.memory : pureMalloc; import core.exception : onOutOfMemoryError; @@ -3799,13 +3801,12 @@ public: } else { - immutable bufferSize = 2*size; - this.pStorage = cast(lookup_t*) pureMalloc(lookup_t.sizeof*bufferSize); - if (!this.pStorage) - { - onOutOfMemoryError(); - } - this(pStorage[0 .. bufferSize]); + immutable nbytes = requiredBufferSize(size); + this.pStorage = (() @trusted => cast(lookup_t*) pureMalloc(nbytes))(); + if (!this.pStorage) { onOutOfMemoryError(); } + immutable nitems = nbytes/lookup_t.sizeof; + lookup_t[] memSpace = (() @trusted => this.pStorage[0 .. nitems])(); + this(memSpace); } } @@ -3827,13 +3828,13 @@ public: @disable this(ref FFT); /// - ~this() @nogc nothrow pure + ~this() @nogc nothrow pure @safe { import core.memory : pureFree; - pureFree(pStorage); + (() @trusted => pureFree(pStorage))(); } - @property size_t size() const @nogc nothrow pure => table.length/2; + @property size_t size() const @nogc nothrow pure @safe => table.length/2; /**Compute the Fourier transform of range using the $(BIGOH N log N) * Cooley-Tukey Algorithm. `range` must be a random-access range with @@ -4117,7 +4118,7 @@ pure @system unittest fft.computeInverse(fft1[], inv[]); } -@nogc nothrow pure @system unittest +@nogc nothrow pure @safe unittest { static struct C { float re, im; } @@ -4126,6 +4127,9 @@ pure @system unittest C[size] fft1; C[size] inv; + fft(arr[], fft1[]); + inverseFft(fft1[], inv[]); + immutable fft = FFT(size); fft.compute(arr[], fft1[]); From 290b907a23d39173f061a5f5df75c764917a8cc1 Mon Sep 17 00:00:00 2001 From: Georgy Markov Date: Mon, 9 Jun 2025 23:34:37 +0300 Subject: [PATCH 25/31] `requiredBufferSize` : check for overflow --- std/numeric.d | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/std/numeric.d b/std/numeric.d index 03ad951b3fe..0493f3068f6 100644 --- a/std/numeric.d +++ b/std/numeric.d @@ -3767,8 +3767,14 @@ private: public: /// Returns the size of a buffer required for the given FFT size. - static size_t requiredBufferSize(size_t fftSize) @nogc nothrow pure @safe { - return lookup_t.sizeof*fftSize*2; + static size_t requiredBufferSize(size_t fftSize) @nogc nothrow pure @safe + { + import core.checkedint : mulu; + immutable nitems = fftSize*2; + bool overflow; + immutable nbytes = mulu(nitems, lookup_t.sizeof, overflow); + if (overflow) assert(false, "multiplication overflowed"); + return nbytes; } /// @@ -3790,7 +3796,7 @@ public: * power of two. */ this(size_t size) @nogc nothrow pure @safe - in (size == 0 || isPowerOf2(size), + in (size == 0 || isPowerOf2(size), "Can only do FFTs on ranges with a size that is a power of two.") { import core.memory : pureMalloc; From 69fa61117937f9617e5cb23480129fbe72125ded Mon Sep 17 00:00:00 2001 From: Georgy Markov Date: Sat, 14 Jun 2025 21:58:46 +0300 Subject: [PATCH 26/31] Fix documentation --- std/numeric.d | 9 +-------- 1 file changed, 1 insertion(+), 8 deletions(-) diff --git a/std/numeric.d b/std/numeric.d index 0493f3068f6..f1d3f0225d7 100644 --- a/std/numeric.d +++ b/std/numeric.d @@ -3463,7 +3463,7 @@ public: } /**A struct for performing fast Fourier transforms of power of two sizes. - * This class encapsulates a large amount of state that is reusable when + * This struct encapsulates a large amount of state that is reusable when * performing multiple FFTs of sizes smaller than or equal to that specified * in the constructor. This results in substantial speedups when performing * multiple FFTs with a known maximum size. However, @@ -3818,8 +3818,6 @@ public: /**This constructor takes a preallocated buffer for a lookup table. * Use `lookupTableSize` to calculate the required buffer size. - * - * Unsafe because the `memSpace` buffer will be cast to `immutable`. */ this(size_t size, return scope void[] buffer) @nogc nothrow pure @system in (size == 0 || isPowerOf2(size), @@ -3981,11 +3979,6 @@ public: /**Convenience functions that create an `Fft` object, run the FFT or inverse * FFT and return the result. Useful for one-off FFTs. - * - * Note: In addition to convenience, these functions are slightly more - * efficient than manually creating an Fft object for a single use, - * as the Fft object is deterministically destroyed before these - * functions return. */ Complex!F[] fft(F = double, R)(R range) { From 08edea909bc9fb4c9f1b89c2d513ef3ba0183429 Mon Sep 17 00:00:00 2001 From: Georgy Markov Date: Sat, 14 Jun 2025 22:00:33 +0300 Subject: [PATCH 27/31] Fix missing identation --- std/numeric.d | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/std/numeric.d b/std/numeric.d index f1d3f0225d7..27a033b1c35 100644 --- a/std/numeric.d +++ b/std/numeric.d @@ -4333,7 +4333,7 @@ in (buf.length == 2) buf[0].im = range[0].im + range[1].im; buf[1].re = range[0].re - range[1].re; buf[1].im = range[0].im - range[1].im; -} + } } // Hard-coded base case for FFT of size 4. Doesn't work as well as the size From ed3ed857a0e390864f23d260c8c3dd5d66479362 Mon Sep 17 00:00:00 2001 From: Georgy Markov Date: Sat, 14 Jun 2025 22:13:30 +0300 Subject: [PATCH 28/31] Make a pointer to immutable data also `immutable` --- std/numeric.d | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/std/numeric.d b/std/numeric.d index 27a033b1c35..89c97060695 100644 --- a/std/numeric.d +++ b/std/numeric.d @@ -3479,7 +3479,7 @@ struct FFT import std.algorithm.iteration : map; private: - lookup_t* pStorage; + immutable(lookup_t)* pStorage; immutable(lookup_t)[] table; /* Create a lookup table of all negative sine values at a resolution of @@ -3808,11 +3808,12 @@ public: else { immutable nbytes = requiredBufferSize(size); - this.pStorage = (() @trusted => cast(lookup_t*) pureMalloc(nbytes))(); - if (!this.pStorage) { onOutOfMemoryError(); } + auto p = (() @trusted => cast(lookup_t*) pureMalloc(nbytes))(); + if (!p) { onOutOfMemoryError(); } immutable nitems = nbytes/lookup_t.sizeof; - lookup_t[] memSpace = (() @trusted => this.pStorage[0 .. nitems])(); + lookup_t[] memSpace = (() @trusted => p[0 .. nitems])(); this(memSpace); + this.pStorage = (() @trusted => cast(immutable) p)(); } } @@ -3835,7 +3836,7 @@ public: ~this() @nogc nothrow pure @safe { import core.memory : pureFree; - (() @trusted => pureFree(pStorage))(); + (() @trusted => pureFree(cast(lookup_t*) pStorage))(); } @property size_t size() const @nogc nothrow pure @safe => table.length/2; From cf48c8c1c182db9191343fc6f693e6e224020eb7 Mon Sep 17 00:00:00 2001 From: Georgy Markov Date: Sat, 14 Jun 2025 22:41:51 +0300 Subject: [PATCH 29/31] Sunset `Fft` --- std/numeric.d | 55 +++------------------------------------------------ 1 file changed, 3 insertions(+), 52 deletions(-) diff --git a/std/numeric.d b/std/numeric.d index 89c97060695..14ca69aa6ce 100644 --- a/std/numeric.d +++ b/std/numeric.d @@ -3359,16 +3359,10 @@ if (!isIntegral!T && // size 2 ^^ 22. private alias lookup_t = float; -/**A class for performing fast Fourier transforms of power of two sizes. - * This class encapsulates a large amount of state that is reusable when - * performing multiple FFTs of sizes smaller than or equal to that specified - * in the constructor. This results in substantial speedups when performing - * multiple FFTs with a known maximum size. However, - * a free function API is provided for convenience if you need to perform a - * one-off FFT. +/**The old version of $(LREF FFT) that uses GC. * - * References: - * $(HTTP en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm) + * It is recommended to use `FFT` in new code. `FFT` allocates + * and deterministically frees non-GC memory, or uses a preallocated buffer. */ final class Fft { @@ -3382,10 +3376,6 @@ private: } public: - /**Create an `Fft` object for computing fast Fourier transforms of - * power of two sizes of `size` or smaller. `size` must be a - * power of two. - */ this(size_t size) { this.impl = FFT(size); @@ -3393,22 +3383,6 @@ public: @property size_t size() const => impl.size; - /**Compute the Fourier transform of range using the $(BIGOH N log N) - * Cooley-Tukey Algorithm. `range` must be a random-access range with - * slicing and a length equal to `size` as provided at the construction of - * this object. The contents of range can be either numeric types, - * which will be interpreted as pure real values, or complex types with - * properties or members `.re` and `.im` that can be read. - * - * Note: Pure real FFTs are automatically detected and the relevant - * optimizations are performed. - * - * Returns: An array of complex numbers representing the transformed data in - * the frequency domain. - * - * Conventions: The exponent is negative and the factor is one, - * i.e., output[j] := sum[ exp(-2 PI i j k / N) input[k] ]. - */ Complex!F[] fft(F = double, R)(R range) const if (isFloatingPoint!F && isRandomAccessRange!R) { @@ -3416,12 +3390,6 @@ public: return impl.compute(range); } - /**Same as the overload, but allows for the results to be stored in a user- - * provided buffer. The buffer must be of the same length as range, must be - * a random-access range, must have slicing, and must contain elements that are - * complex-like. This means that they must have a .re and a .im member or - * property that can be both read and written and are floating point numbers. - */ void fft(Ret, R)(R range, Ret buf) const if (isRandomAccessRange!Ret && isComplexLike!(ElementType!Ret) && hasSlicing!Ret) in (buf.length == range.length) @@ -3430,18 +3398,6 @@ public: impl.compute(range, buf); } - /** - * Computes the inverse Fourier transform of a range. The range must be a - * random access range with slicing, have a length equal to the size - * provided at construction of this object, and contain elements that are - * either of type std.complex.Complex or have essentially - * the same compile-time interface. - * - * Returns: The time-domain signal. - * - * Conventions: The exponent is positive and the factor is 1/N, i.e., - * output[j] := (1 / N) sum[ exp(+2 PI i j k / N) input[k] ]. - */ Complex!F[] inverseFft(F = double, R)(R range) const if (isRandomAccessRange!R && isComplexLike!(ElementType!R) && isFloatingPoint!F) { @@ -3449,11 +3405,6 @@ public: return impl.computeInverse(range); } - /** - * Inverse FFT that allows a user-supplied buffer to be provided. The buffer - * must be a random access range with slicing, and its elements - * must be some complex-like type. - */ void inverseFft(Ret, R)(R range, Ret buf) const if (isRandomAccessRange!Ret && isComplexLike!(ElementType!Ret) && hasSlicing!Ret) { From 8edebe51d0743ef806483fd63763e1cd880db685 Mon Sep 17 00:00:00 2001 From: Georgy Markov Date: Sat, 14 Jun 2025 22:57:25 +0300 Subject: [PATCH 30/31] Mark old unittests `@safe` --- std/numeric.d | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/std/numeric.d b/std/numeric.d index 14ca69aa6ce..bdfdb9a0837 100644 --- a/std/numeric.d +++ b/std/numeric.d @@ -3959,7 +3959,7 @@ void inverseFft(Ret, R)(R range, Ret buf) return fftObj.computeInverse!(Ret, R)(range, buf); } -pure @system unittest +pure @safe unittest { import std.algorithm; import std.conv; @@ -4030,7 +4030,7 @@ pure @system unittest } // https://github.com/dlang/phobos/issues/10796 -pure @system unittest +pure @safe unittest { import std.algorithm; import std.range; From f941d3d1eabffb742a2e1b18ae2e932cf1aac737 Mon Sep 17 00:00:00 2001 From: Georgy Markov Date: Sun, 15 Jun 2025 14:26:46 +0300 Subject: [PATCH 31/31] Add changelog entry --- changelog/add_nogc_fft.dd | 45 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 45 insertions(+) create mode 100644 changelog/add_nogc_fft.dd diff --git a/changelog/add_nogc_fft.dd b/changelog/add_nogc_fft.dd new file mode 100644 index 00000000000..bd214f729da --- /dev/null +++ b/changelog/add_nogc_fft.dd @@ -0,0 +1,45 @@ +Add `FFT` that doesn’t use GC. + +`FFT` is an FFT solver compatible with `@nogc`, `nothrow`, `@safe`, and `pure`. +```D +@nogc nothrow pure @safe void fun() +{ + import std.complex; + + enum size = 8; + float[size] arr = [1,2,3,4,5,6,7,8]; + Complex!float[size] fft; + + immutable fftSolver = FFT(size); + + fftSolver.compute(arr[], fft[]); +} +``` + +Convenience functions now use `FFT`. In effect, `fft(R, Ret)` and +`inverseFft(R, Ret)` also can be used in `@nogc`, `nothrow`, `pure`, +and `@safe` code. + +`FFT` can also avoid dynamic memory allocation by using a preallocated buffer. +However, this is not `@safe`. The buffer must not be modified or freed within +the lifetime of the `FFT` instance. + +```D +@nogc nothrow pure @system void unsafeFun() +{ + import std.complex; + + enum size = 8; + float[size] arr = [1,2,3,4,5,6,7,8]; + Complex!float[size] fft; + + enum bufferSize = FFT.requiredBufferSize(size); + ubyte[bufferSize] buff; + immutable fftSolver = FFT(size, buff); + + fftSolver.compute(arr[], fft[]); +} +``` + +`Fft` class is still available for backwards compatibility. It is not recommended +for new code.