ejmahler / RustFFT

RustFFT is a high-performance FFT library written in pure Rust.
Apache License 2.0
692 stars 48 forks source link

support for higher precision types #62

Closed roosephu closed 3 years ago

roosephu commented 3 years ago

Is it possible to support types with higher precision than f64? Converting from f64 to custom type might lose precision, as done in https://github.com/ejmahler/RustFFT/blob/master/src/twiddles.rs#L13.

I'm thinking only converting from usize to FftNum and using num::traits::{Float, FloatConsts} to compute pi and do arithmetic operations (e.g., sqrt, cos, sin).

ejmahler commented 3 years ago

That's an interesting thought. The computation of twiddle factors is the only thing that constrains precision like that, so there isn't much that we'd have to update.

Unfortunately I think we're limited by specialization. I'm very open to supporting this, but it'd have to be done in a way that compiles on stable.

roosephu commented 3 years ago

I made a pull request trying to resolve the issue. The method doesn't rely on specialization but on FftNum implementing num_traits::{Float, FloatConst} so probably it compiles on stable. Could you please take a look?

ejmahler commented 3 years ago

I'm sorry to say that after looking into this, I don't think RustFFT can support this until Rust supports specialization on stable.

I implemented the two fixes to precision problems in rader's algorithm and bluestein's algorithm, but we're fundamentally stuck with f64 precision for twiddles. We can't add trait bounds to FftNum, and specialization will require specialization-aware changes to the definition of FftNum, neither of which are on the table.

I'm closing this for now, but if you have any bursts of inspiration on how to solve this problem, please bring it up.

Gpotier commented 3 years ago

Would it be possible to replace the Copy requirement on rustfft::FftNum by Clone so rug_maths::Float can be used, enabling arbitrary precision FFT ?

ejmahler commented 3 years ago

Unfortunately not, for pretty much the same reason: breaking changes.

Imagine a user of RustFFT made a function generic over FftNum:

fn compute_fft<T: FftNum>(data: &mut [T]) {
   let local = data[0]; // copy happens here
}

In this function, right now they can rely on T being able to Copy. If we removed the copy bound, suddenly their code would stop compiling.

I don't know of anyone who does this, but i don't know for sure that nobody does, and I'd rather not take the chance.

Gpotier commented 3 years ago

I get it, Thanks.

HEnquist commented 3 years ago

Would it be possible to get Copy added to rug_maths?

Gpotier commented 3 years ago

It does not seems to be possible. From the author of rug on rug::Float :

They are not Copy for the same reason Vec and String are not Copy. They have an allocation which is freed when they are dropped, and their implementations can modify the allocated memory. By contrast, the lower level gmp_mpfr_sys::mpfr::mpfr_t is Copy, but that is for the same reason a raw pointer is Copy. There are no destructors for mpfr_t, the mpfr::clear function would need to be called explicitly. Also, writing to the allocated data would require unsafe code where the caller has to take care of aliasing and such, things that Float hides.

So rug_maths::Float which is only a simple wrapper over rug::Float cannot be copy.

HEnquist commented 3 years ago

I took a look at rug_maths, and I think it's not suitable anyway. I looked quickly so might be missing something. But I don't see the things needed to calculate twiddle factors. We need pi, and sin(). If we can't calculate higher precision twiddle factors, then the result won't be much better than a normal f64.

Gpotier commented 3 years ago

Yes you are right rug_maths is missing lots of things. It is probably better to go with rug directly which provides sin() and gmp-mpfr-sys (dependency of rug) which provides binding to mpfr_const_pi() to compute pi in arbitrary precision. I think I will give it a try that way.

overdrivenpotato commented 3 years ago

@ejmahler Is it possible for the Clone bound to be put behind a feature flag? This would avoid breaking the API for users unless e.g. a clone feature is enabled (it can be off by default).