On Fri, Nov 18, 2011 at 11:42 AM, Robert Kern wrote:

On Fri, Nov 18, 2011 at 11:19, Dag Sverre Seljebotn wrote:

...

I've been in touch with Martin Reinecke, author of the libpsht code for spherical harmonic transforms, about licensing issues.

libpsht itself will remain under the GPL, but he is likely to release his C port of FFTPACK under BSD in the near future, as it is based on the public domain FFTPACK.

I'm grateful for this change for my own purposes (allows releasing my own competing SHT library under the BSD) -- but it could perhaps be useful for NumPy or SciPy as well, depending on how complete the port is? E.g., perhaps make numpy.fft more complete (is the numpy.fft/scipy.fftpack split simply because of the Fortran dependency?).

It used to be the case that scipy.fftpack allowed one to build against multiple different, usually faster, FFT libraries like FFTW. I think we have backed away from that since the cost of maintaining the build configuration for all of those different backends was so high. It's worth noting that numpy.fft is already using a C translation of FFTPACK. I'm not sure what the differences are between this translation and Martin's.

Having a Bluestein transformation alone would be worthwhile, as it would avoid the N^2 penalty for prime sizes. I am wondering about precision issues, though (when I tried implementing bluestein transforms on top of fftpack, it gave very bad results numerically-wise). A comparison with fftw would be good here. regards, David

On Fri, Nov 18, 2011 at 5:18 AM, Dag Sverre Seljebotn < d.s.seljebotn@astro.uio.no> wrote:

On 11/18/2011 12:58 PM, David Cournapeau wrote:

...

On Fri, Nov 18, 2011 at 11:42 AM, Robert Kern wrote:

...

On Fri, Nov 18, 2011 at 11:19, Dag Sverre Seljebotn wrote:

...

I've been in touch with Martin Reinecke, author of the libpsht code for spherical harmonic transforms, about licensing issues.

libpsht itself will remain under the GPL, but he is likely to release his C port of FFTPACK under BSD in the near future, as it is based on the public domain FFTPACK.

I'm grateful for this change for my own purposes (allows releasing my own competing SHT library under the BSD) -- but it could perhaps be useful for NumPy or SciPy as well, depending on how complete the port is? E.g., perhaps make numpy.fft more complete (is the numpy.fft/scipy.fftpack split simply because of the Fortran dependency?).

It used to be the case that scipy.fftpack allowed one to build against multiple different, usually faster, FFT libraries like FFTW. I think we have backed away from that since the cost of maintaining the build configuration for all of those different backends was so high. It's worth noting that numpy.fft is already using a C translation of FFTPACK. I'm not sure what the differences are between this translation and Martin's.

Here's some more info forwarded from Martin:

""" - only FFTs are supported (no DCTs/DSTs) - only double precision is supported (extension to single precision might not be much work, though) - both complex and real FFTs are supported - real FFTs allow various storage schemes for the (half)complex frequency domain data (classic FFTPACK scheme, FFTW or halfcomplex scheme, uncompressed complex storage) - precision of transforms involving large prime factors should be slightly better than with original FFTPACK - Bluestein's algorithm is automatically selected if considered profitable - small accuracy self-testing code is provided.

Fairly complete interface documentation is available in Doxygen format. I'll prepare a source package later in the afternoon and send it around.

Best regards, Martin """

...

Having a Bluestein transformation alone would be worthwhile, as it would avoid the N^2 penalty for prime sizes.

I am wondering about precision issues, though (when I tried implementing bluestein transforms on top of fftpack, it gave very bad results numerically-wise). A comparison with fftw would be good here.

Well, there's an indirect comparison: My SHT code currently uses FFTW3, and it manages to agree with Reinecke's SHT code to a precision of better than ~1e-12 for full SHTs. That includes several other sources of errors though.

(That's an average over several different-sized FFTs, of which half has n=8192 and the other half all have different size, starting from 4 and increasing up to 8192 in steps of 4 -- meaning prime factors on the order of 1000).

I agree, a more direct comparison with FFTW would be good.

In more detail from the README:

I replaced the iterative sine and cosine calculations in radfg() and radbg() 17 by an exact calculation, which slightly improves the transform accuracy for 18 real FFTs with lengths containing large prime factors.

The current numpy FFT code could use a review in any case, IIRC, there were some possible integer related issues in the arguments. Having clean documented C code would probably be an improvement and also help with maintenance. I suspect there were iterative computations of sin/cos that were spoiling the precision of the single precision versions (scipy), so fixing those might also be useful. Chuck

On Fri, Nov 18, 2011 at 9:22 AM, Dag Sverre Seljebotn < d.s.seljebotn@astro.uio.no> wrote:

...

On 11/18/2011 12:58 PM, David Cournapeau wrote:

...

On Fri, Nov 18, 2011 at 11:42 AM, Robert Kern wrote:

...

On Fri, Nov 18, 2011 at 11:19, Dag Sverre Seljebotn wrote:

...

I've been in touch with Martin Reinecke, author of the libpsht code for spherical harmonic transforms, about licensing issues.

libpsht itself will remain under the GPL, but he is likely to release his C port of FFTPACK under BSD in the near future, as it is based on the public domain FFTPACK.

I'm grateful for this change for my own purposes (allows releasing my own competing SHT library under the BSD) -- but it could perhaps be useful for NumPy or SciPy as well, depending on how complete the port is? E.g., perhaps make numpy.fft more complete (is the numpy.fft/scipy.fftpack split simply because of the Fortran dependency?).

It used to be the case that scipy.fftpack allowed one to build against multiple different, usually faster, FFT libraries like FFTW. I think we have backed away from that since the cost of maintaining the build configuration for all of those different backends was so high. It's worth noting that numpy.fft is already using a C translation of FFTPACK. I'm not sure what the differences are between this translation and Martin's.

Here's some more info forwarded from Martin:

""" - only FFTs are supported (no DCTs/DSTs) - only double precision is supported (extension to single precision might not be much work, though) - both complex and real FFTs are supported - real FFTs allow various storage schemes for the (half)complex frequency domain data (classic FFTPACK scheme, FFTW or halfcomplex scheme, uncompressed complex storage) - precision of transforms involving large prime factors should be slightly better than with original FFTPACK - Bluestein's algorithm is automatically selected if considered

On 11/18/2011 01:18 PM, Dag Sverre Seljebotn wrote: profitable

...

- small accuracy self-testing code is provided.

Fairly complete interface documentation is available in Doxygen format. I'll prepare a source package later in the afternoon and send it around.

OK, Martin sent me a BSD-ed version of his libfftpack, and I stuck it on github:

https://github.com/dagss/libfftpack

Note that ls_fft.h contains Martin's API for it (plan construction and execution etc.)

#define CC(a,b,c) cc[(a)+ido*((b)+cdim*(c))] Hmm... ;) All in all, the code looks pretty readable, although the coding style needs fixups and some of the macros using complex types are never used and should be removed. The size_t type also needs to be replaced by npy_intp for numpy. I think it is worth taking a closer look at integrating the code. Chuck