diff options
Diffstat (limited to 'numpy')
| -rw-r--r-- | numpy/fft/README.md | 53 | ||||
| -rw-r--r-- | numpy/fft/__init__.py | 2 | ||||
| -rw-r--r-- | numpy/fft/fftpack.c | 1536 | ||||
| -rw-r--r-- | numpy/fft/fftpack.h | 28 | ||||
| -rw-r--r-- | numpy/fft/fftpack_litemodule.c | 366 | ||||
| -rw-r--r-- | numpy/fft/helper.py | 101 | ||||
| -rw-r--r-- | numpy/fft/pocketfft.c | 2398 | ||||
| -rw-r--r-- | numpy/fft/pocketfft.py (renamed from numpy/fft/fftpack.py) | 104 | ||||
| -rw-r--r-- | numpy/fft/setup.py | 6 | ||||
| -rw-r--r-- | numpy/fft/tests/test_helper.py | 79 | ||||
| -rw-r--r-- | numpy/fft/tests/test_pocketfft.py (renamed from numpy/fft/tests/test_fftpack.py) | 10 |
11 files changed, 2504 insertions, 2179 deletions
diff --git a/numpy/fft/README.md b/numpy/fft/README.md new file mode 100644 index 000000000..7040a2e9b --- /dev/null +++ b/numpy/fft/README.md @@ -0,0 +1,53 @@ +PocketFFT +--------- + +This is a heavily modified implementation of FFTPack [1,2], with the following +advantages: + +- strictly C99 compliant +- more accurate twiddle factor computation +- very fast plan generation +- worst case complexity for transform sizes with large prime factors is + `N*log(N)`, because Bluestein's algorithm [3] is used for these cases. + +License +------- + +3-clause BSD (see LICENSE.md) + + +Some code details +----------------- + +Twiddle factor computation: + +- making use of symmetries to reduce number of sin/cos evaluations +- all angles are reduced to the range `[0; pi/4]` for higher accuracy +- an adapted implementation of `sincospi()` is used, which actually computes + `sin(x)` and `(cos(x)-1)`. +- if `n` sin/cos pairs are required, the adjusted `sincospi()` is only called + `2*sqrt(n)` times; the remaining values are obtained by evaluating the + angle addition theorems in a numerically accurate way. + +Parallel invocation: + +- Plans only contain read-only data; all temporary arrays are allocated and + deallocated during an individual FFT execution. This means that a single plan + can be used in several threads at the same time. + +Efficient codelets are available for the factors: + +- 2, 3, 4, 5, 7, 11 for complex-valued FFTs +- 2, 3, 4, 5 for real-valued FFTs + +Larger prime factors are handled by somewhat less efficient, generic routines. + +For lengths with very large prime factors, Bluestein's algorithm is used, and +instead of an FFT of length `n`, a convolution of length `n2 >= 2*n-1` +is performed, where `n2` is chosen to be highly composite. + + +[1] Swarztrauber, P. 1982, Vectorizing the Fast Fourier Transforms + (New York: Academic Press), 51 +[2] https://www.netlib.org/fftpack/ +[3] https://en.wikipedia.org/wiki/Chirp_Z-transform diff --git a/numpy/fft/__init__.py b/numpy/fft/__init__.py index 44243b483..64b35bc19 100644 --- a/numpy/fft/__init__.py +++ b/numpy/fft/__init__.py @@ -3,7 +3,7 @@ from __future__ import division, absolute_import, print_function # To get sub-modules from .info import __doc__ -from .fftpack import * +from .pocketfft import * from .helper import * from numpy._pytesttester import PytestTester diff --git a/numpy/fft/fftpack.c b/numpy/fft/fftpack.c deleted file mode 100644 index 07fa2bf4c..000000000 --- a/numpy/fft/fftpack.c +++ /dev/null @@ -1,1536 +0,0 @@ -/* - * fftpack.c : A set of FFT routines in C. - * Algorithmically based on Fortran-77 FFTPACK by Paul N. Swarztrauber (Version 4, 1985). -*/ -#define NPY_NO_DEPRECATED_API NPY_API_VERSION - -#include <Python.h> -#include <math.h> -#include <stdio.h> -#include <numpy/ndarraytypes.h> - -#define DOUBLE -#ifdef DOUBLE -#define Treal double -#else -#define Treal float -#endif - -#define ref(u,a) u[a] - -/* Macros for accurate calculation of the twiddle factors. */ -#define TWOPI 6.283185307179586476925286766559005768391 -#define cos2pi(m, n) cos((TWOPI * (m)) / (n)) -#define sin2pi(m, n) sin((TWOPI * (m)) / (n)) - -#define MAXFAC 13 /* maximum number of factors in factorization of n */ -#define NSPECIAL 4 /* number of factors for which we have special-case routines */ - -#ifdef __cplusplus -extern "C" { -#endif - -static void sincos2pi(int m, int n, Treal* si, Treal* co) -/* Calculates sin(2pi * m/n) and cos(2pi * m/n). It is more accurate - * than the naive calculation as the fraction m/n is reduced to [0, 1/8) first. - * Due to the symmetry of sin(x) and cos(x) the values for all x can be - * determined from the function values of the reduced argument in the first - * octant. - */ - { - int n8, m8, octant; - n8 = 8 * n; - m8 = (8 * m) % n8; - octant = m8 / n; - m8 = m8 % n; - switch(octant) { - case 0: - *co = cos2pi(m8, n8); - *si = sin2pi(m8, n8); - break; - case 1: - *co = sin2pi(n-m8, n8); - *si = cos2pi(n-m8, n8); - break; - case 2: - *co = -sin2pi(m8, n8); - *si = cos2pi(m8, n8); - break; - case 3: - *co = -cos2pi(n-m8, n8); - *si = sin2pi(n-m8, n8); - break; - case 4: - *co = -cos2pi(m8, n8); - *si = -sin2pi(m8, n8); - break; - case 5: - *co = -sin2pi(n-m8, n8); - *si = -cos2pi(n-m8, n8); - break; - case 6: - *co = sin2pi(m8, n8); - *si = -cos2pi(m8, n8); - break; - case 7: - *co = cos2pi(n-m8, n8); - *si = -sin2pi(n-m8, n8); - break; - } - } - -/* ---------------------------------------------------------------------- - passf2, passf3, passf4, passf5, passf. Complex FFT passes fwd and bwd. ------------------------------------------------------------------------ */ - -static void passf2(int ido, int l1, const Treal cc[], Treal ch[], const Treal wa1[], int isign) - /* isign==+1 for backward transform */ - { - int i, k, ah, ac; - Treal ti2, tr2; - if (ido <= 2) { - for (k=0; k<l1; k++) { - ah = k*ido; - ac = 2*k*ido; - ch[ah] = ref(cc,ac) + ref(cc,ac + ido); - ch[ah + ido*l1] = ref(cc,ac) - ref(cc,ac + ido); - ch[ah+1] = ref(cc,ac+1) + ref(cc,ac + ido + 1); - ch[ah + ido*l1 + 1] = ref(cc,ac+1) - ref(cc,ac + ido + 1); - } - } else { - for (k=0; k<l1; k++) { - for (i=0; i<ido-1; i+=2) { - ah = i + k*ido; - ac = i + 2*k*ido; - ch[ah] = ref(cc,ac) + ref(cc,ac + ido); - tr2 = ref(cc,ac) - ref(cc,ac + ido); - ch[ah+1] = ref(cc,ac+1) + ref(cc,ac + 1 + ido); - ti2 = ref(cc,ac+1) - ref(cc,ac + 1 + ido); - ch[ah+l1*ido+1] = wa1[i]*ti2 + isign*wa1[i+1]*tr2; - ch[ah+l1*ido] = wa1[i]*tr2 - isign*wa1[i+1]*ti2; - } - } - } - } /* passf2 */ - - -static void passf3(int ido, int l1, const Treal cc[], Treal ch[], - const Treal wa1[], const Treal wa2[], int isign) - /* isign==+1 for backward transform */ - { - static const Treal taur = -0.5; - static const Treal taui = 0.86602540378443864676; - int i, k, ac, ah; - Treal ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2; - if (ido == 2) { - for (k=1; k<=l1; k++) { - ac = (3*k - 2)*ido; - tr2 = ref(cc,ac) + ref(cc,ac + ido); - cr2 = ref(cc,ac - ido) + taur*tr2; - ah = (k - 1)*ido; - ch[ah] = ref(cc,ac - ido) + tr2; - - ti2 = ref(cc,ac + 1) + ref(cc,ac + ido + 1); - ci2 = ref(cc,ac - ido + 1) + taur*ti2; - ch[ah + 1] = ref(cc,ac - ido + 1) + ti2; - - cr3 = isign*taui*(ref(cc,ac) - ref(cc,ac + ido)); - ci3 = isign*taui*(ref(cc,ac + 1) - ref(cc,ac + ido + 1)); - ch[ah + l1*ido] = cr2 - ci3; - ch[ah + 2*l1*ido] = cr2 + ci3; - ch[ah + l1*ido + 1] = ci2 + cr3; - ch[ah + 2*l1*ido + 1] = ci2 - cr3; - } - } else { - for (k=1; k<=l1; k++) { - for (i=0; i<ido-1; i+=2) { - ac = i + (3*k - 2)*ido; - tr2 = ref(cc,ac) + ref(cc,ac + ido); - cr2 = ref(cc,ac - ido) + taur*tr2; - ah = i + (k-1)*ido; - ch[ah] = ref(cc,ac - ido) + tr2; - ti2 = ref(cc,ac + 1) + ref(cc,ac + ido + 1); - ci2 = ref(cc,ac - ido + 1) + taur*ti2; - ch[ah + 1] = ref(cc,ac - ido + 1) + ti2; - cr3 = isign*taui*(ref(cc,ac) - ref(cc,ac + ido)); - ci3 = isign*taui*(ref(cc,ac + 1) - ref(cc,ac + ido + 1)); - dr2 = cr2 - ci3; - dr3 = cr2 + ci3; - di2 = ci2 + cr3; - di3 = ci2 - cr3; - ch[ah + l1*ido + 1] = wa1[i]*di2 + isign*wa1[i+1]*dr2; - ch[ah + l1*ido] = wa1[i]*dr2 - isign*wa1[i+1]*di2; - ch[ah + 2*l1*ido + 1] = wa2[i]*di3 + isign*wa2[i+1]*dr3; - ch[ah + 2*l1*ido] = wa2[i]*dr3 - isign*wa2[i+1]*di3; - } - } - } - } /* passf3 */ - - -static void passf4(int ido, int l1, const Treal cc[], Treal ch[], - const Treal wa1[], const Treal wa2[], const Treal wa3[], int isign) - /* isign == -1 for forward transform and +1 for backward transform */ - { - int i, k, ac, ah; - Treal ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4; - if (ido == 2) { - for (k=0; k<l1; k++) { - ac = 4*k*ido + 1; - ti1 = ref(cc,ac) - ref(cc,ac + 2*ido); - ti2 = ref(cc,ac) + ref(cc,ac + 2*ido); - tr4 = ref(cc,ac + 3*ido) - ref(cc,ac + ido); - ti3 = ref(cc,ac + ido) + ref(cc,ac + 3*ido); - tr1 = ref(cc,ac - 1) - ref(cc,ac + 2*ido - 1); - tr2 = ref(cc,ac - 1) + ref(cc,ac + 2*ido - 1); - ti4 = ref(cc,ac + ido - 1) - ref(cc,ac + 3*ido - 1); - tr3 = ref(cc,ac + ido - 1) + ref(cc,ac + 3*ido - 1); - ah = k*ido; - ch[ah] = tr2 + tr3; - ch[ah + 2*l1*ido] = tr2 - tr3; - ch[ah + 1] = ti2 + ti3; - ch[ah + 2*l1*ido + 1] = ti2 - ti3; - ch[ah + l1*ido] = tr1 + isign*tr4; - ch[ah + 3*l1*ido] = tr1 - isign*tr4; - ch[ah + l1*ido + 1] = ti1 + isign*ti4; - ch[ah + 3*l1*ido + 1] = ti1 - isign*ti4; - } - } else { - for (k=0; k<l1; k++) { - for (i=0; i<ido-1; i+=2) { - ac = i + 1 + 4*k*ido; - ti1 = ref(cc,ac) - ref(cc,ac + 2*ido); - ti2 = ref(cc,ac) + ref(cc,ac + 2*ido); - ti3 = ref(cc,ac + ido) + ref(cc,ac + 3*ido); - tr4 = ref(cc,ac + 3*ido) - ref(cc,ac + ido); - tr1 = ref(cc,ac - 1) - ref(cc,ac + 2*ido - 1); - tr2 = ref(cc,ac - 1) + ref(cc,ac + 2*ido - 1); - ti4 = ref(cc,ac + ido - 1) - ref(cc,ac + 3*ido - 1); - tr3 = ref(cc,ac + ido - 1) + ref(cc,ac + 3*ido - 1); - ah = i + k*ido; - ch[ah] = tr2 + tr3; - cr3 = tr2 - tr3; - ch[ah + 1] = ti2 + ti3; - ci3 = ti2 - ti3; - cr2 = tr1 + isign*tr4; - cr4 = tr1 - isign*tr4; - ci2 = ti1 + isign*ti4; - ci4 = ti1 - isign*ti4; - ch[ah + l1*ido] = wa1[i]*cr2 - isign*wa1[i + 1]*ci2; - ch[ah + l1*ido + 1] = wa1[i]*ci2 + isign*wa1[i + 1]*cr2; - ch[ah + 2*l1*ido] = wa2[i]*cr3 - isign*wa2[i + 1]*ci3; - ch[ah + 2*l1*ido + 1] = wa2[i]*ci3 + isign*wa2[i + 1]*cr3; - ch[ah + 3*l1*ido] = wa3[i]*cr4 -isign*wa3[i + 1]*ci4; - ch[ah + 3*l1*ido + 1] = wa3[i]*ci4 + isign*wa3[i + 1]*cr4; - } - } - } - } /* passf4 */ - - -static void passf5(int ido, int l1, const Treal cc[], Treal ch[], - const Treal wa1[], const Treal wa2[], const Treal wa3[], const Treal wa4[], int isign) - /* isign == -1 for forward transform and +1 for backward transform */ - { - static const Treal tr11 = 0.3090169943749474241; - static const Treal ti11 = 0.95105651629515357212; - static const Treal tr12 = -0.8090169943749474241; - static const Treal ti12 = 0.58778525229247312917; - int i, k, ac, ah; - Treal ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5, cr4, ti2, ti3, - ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5; - if (ido == 2) { - for (k = 1; k <= l1; ++k) { - ac = (5*k - 4)*ido + 1; - ti5 = ref(cc,ac) - ref(cc,ac + 3*ido); - ti2 = ref(cc,ac) + ref(cc,ac + 3*ido); - ti4 = ref(cc,ac + ido) - ref(cc,ac + 2*ido); - ti3 = ref(cc,ac + ido) + ref(cc,ac + 2*ido); - tr5 = ref(cc,ac - 1) - ref(cc,ac + 3*ido - 1); - tr2 = ref(cc,ac - 1) + ref(cc,ac + 3*ido - 1); - tr4 = ref(cc,ac + ido - 1) - ref(cc,ac + 2*ido - 1); - tr3 = ref(cc,ac + ido - 1) + ref(cc,ac + 2*ido - 1); - ah = (k - 1)*ido; - ch[ah] = ref(cc,ac - ido - 1) + tr2 + tr3; - ch[ah + 1] = ref(cc,ac - ido) + ti2 + ti3; - cr2 = ref(cc,ac - ido - 1) + tr11*tr2 + tr12*tr3; - ci2 = ref(cc,ac - ido) + tr11*ti2 + tr12*ti3; - cr3 = ref(cc,ac - ido - 1) + tr12*tr2 + tr11*tr3; - ci3 = ref(cc,ac - ido) + tr12*ti2 + tr11*ti3; - cr5 = isign*(ti11*tr5 + ti12*tr4); - ci5 = isign*(ti11*ti5 + ti12*ti4); - cr4 = isign*(ti12*tr5 - ti11*tr4); - ci4 = isign*(ti12*ti5 - ti11*ti4); - ch[ah + l1*ido] = cr2 - ci5; - ch[ah + 4*l1*ido] = cr2 + ci5; - ch[ah + l1*ido + 1] = ci2 + cr5; - ch[ah + 2*l1*ido + 1] = ci3 + cr4; - ch[ah + 2*l1*ido] = cr3 - ci4; - ch[ah + 3*l1*ido] = cr3 + ci4; - ch[ah + 3*l1*ido + 1] = ci3 - cr4; - ch[ah + 4*l1*ido + 1] = ci2 - cr5; - } - } else { - for (k=1; k<=l1; k++) { - for (i=0; i<ido-1; i+=2) { - ac = i + 1 + (k*5 - 4)*ido; - ti5 = ref(cc,ac) - ref(cc,ac + 3*ido); - ti2 = ref(cc,ac) + ref(cc,ac + 3*ido); - ti4 = ref(cc,ac + ido) - ref(cc,ac + 2*ido); - ti3 = ref(cc,ac + ido) + ref(cc,ac + 2*ido); - tr5 = ref(cc,ac - 1) - ref(cc,ac + 3*ido - 1); - tr2 = ref(cc,ac - 1) + ref(cc,ac + 3*ido - 1); - tr4 = ref(cc,ac + ido - 1) - ref(cc,ac + 2*ido - 1); - tr3 = ref(cc,ac + ido - 1) + ref(cc,ac + 2*ido - 1); - ah = i + (k - 1)*ido; - ch[ah] = ref(cc,ac - ido - 1) + tr2 + tr3; - ch[ah + 1] = ref(cc,ac - ido) + ti2 + ti3; - cr2 = ref(cc,ac - ido - 1) + tr11*tr2 + tr12*tr3; - - ci2 = ref(cc,ac - ido) + tr11*ti2 + tr12*ti3; - cr3 = ref(cc,ac - ido - 1) + tr12*tr2 + tr11*tr3; - - ci3 = ref(cc,ac - ido) + tr12*ti2 + tr11*ti3; - cr5 = isign*(ti11*tr5 + ti12*tr4); - ci5 = isign*(ti11*ti5 + ti12*ti4); - cr4 = isign*(ti12*tr5 - ti11*tr4); - ci4 = isign*(ti12*ti5 - ti11*ti4); - dr3 = cr3 - ci4; - dr4 = cr3 + ci4; - di3 = ci3 + cr4; - di4 = ci3 - cr4; - dr5 = cr2 + ci5; - dr2 = cr2 - ci5; - di5 = ci2 - cr5; - di2 = ci2 + cr5; - ch[ah + l1*ido] = wa1[i]*dr2 - isign*wa1[i+1]*di2; - ch[ah + l1*ido + 1] = wa1[i]*di2 + isign*wa1[i+1]*dr2; - ch[ah + 2*l1*ido] = wa2[i]*dr3 - isign*wa2[i+1]*di3; - ch[ah + 2*l1*ido + 1] = wa2[i]*di3 + isign*wa2[i+1]*dr3; - ch[ah + 3*l1*ido] = wa3[i]*dr4 - isign*wa3[i+1]*di4; - ch[ah + 3*l1*ido + 1] = wa3[i]*di4 + isign*wa3[i+1]*dr4; - ch[ah + 4*l1*ido] = wa4[i]*dr5 - isign*wa4[i+1]*di5; - ch[ah + 4*l1*ido + 1] = wa4[i]*di5 + isign*wa4[i+1]*dr5; - } - } - } - } /* passf5 */ - - -static void passf(int *nac, int ido, int ip, int l1, int idl1, - Treal cc[], Treal ch[], - const Treal wa[], int isign) - /* isign is -1 for forward transform and +1 for backward transform */ - { - int idij, idlj, idot, ipph, i, j, k, l, jc, lc, ik, idj, idl, inc,idp; - Treal wai, war; - - idot = ido / 2; - /* nt = ip*idl1;*/ - ipph = (ip + 1) / 2; - idp = ip*ido; - if (ido >= l1) { - for (j=1; j<ipph; j++) { - jc = ip - j; - for (k=0; k<l1; k++) { - for (i=0; i<ido; i++) { - ch[i + (k + j*l1)*ido] = - ref(cc,i + (j + k*ip)*ido) + ref(cc,i + (jc + k*ip)*ido); - ch[i + (k + jc*l1)*ido] = - ref(cc,i + (j + k*ip)*ido) - ref(cc,i + (jc + k*ip)*ido); - } - } - } - for (k=0; k<l1; k++) - for (i=0; i<ido; i++) - ch[i + k*ido] = ref(cc,i + k*ip*ido); - } else { - for (j=1; j<ipph; j++) { - jc = ip - j; - for (i=0; i<ido; i++) { - for (k=0; k<l1; k++) { - ch[i + (k + j*l1)*ido] = ref(cc,i + (j + k*ip)*ido) + ref(cc,i + (jc + k* - ip)*ido); - ch[i + (k + jc*l1)*ido] = ref(cc,i + (j + k*ip)*ido) - ref(cc,i + (jc + k* - ip)*ido); - } - } - } - for (i=0; i<ido; i++) - for (k=0; k<l1; k++) - ch[i + k*ido] = ref(cc,i + k*ip*ido); - } - - idl = 2 - ido; - inc = 0; - for (l=1; l<ipph; l++) { - lc = ip - l; - idl += ido; - for (ik=0; ik<idl1; ik++) { - cc[ik + l*idl1] = ch[ik] + wa[idl - 2]*ch[ik + idl1]; - cc[ik + lc*idl1] = isign*wa[idl-1]*ch[ik + (ip-1)*idl1]; - } - idlj = idl; - inc += ido; - for (j=2; j<ipph; j++) { - jc = ip - j; - idlj += inc; - if (idlj > idp) idlj -= idp; - war = wa[idlj - 2]; - wai = wa[idlj-1]; - for (ik=0; ik<idl1; ik++) { - cc[ik + l*idl1] += war*ch[ik + j*idl1]; - cc[ik + lc*idl1] += isign*wai*ch[ik + jc*idl1]; - } - } - } - for (j=1; j<ipph; j++) - for (ik=0; ik<idl1; ik++) - ch[ik] += ch[ik + j*idl1]; - for (j=1; j<ipph; j++) { - jc = ip - j; - for (ik=1; ik<idl1; ik+=2) { - ch[ik - 1 + j*idl1] = cc[ik - 1 + j*idl1] - cc[ik + jc*idl1]; - ch[ik - 1 + jc*idl1] = cc[ik - 1 + j*idl1] + cc[ik + jc*idl1]; - ch[ik + j*idl1] = cc[ik + j*idl1] + cc[ik - 1 + jc*idl1]; - ch[ik + jc*idl1] = cc[ik + j*idl1] - cc[ik - 1 + jc*idl1]; - } - } - *nac = 1; - if (ido == 2) return; - *nac = 0; - for (ik=0; ik<idl1; ik++) - cc[ik] = ch[ik]; - for (j=1; j<ip; j++) { - for (k=0; k<l1; k++) { - cc[(k + j*l1)*ido + 0] = ch[(k + j*l1)*ido + 0]; - cc[(k + j*l1)*ido + 1] = ch[(k + j*l1)*ido + 1]; - } - } - if (idot <= l1) { - idij = 0; - for (j=1; j<ip; j++) { - idij += 2; - for (i=3; i<ido; i+=2) { - idij += 2; - for (k=0; k<l1; k++) { - cc[i - 1 + (k + j*l1)*ido] = - wa[idij - 2]*ch[i - 1 + (k + j*l1)*ido] - - isign*wa[idij-1]*ch[i + (k + j*l1)*ido]; - cc[i + (k + j*l1)*ido] = - wa[idij - 2]*ch[i + (k + j*l1)*ido] + - isign*wa[idij-1]*ch[i - 1 + (k + j*l1)*ido]; - } - } - } - } else { - idj = 2 - ido; - for (j=1; j<ip; j++) { - idj += ido; - for (k = 0; k < l1; k++) { - idij = idj; - for (i=3; i<ido; i+=2) { - idij += 2; - cc[i - 1 + (k + j*l1)*ido] = - wa[idij - 2]*ch[i - 1 + (k + j*l1)*ido] - - isign*wa[idij-1]*ch[i + (k + j*l1)*ido]; - cc[i + (k + j*l1)*ido] = - wa[idij - 2]*ch[i + (k + j*l1)*ido] + - isign*wa[idij-1]*ch[i - 1 + (k + j*l1)*ido]; - } - } - } - } - } /* passf */ - - - /* ---------------------------------------------------------------------- -radf2,radb2, radf3,radb3, radf4,radb4, radf5,radb5, radfg,radbg. -Treal FFT passes fwd and bwd. ----------------------------------------------------------------------- */ - -static void radf2(int ido, int l1, const Treal cc[], Treal ch[], const Treal wa1[]) - { - int i, k, ic; - Treal ti2, tr2; - for (k=0; k<l1; k++) { - ch[2*k*ido] = - ref(cc,k*ido) + ref(cc,(k + l1)*ido); - ch[(2*k+1)*ido + ido-1] = - ref(cc,k*ido) - ref(cc,(k + l1)*ido); - } - if (ido < 2) return; - if (ido != 2) { - for (k=0; k<l1; k++) { - for (i=2; i<ido; i+=2) { - ic = ido - i; - tr2 = wa1[i - 2]*ref(cc, i-1 + (k + l1)*ido) + wa1[i - 1]*ref(cc, i + (k + l1)*ido); - ti2 = wa1[i - 2]*ref(cc, i + (k + l1)*ido) - wa1[i - 1]*ref(cc, i-1 + (k + l1)*ido); - ch[i + 2*k*ido] = ref(cc,i + k*ido) + ti2; - ch[ic + (2*k+1)*ido] = ti2 - ref(cc,i + k*ido); - ch[i - 1 + 2*k*ido] = ref(cc,i - 1 + k*ido) + tr2; - ch[ic - 1 + (2*k+1)*ido] = ref(cc,i - 1 + k*ido) - tr2; - } - } - if (ido % 2 == 1) return; - } - for (k=0; k<l1; k++) { - ch[(2*k+1)*ido] = -ref(cc,ido-1 + (k + l1)*ido); - ch[ido-1 + 2*k*ido] = ref(cc,ido-1 + k*ido); - } - } /* radf2 */ - - -static void radb2(int ido, int l1, const Treal cc[], Treal ch[], const Treal wa1[]) - { - int i, k, ic; - Treal ti2, tr2; - for (k=0; k<l1; k++) { - ch[k*ido] = - ref(cc,2*k*ido) + ref(cc,ido-1 + (2*k+1)*ido); - ch[(k + l1)*ido] = - ref(cc,2*k*ido) - ref(cc,ido-1 + (2*k+1)*ido); - } - if (ido < 2) return; - if (ido != 2) { - for (k = 0; k < l1; ++k) { - for (i = 2; i < ido; i += 2) { - ic = ido - i; - ch[i-1 + k*ido] = - ref(cc,i-1 + 2*k*ido) + ref(cc,ic-1 + (2*k+1)*ido); - tr2 = ref(cc,i-1 + 2*k*ido) - ref(cc,ic-1 + (2*k+1)*ido); - ch[i + k*ido] = - ref(cc,i + 2*k*ido) - ref(cc,ic + (2*k+1)*ido); - ti2 = ref(cc,i + (2*k)*ido) + ref(cc,ic + (2*k+1)*ido); - ch[i-1 + (k + l1)*ido] = - wa1[i - 2]*tr2 - wa1[i - 1]*ti2; - ch[i + (k + l1)*ido] = - wa1[i - 2]*ti2 + wa1[i - 1]*tr2; - } - } - if (ido % 2 == 1) return; - } - for (k = 0; k < l1; k++) { - ch[ido-1 + k*ido] = 2*ref(cc,ido-1 + 2*k*ido); - ch[ido-1 + (k + l1)*ido] = -2*ref(cc,(2*k+1)*ido); - } - } /* radb2 */ - - -static void radf3(int ido, int l1, const Treal cc[], Treal ch[], - const Treal wa1[], const Treal wa2[]) - { - static const Treal taur = -0.5; - static const Treal taui = 0.86602540378443864676; - int i, k, ic; - Treal ci2, di2, di3, cr2, dr2, dr3, ti2, ti3, tr2, tr3; - for (k=0; k<l1; k++) { - cr2 = ref(cc,(k + l1)*ido) + ref(cc,(k + 2*l1)*ido); - ch[3*k*ido] = ref(cc,k*ido) + cr2; - ch[(3*k+2)*ido] = taui*(ref(cc,(k + l1*2)*ido) - ref(cc,(k + l1)*ido)); - ch[ido-1 + (3*k + 1)*ido] = ref(cc,k*ido) + taur*cr2; - } - if (ido == 1) return; - for (k=0; k<l1; k++) { - for (i=2; i<ido; i+=2) { - ic = ido - i; - dr2 = wa1[i - 2]*ref(cc,i - 1 + (k + l1)*ido) + - wa1[i - 1]*ref(cc,i + (k + l1)*ido); - di2 = wa1[i - 2]*ref(cc,i + (k + l1)*ido) - wa1[i - 1]*ref(cc,i - 1 + (k + l1)*ido); - dr3 = wa2[i - 2]*ref(cc,i - 1 + (k + l1*2)*ido) + wa2[i - 1]*ref(cc,i + (k + l1*2)*ido); - di3 = wa2[i - 2]*ref(cc,i + (k + l1*2)*ido) - wa2[i - 1]*ref(cc,i - 1 + (k + l1*2)*ido); - cr2 = dr2 + dr3; - ci2 = di2 + di3; - ch[i - 1 + 3*k*ido] = ref(cc,i - 1 + k*ido) + cr2; - ch[i + 3*k*ido] = ref(cc,i + k*ido) + ci2; - tr2 = ref(cc,i - 1 + k*ido) + taur*cr2; - ti2 = ref(cc,i + k*ido) + taur*ci2; - tr3 = taui*(di2 - di3); - ti3 = taui*(dr3 - dr2); - ch[i - 1 + (3*k + 2)*ido] = tr2 + tr3; - ch[ic - 1 + (3*k + 1)*ido] = tr2 - tr3; - ch[i + (3*k + 2)*ido] = ti2 + ti3; - ch[ic + (3*k + 1)*ido] = ti3 - ti2; - } - } - } /* radf3 */ - - -static void radb3(int ido, int l1, const Treal cc[], Treal ch[], - const Treal wa1[], const Treal wa2[]) - { - static const Treal taur = -0.5; - static const Treal taui = 0.86602540378443864676; - int i, k, ic; - Treal ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2; - for (k=0; k<l1; k++) { - tr2 = 2*ref(cc,ido-1 + (3*k + 1)*ido); - cr2 = ref(cc,3*k*ido) + taur*tr2; - ch[k*ido] = ref(cc,3*k*ido) + tr2; - ci3 = 2*taui*ref(cc,(3*k + 2)*ido); - ch[(k + l1)*ido] = cr2 - ci3; - ch[(k + 2*l1)*ido] = cr2 + ci3; - } - if (ido == 1) return; - for (k=0; k<l1; k++) { - for (i=2; i<ido; i+=2) { - ic = ido - i; - tr2 = ref(cc,i - 1 + (3*k + 2)*ido) + ref(cc,ic - 1 + (3*k + 1)*ido); - cr2 = ref(cc,i - 1 + 3*k*ido) + taur*tr2; - ch[i - 1 + k*ido] = ref(cc,i - 1 + 3*k*ido) + tr2; - ti2 = ref(cc,i + (3*k + 2)*ido) - ref(cc,ic + (3*k + 1)*ido); - ci2 = ref(cc,i + 3*k*ido) + taur*ti2; - ch[i + k*ido] = ref(cc,i + 3*k*ido) + ti2; - cr3 = taui*(ref(cc,i - 1 + (3*k + 2)*ido) - ref(cc,ic - 1 + (3*k + 1)*ido)); - ci3 = taui*(ref(cc,i + (3*k + 2)*ido) + ref(cc,ic + (3*k + 1)*ido)); - dr2 = cr2 - ci3; - dr3 = cr2 + ci3; - di2 = ci2 + cr3; - di3 = ci2 - cr3; - ch[i - 1 + (k + l1)*ido] = wa1[i - 2]*dr2 - wa1[i - 1]*di2; - ch[i + (k + l1)*ido] = wa1[i - 2]*di2 + wa1[i - 1]*dr2; - ch[i - 1 + (k + 2*l1)*ido] = wa2[i - 2]*dr3 - wa2[i - 1]*di3; - ch[i + (k + 2*l1)*ido] = wa2[i - 2]*di3 + wa2[i - 1]*dr3; - } - } - } /* radb3 */ - - -static void radf4(int ido, int l1, const Treal cc[], Treal ch[], - const Treal wa1[], const Treal wa2[], const Treal wa3[]) - { - static const Treal hsqt2 = 0.70710678118654752440; - int i, k, ic; - Treal ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4; - for (k=0; k<l1; k++) { - tr1 = ref(cc,(k + l1)*ido) + ref(cc,(k + 3*l1)*ido); - tr2 = ref(cc,k*ido) + ref(cc,(k + 2*l1)*ido); - ch[4*k*ido] = tr1 + tr2; - ch[ido-1 + (4*k + 3)*ido] = tr2 - tr1; - ch[ido-1 + (4*k + 1)*ido] = ref(cc,k*ido) - ref(cc,(k + 2*l1)*ido); - ch[(4*k + 2)*ido] = ref(cc,(k + 3*l1)*ido) - ref(cc,(k + l1)*ido); - } - if (ido < 2) return; - if (ido != 2) { - for (k=0; k<l1; k++) { - for (i=2; i<ido; i += 2) { - ic = ido - i; - cr2 = wa1[i - 2]*ref(cc,i - 1 + (k + l1)*ido) + wa1[i - 1]*ref(cc,i + (k + l1)*ido); - ci2 = wa1[i - 2]*ref(cc,i + (k + l1)*ido) - wa1[i - 1]*ref(cc,i - 1 + (k + l1)*ido); - cr3 = wa2[i - 2]*ref(cc,i - 1 + (k + 2*l1)*ido) + wa2[i - 1]*ref(cc,i + (k + 2*l1)* - ido); - ci3 = wa2[i - 2]*ref(cc,i + (k + 2*l1)*ido) - wa2[i - 1]*ref(cc,i - 1 + (k + 2*l1)* - ido); - cr4 = wa3[i - 2]*ref(cc,i - 1 + (k + 3*l1)*ido) + wa3[i - 1]*ref(cc,i + (k + 3*l1)* - ido); - ci4 = wa3[i - 2]*ref(cc,i + (k + 3*l1)*ido) - wa3[i - 1]*ref(cc,i - 1 + (k + 3*l1)* - ido); - tr1 = cr2 + cr4; - tr4 = cr4 - cr2; - ti1 = ci2 + ci4; - ti4 = ci2 - ci4; - ti2 = ref(cc,i + k*ido) + ci3; - ti3 = ref(cc,i + k*ido) - ci3; - tr2 = ref(cc,i - 1 + k*ido) + cr3; - tr3 = ref(cc,i - 1 + k*ido) - cr3; - ch[i - 1 + 4*k*ido] = tr1 + tr2; - ch[ic - 1 + (4*k + 3)*ido] = tr2 - tr1; - ch[i + 4*k*ido] = ti1 + ti2; - ch[ic + (4*k + 3)*ido] = ti1 - ti2; - ch[i - 1 + (4*k + 2)*ido] = ti4 + tr3; - ch[ic - 1 + (4*k + 1)*ido] = tr3 - ti4; - ch[i + (4*k + 2)*ido] = tr4 + ti3; - ch[ic + (4*k + 1)*ido] = tr4 - ti3; - } - } - if (ido % 2 == 1) return; - } - for (k=0; k<l1; k++) { - ti1 = -hsqt2*(ref(cc,ido-1 + (k + l1)*ido) + ref(cc,ido-1 + (k + 3*l1)*ido)); - tr1 = hsqt2*(ref(cc,ido-1 + (k + l1)*ido) - ref(cc,ido-1 + (k + 3*l1)*ido)); - ch[ido-1 + 4*k*ido] = tr1 + ref(cc,ido-1 + k*ido); - ch[ido-1 + (4*k + 2)*ido] = ref(cc,ido-1 + k*ido) - tr1; - ch[(4*k + 1)*ido] = ti1 - ref(cc,ido-1 + (k + 2*l1)*ido); - ch[(4*k + 3)*ido] = ti1 + ref(cc,ido-1 + (k + 2*l1)*ido); - } - } /* radf4 */ - - -static void radb4(int ido, int l1, const Treal cc[], Treal ch[], - const Treal wa1[], const Treal wa2[], const Treal wa3[]) - { - static const Treal sqrt2 = 1.41421356237309504880; - int i, k, ic; - Treal ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4; - for (k = 0; k < l1; k++) { - tr1 = ref(cc,4*k*ido) - ref(cc,ido-1 + (4*k + 3)*ido); - tr2 = ref(cc,4*k*ido) + ref(cc,ido-1 + (4*k + 3)*ido); - tr3 = ref(cc,ido-1 + (4*k + 1)*ido) + ref(cc,ido-1 + (4*k + 1)*ido); - tr4 = ref(cc,(4*k + 2)*ido) + ref(cc,(4*k + 2)*ido); - ch[k*ido] = tr2 + tr3; - ch[(k + l1)*ido] = tr1 - tr4; - ch[(k + 2*l1)*ido] = tr2 - tr3; - ch[(k + 3*l1)*ido] = tr1 + tr4; - } - if (ido < 2) return; - if (ido != 2) { - for (k = 0; k < l1; ++k) { - for (i = 2; i < ido; i += 2) { - ic = ido - i; - ti1 = ref(cc,i + 4*k*ido) + ref(cc,ic + (4*k + 3)*ido); - ti2 = ref(cc,i + 4*k*ido) - ref(cc,ic + (4*k + 3)*ido); - ti3 = ref(cc,i + (4*k + 2)*ido) - ref(cc,ic + (4*k + 1)*ido); - tr4 = ref(cc,i + (4*k + 2)*ido) + ref(cc,ic + (4*k + 1)*ido); - tr1 = ref(cc,i - 1 + 4*k*ido) - ref(cc,ic - 1 + (4*k + 3)*ido); - tr2 = ref(cc,i - 1 + 4*k*ido) + ref(cc,ic - 1 + (4*k + 3)*ido); - ti4 = ref(cc,i - 1 + (4*k + 2)*ido) - ref(cc,ic - 1 + (4*k + 1)*ido); - tr3 = ref(cc,i - 1 + (4*k + 2)*ido) + ref(cc,ic - 1 + (4*k + 1)*ido); - ch[i - 1 + k*ido] = tr2 + tr3; - cr3 = tr2 - tr3; - ch[i + k*ido] = ti2 + ti3; - ci3 = ti2 - ti3; - cr2 = tr1 - tr4; - cr4 = tr1 + tr4; - ci2 = ti1 + ti4; - ci4 = ti1 - ti4; - ch[i - 1 + (k + l1)*ido] = wa1[i - 2]*cr2 - wa1[i - 1]*ci2; - ch[i + (k + l1)*ido] = wa1[i - 2]*ci2 + wa1[i - 1]*cr2; - ch[i - 1 + (k + 2*l1)*ido] = wa2[i - 2]*cr3 - wa2[i - 1]*ci3; - ch[i + (k + 2*l1)*ido] = wa2[i - 2]*ci3 + wa2[i - 1]*cr3; - ch[i - 1 + (k + 3*l1)*ido] = wa3[i - 2]*cr4 - wa3[i - 1]*ci4; - ch[i + (k + 3*l1)*ido] = wa3[i - 2]*ci4 + wa3[i - 1]*cr4; - } - } - if (ido % 2 == 1) return; - } - for (k = 0; k < l1; k++) { - ti1 = ref(cc,(4*k + 1)*ido) + ref(cc,(4*k + 3)*ido); - ti2 = ref(cc,(4*k + 3)*ido) - ref(cc,(4*k + 1)*ido); - tr1 = ref(cc,ido-1 + 4*k*ido) - ref(cc,ido-1 + (4*k + 2)*ido); - tr2 = ref(cc,ido-1 + 4*k*ido) + ref(cc,ido-1 + (4*k + 2)*ido); - ch[ido-1 + k*ido] = tr2 + tr2; - ch[ido-1 + (k + l1)*ido] = sqrt2*(tr1 - ti1); - ch[ido-1 + (k + 2*l1)*ido] = ti2 + ti2; - ch[ido-1 + (k + 3*l1)*ido] = -sqrt2*(tr1 + ti1); - } - } /* radb4 */ - - -static void radf5(int ido, int l1, const Treal cc[], Treal ch[], - const Treal wa1[], const Treal wa2[], const Treal wa3[], const Treal wa4[]) - { - static const Treal tr11 = 0.3090169943749474241; - static const Treal ti11 = 0.95105651629515357212; - static const Treal tr12 = -0.8090169943749474241; - static const Treal ti12 = 0.58778525229247312917; - int i, k, ic; - Treal ci2, di2, ci4, ci5, di3, di4, di5, ci3, cr2, cr3, dr2, dr3, dr4, dr5, - cr5, cr4, ti2, ti3, ti5, ti4, tr2, tr3, tr4, tr5; - for (k = 0; k < l1; k++) { - cr2 = ref(cc,(k + 4*l1)*ido) + ref(cc,(k + l1)*ido); - ci5 = ref(cc,(k + 4*l1)*ido) - ref(cc,(k + l1)*ido); - cr3 = ref(cc,(k + 3*l1)*ido) + ref(cc,(k + 2*l1)*ido); - ci4 = ref(cc,(k + 3*l1)*ido) - ref(cc,(k + 2*l1)*ido); - ch[5*k*ido] = ref(cc,k*ido) + cr2 + cr3; - ch[ido-1 + (5*k + 1)*ido] = ref(cc,k*ido) + tr11*cr2 + tr12*cr3; - ch[(5*k + 2)*ido] = ti11*ci5 + ti12*ci4; - ch[ido-1 + (5*k + 3)*ido] = ref(cc,k*ido) + tr12*cr2 + tr11*cr3; - ch[(5*k + 4)*ido] = ti12*ci5 - ti11*ci4; - } - if (ido == 1) return; - for (k = 0; k < l1; ++k) { - for (i = 2; i < ido; i += 2) { - ic = ido - i; - dr2 = wa1[i - 2]*ref(cc,i - 1 + (k + l1)*ido) + wa1[i - 1]*ref(cc,i + (k + l1)*ido); - di2 = wa1[i - 2]*ref(cc,i + (k + l1)*ido) - wa1[i - 1]*ref(cc,i - 1 + (k + l1)*ido); - dr3 = wa2[i - 2]*ref(cc,i - 1 + (k + 2*l1)*ido) + wa2[i - 1]*ref(cc,i + (k + 2*l1)*ido); - di3 = wa2[i - 2]*ref(cc,i + (k + 2*l1)*ido) - wa2[i - 1]*ref(cc,i - 1 + (k + 2*l1)*ido); - dr4 = wa3[i - 2]*ref(cc,i - 1 + (k + 3*l1)*ido) + wa3[i - 1]*ref(cc,i + (k + 3*l1)*ido); - di4 = wa3[i - 2]*ref(cc,i + (k + 3*l1)*ido) - wa3[i - 1]*ref(cc,i - 1 + (k + 3*l1)*ido); - dr5 = wa4[i - 2]*ref(cc,i - 1 + (k + 4*l1)*ido) + wa4[i - 1]*ref(cc,i + (k + 4*l1)*ido); - di5 = wa4[i - 2]*ref(cc,i + (k + 4*l1)*ido) - wa4[i - 1]*ref(cc,i - 1 + (k + 4*l1)*ido); - cr2 = dr2 + dr5; - ci5 = dr5 - dr2; - cr5 = di2 - di5; - ci2 = di2 + di5; - cr3 = dr3 + dr4; - ci4 = dr4 - dr3; - cr4 = di3 - di4; - ci3 = di3 + di4; - ch[i - 1 + 5*k*ido] = ref(cc,i - 1 + k*ido) + cr2 + cr3; - ch[i + 5*k*ido] = ref(cc,i + k*ido) + ci2 + ci3; - tr2 = ref(cc,i - 1 + k*ido) + tr11*cr2 + tr12*cr3; - ti2 = ref(cc,i + k*ido) + tr11*ci2 + tr12*ci3; - tr3 = ref(cc,i - 1 + k*ido) + tr12*cr2 + tr11*cr3; - ti3 = ref(cc,i + k*ido) + tr12*ci2 + tr11*ci3; - tr5 = ti11*cr5 + ti12*cr4; - ti5 = ti11*ci5 + ti12*ci4; - tr4 = ti12*cr5 - ti11*cr4; - ti4 = ti12*ci5 - ti11*ci4; - ch[i - 1 + (5*k + 2)*ido] = tr2 + tr5; - ch[ic - 1 + (5*k + 1)*ido] = tr2 - tr5; - ch[i + (5*k + 2)*ido] = ti2 + ti5; - ch[ic + (5*k + 1)*ido] = ti5 - ti2; - ch[i - 1 + (5*k + 4)*ido] = tr3 + tr4; - ch[ic - 1 + (5*k + 3)*ido] = tr3 - tr4; - ch[i + (5*k + 4)*ido] = ti3 + ti4; - ch[ic + (5*k + 3)*ido] = ti4 - ti3; - } - } - } /* radf5 */ - - -static void radb5(int ido, int l1, const Treal cc[], Treal ch[], - const Treal wa1[], const Treal wa2[], const Treal wa3[], const Treal wa4[]) - { - static const Treal tr11 = 0.3090169943749474241; - static const Treal ti11 = 0.95105651629515357212; - static const Treal tr12 = -0.8090169943749474241; - static const Treal ti12 = 0.58778525229247312917; - int i, k, ic; - Treal ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5, cr4, ti2, ti3, - ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5; - for (k = 0; k < l1; k++) { - ti5 = 2*ref(cc,(5*k + 2)*ido); - ti4 = 2*ref(cc,(5*k + 4)*ido); - tr2 = 2*ref(cc,ido-1 + (5*k + 1)*ido); - tr3 = 2*ref(cc,ido-1 + (5*k + 3)*ido); - ch[k*ido] = ref(cc,5*k*ido) + tr2 + tr3; - cr2 = ref(cc,5*k*ido) + tr11*tr2 + tr12*tr3; - cr3 = ref(cc,5*k*ido) + tr12*tr2 + tr11*tr3; - ci5 = ti11*ti5 + ti12*ti4; - ci4 = ti12*ti5 - ti11*ti4; - ch[(k + l1)*ido] = cr2 - ci5; - ch[(k + 2*l1)*ido] = cr3 - ci4; - ch[(k + 3*l1)*ido] = cr3 + ci4; - ch[(k + 4*l1)*ido] = cr2 + ci5; - } - if (ido == 1) return; - for (k = 0; k < l1; ++k) { - for (i = 2; i < ido; i += 2) { - ic = ido - i; - ti5 = ref(cc,i + (5*k + 2)*ido) + ref(cc,ic + (5*k + 1)*ido); - ti2 = ref(cc,i + (5*k + 2)*ido) - ref(cc,ic + (5*k + 1)*ido); - ti4 = ref(cc,i + (5*k + 4)*ido) + ref(cc,ic + (5*k + 3)*ido); - ti3 = ref(cc,i + (5*k + 4)*ido) - ref(cc,ic + (5*k + 3)*ido); - tr5 = ref(cc,i - 1 + (5*k + 2)*ido) - ref(cc,ic - 1 + (5*k + 1)*ido); - tr2 = ref(cc,i - 1 + (5*k + 2)*ido) + ref(cc,ic - 1 + (5*k + 1)*ido); - tr4 = ref(cc,i - 1 + (5*k + 4)*ido) - ref(cc,ic - 1 + (5*k + 3)*ido); - tr3 = ref(cc,i - 1 + (5*k + 4)*ido) + ref(cc,ic - 1 + (5*k + 3)*ido); - ch[i - 1 + k*ido] = ref(cc,i - 1 + 5*k*ido) + tr2 + tr3; - ch[i + k*ido] = ref(cc,i + 5*k*ido) + ti2 + ti3; - cr2 = ref(cc,i - 1 + 5*k*ido) + tr11*tr2 + tr12*tr3; - - ci2 = ref(cc,i + 5*k*ido) + tr11*ti2 + tr12*ti3; - cr3 = ref(cc,i - 1 + 5*k*ido) + tr12*tr2 + tr11*tr3; - - ci3 = ref(cc,i + 5*k*ido) + tr12*ti2 + tr11*ti3; - cr5 = ti11*tr5 + ti12*tr4; - ci5 = ti11*ti5 + ti12*ti4; - cr4 = ti12*tr5 - ti11*tr4; - ci4 = ti12*ti5 - ti11*ti4; - dr3 = cr3 - ci4; - dr4 = cr3 + ci4; - di3 = ci3 + cr4; - di4 = ci3 - cr4; - dr5 = cr2 + ci5; - dr2 = cr2 - ci5; - di5 = ci2 - cr5; - di2 = ci2 + cr5; - ch[i - 1 + (k + l1)*ido] = wa1[i - 2]*dr2 - wa1[i - 1]*di2; - ch[i + (k + l1)*ido] = wa1[i - 2]*di2 + wa1[i - 1]*dr2; - ch[i - 1 + (k + 2*l1)*ido] = wa2[i - 2]*dr3 - wa2[i - 1]*di3; - ch[i + (k + 2*l1)*ido] = wa2[i - 2]*di3 + wa2[i - 1]*dr3; - ch[i - 1 + (k + 3*l1)*ido] = wa3[i - 2]*dr4 - wa3[i - 1]*di4; - ch[i + (k + 3*l1)*ido] = wa3[i - 2]*di4 + wa3[i - 1]*dr4; - ch[i - 1 + (k + 4*l1)*ido] = wa4[i - 2]*dr5 - wa4[i - 1]*di5; - ch[i + (k + 4*l1)*ido] = wa4[i - 2]*di5 + wa4[i - 1]*dr5; - } - } - } /* radb5 */ - - -static void radfg(int ido, int ip, int l1, int idl1, - Treal cc[], Treal ch[], const Treal wa[]) - { - int idij, ipph, i, j, k, l, j2, ic, jc, lc, ik, is, nbd; - Treal dc2, ai1, ai2, ar1, ar2, ds2, dcp, dsp, ar1h, ar2h; - sincos2pi(1, ip, &dsp, &dcp); - ipph = (ip + 1) / 2; - nbd = (ido - 1) / 2; - if (ido != 1) { - for (ik=0; ik<idl1; ik++) ch[ik] = cc[ik]; - for (j=1; j<ip; j++) - for (k=0; k<l1; k++) - ch[(k + j*l1)*ido] = cc[(k + j*l1)*ido]; - if (nbd <= l1) { - is = -ido; - for (j=1; j<ip; j++) { - is += ido; - idij = is-1; - for (i=2; i<ido; i+=2) { - idij += 2; - for (k=0; k<l1; k++) { - ch[i - 1 + (k + j*l1)*ido] = - wa[idij - 1]*cc[i - 1 + (k + j*l1)*ido] + wa[idij]*cc[i + (k + j*l1)*ido]; - ch[i + (k + j*l1)*ido] = - wa[idij - 1]*cc[i + (k + j*l1)*ido] - wa[idij]*cc[i - 1 + (k + j*l1)*ido]; - } - } - } - } else { - is = -ido; - for (j=1; j<ip; j++) { - is += ido; - for (k=0; k<l1; k++) { - idij = is-1; - for (i=2; i<ido; i+=2) { - idij += 2; - ch[i - 1 + (k + j*l1)*ido] = - wa[idij - 1]*cc[i - 1 + (k + j*l1)*ido] + wa[idij]*cc[i + (k + j*l1)*ido]; - ch[i + (k + j*l1)*ido] = - wa[idij - 1]*cc[i + (k + j*l1)*ido] - wa[idij]*cc[i - 1 + (k + j*l1)*ido]; - } - } - } - } - if (nbd >= l1) { - for (j=1; j<ipph; j++) { - jc = ip - j; - for (k=0; k<l1; k++) { - for (i=2; i<ido; i+=2) { - cc[i - 1 + (k + j*l1)*ido] = ch[i - 1 + (k + j*l1)*ido] + ch[i - 1 + (k + jc*l1)*ido]; - cc[i - 1 + (k + jc*l1)*ido] = ch[i + (k + j*l1)*ido] - ch[i + (k + jc*l1)*ido]; - cc[i + (k + j*l1)*ido] = ch[i + (k + j*l1)*ido] + ch[i + (k + jc*l1)*ido]; - cc[i + (k + jc*l1)*ido] = ch[i - 1 + (k + jc*l1)*ido] - ch[i - 1 + (k + j*l1)*ido]; - } - } - } - } else { - for (j=1; j<ipph; j++) { - jc = ip - j; - for (i=2; i<ido; i+=2) { - for (k=0; k<l1; k++) { - cc[i - 1 + (k + j*l1)*ido] = - ch[i - 1 + (k + j*l1)*ido] + ch[i - 1 + (k + jc*l1)*ido]; - cc[i - 1 + (k + jc*l1)*ido] = ch[i + (k + j*l1)*ido] - ch[i + (k + jc*l1)*ido]; - cc[i + (k + j*l1)*ido] = ch[i + (k + j*l1)*ido] + ch[i + (k + jc*l1)*ido]; - cc[i + (k + jc*l1)*ido] = ch[i - 1 + (k + jc*l1)*ido] - ch[i - 1 + (k + j*l1)*ido]; - } - } - } - } - } else { /* now ido == 1 */ - for (ik=0; ik<idl1; ik++) cc[ik] = ch[ik]; - } - for (j=1; j<ipph; j++) { - jc = ip - j; - for (k=0; k<l1; k++) { - cc[(k + j*l1)*ido] = ch[(k + j*l1)*ido] + ch[(k + jc*l1)*ido]; - cc[(k + jc*l1)*ido] = ch[(k + jc*l1)*ido] - ch[(k + j*l1)*ido]; - } - } - - ar1 = 1; - ai1 = 0; - for (l=1; l<ipph; l++) { - lc = ip - l; - ar1h = dcp*ar1 - dsp*ai1; - ai1 = dcp*ai1 + dsp*ar1; - ar1 = ar1h; - for (ik=0; ik<idl1; ik++) { - ch[ik + l*idl1] = cc[ik] + ar1*cc[ik + idl1]; - ch[ik + lc*idl1] = ai1*cc[ik + (ip-1)*idl1]; - } - dc2 = ar1; - ds2 = ai1; - ar2 = ar1; - ai2 = ai1; - for (j=2; j<ipph; j++) { - jc = ip - j; - ar2h = dc2*ar2 - ds2*ai2; - ai2 = dc2*ai2 + ds2*ar2; - ar2 = ar2h; - for (ik=0; ik<idl1; ik++) { - ch[ik + l*idl1] += ar2*cc[ik + j*idl1]; - ch[ik + lc*idl1] += ai2*cc[ik + jc*idl1]; - } - } - } - - for (j=1; j<ipph; j++) - for (ik=0; ik<idl1; ik++) - ch[ik] += cc[ik + j*idl1]; - - if (ido >= l1) { - for (k=0; k<l1; k++) { - for (i=0; i<ido; i++) { - ref(cc,i + k*ip*ido) = ch[i + k*ido]; - } - } - } else { - for (i=0; i<ido; i++) { - for (k=0; k<l1; k++) { - ref(cc,i + k*ip*ido) = ch[i + k*ido]; - } - } - } - for (j=1; j<ipph; j++) { - jc = ip - j; - j2 = 2*j; - for (k=0; k<l1; k++) { - ref(cc,ido-1 + (j2 - 1 + k*ip)*ido) = - ch[(k + j*l1)*ido]; - ref(cc,(j2 + k*ip)*ido) = - ch[(k + jc*l1)*ido]; - } - } - if (ido == 1) return; - if (nbd >= l1) { - for (j=1; j<ipph; j++) { - jc = ip - j; - j2 = 2*j; - for (k=0; k<l1; k++) { - for (i=2; i<ido; i+=2) { - ic = ido - i; - ref(cc,i - 1 + (j2 + k*ip)*ido) = ch[i - 1 + (k + j*l1)*ido] + ch[i - 1 + (k + jc*l1)*ido]; - ref(cc,ic - 1 + (j2 - 1 + k*ip)*ido) = ch[i - 1 + (k + j*l1)*ido] - ch[i - 1 + (k + jc*l1)*ido]; - ref(cc,i + (j2 + k*ip)*ido) = ch[i + (k + j*l1)*ido] + ch[i + (k + jc*l1)*ido]; - ref(cc,ic + (j2 - 1 + k*ip)*ido) = ch[i + (k + jc*l1)*ido] - ch[i + (k + j*l1)*ido]; - } - } - } - } else { - for (j=1; j<ipph; j++) { - jc = ip - j; - j2 = 2*j; - for (i=2; i<ido; i+=2) { - ic = ido - i; - for (k=0; k<l1; k++) { - ref(cc,i - 1 + (j2 + k*ip)*ido) = ch[i - 1 + (k + j*l1)*ido] + ch[i - 1 + (k + jc*l1)*ido]; - ref(cc,ic - 1 + (j2 - 1 + k*ip)*ido) = ch[i - 1 + (k + j*l1)*ido] - ch[i - 1 + (k + jc*l1)*ido]; - ref(cc,i + (j2 + k*ip)*ido) = ch[i + (k + j*l1)*ido] + ch[i + (k + jc*l1)*ido]; - ref(cc,ic + (j2 - 1 + k*ip)*ido) = ch[i + (k + jc*l1)*ido] - ch[i + (k + j*l1)*ido]; - } - } - } - } - } /* radfg */ - - -static void radbg(int ido, int ip, int l1, int idl1, - Treal cc[], Treal ch[], const Treal wa[]) - { - int idij, ipph, i, j, k, l, j2, ic, jc, lc, ik, is; - Treal dc2, ai1, ai2, ar1, ar2, ds2; - int nbd; - Treal dcp, dsp, ar1h, ar2h; - sincos2pi(1, ip, &dsp, &dcp); - nbd = (ido - 1) / 2; - ipph = (ip + 1) / 2; - if (ido >= l1) { - for (k=0; k<l1; k++) { - for (i=0; i<ido; i++) { - ch[i + k*ido] = ref(cc,i + k*ip*ido); - } - } - } else { - for (i=0; i<ido; i++) { - for (k=0; k<l1; k++) { - ch[i + k*ido] = ref(cc,i + k*ip*ido); - } - } - } - for (j=1; j<ipph; j++) { - jc = ip - j; - j2 = 2*j; - for (k=0; k<l1; k++) { - ch[(k + j*l1)*ido] = ref(cc,ido-1 + (j2 - 1 + k*ip)*ido) + ref(cc,ido-1 + (j2 - 1 + k*ip)* - ido); - ch[(k + jc*l1)*ido] = ref(cc,(j2 + k*ip)*ido) + ref(cc,(j2 + k*ip)*ido); - } - } - - if (ido != 1) { - if (nbd >= l1) { - for (j=1; j<ipph; j++) { - jc = ip - j; - for (k=0; k<l1; k++) { - for (i=2; i<ido; i+=2) { - ic = ido - i; - ch[i - 1 + (k + j*l1)*ido] = ref(cc,i - 1 + (2*j + k*ip)*ido) + ref(cc, - ic - 1 + (2*j - 1 + k*ip)*ido); - ch[i - 1 + (k + jc*l1)*ido] = ref(cc,i - 1 + (2*j + k*ip)*ido) - - ref(cc,ic - 1 + (2*j - 1 + k*ip)*ido); - ch[i + (k + j*l1)*ido] = ref(cc,i + (2*j + k*ip)*ido) - ref(cc,ic - + (2*j - 1 + k*ip)*ido); - ch[i + (k + jc*l1)*ido] = ref(cc,i + (2*j + k*ip)*ido) + ref(cc,ic - + (2*j - 1 + k*ip)*ido); - } - } - } - } else { - for (j=1; j<ipph; j++) { - jc = ip - j; - for (i=2; i<ido; i+=2) { - ic = ido - i; - for (k=0; k<l1; k++) { - ch[i - 1 + (k + j*l1)*ido] = ref(cc,i - 1 + (2*j + k*ip)*ido) + ref(cc, - ic - 1 + (2*j - 1 + k*ip)*ido); - ch[i - 1 + (k + jc*l1)*ido] = ref(cc,i - 1 + (2*j + k*ip)*ido) - - ref(cc,ic - 1 + (2*j - 1 + k*ip)*ido); - ch[i + (k + j*l1)*ido] = ref(cc,i + (2*j + k*ip)*ido) - ref(cc,ic - + (2*j - 1 + k*ip)*ido); - ch[i + (k + jc*l1)*ido] = ref(cc,i + (2*j + k*ip)*ido) + ref(cc,ic - + (2*j - 1 + k*ip)*ido); - } - } - } - } - } - - ar1 = 1; - ai1 = 0; - for (l=1; l<ipph; l++) { - lc = ip - l; - ar1h = dcp*ar1 - dsp*ai1; - ai1 = dcp*ai1 + dsp*ar1; - ar1 = ar1h; - for (ik=0; ik<idl1; ik++) { - cc[ik + l*idl1] = ch[ik] + ar1*ch[ik + idl1]; - cc[ik + lc*idl1] = ai1*ch[ik + (ip-1)*idl1]; - } - dc2 = ar1; - ds2 = ai1; - ar2 = ar1; - ai2 = ai1; - for (j=2; j<ipph; j++) { - jc = ip - j; - ar2h = dc2*ar2 - ds2*ai2; - ai2 = dc2*ai2 + ds2*ar2; - ar2 = ar2h; - for (ik=0; ik<idl1; ik++) { - cc[ik + l*idl1] += ar2*ch[ik + j*idl1]; - cc[ik + lc*idl1] += ai2*ch[ik + jc*idl1]; - } - } - } - for (j=1; j<ipph; j++) { - for (ik=0; ik<idl1; ik++) { - ch[ik] += ch[ik + j*idl1]; - } - } - for (j=1; j<ipph; j++) { - jc = ip - j; - for (k=0; k<l1; k++) { - ch[(k + j*l1)*ido] = cc[(k + j*l1)*ido] - cc[(k + jc*l1)*ido]; - ch[(k + jc*l1)*ido] = cc[(k + j*l1)*ido] + cc[(k + jc*l1)*ido]; - } - } - - if (ido == 1) return; - if (nbd >= l1) { - for (j=1; j<ipph; j++) { - jc = ip - j; - for (k=0; k<l1; k++) { - for (i=2; i<ido; i+=2) { - ch[i - 1 + (k + j*l1)*ido] = cc[i - 1 + (k + j*l1)*ido] - cc[i + (k + jc*l1)*ido]; - ch[i - 1 + (k + jc*l1)*ido] = cc[i - 1 + (k + j*l1)*ido] + cc[i + (k + jc*l1)*ido]; - ch[i + (k + j*l1)*ido] = cc[i + (k + j*l1)*ido] + cc[i - 1 + (k + jc*l1)*ido]; - ch[i + (k + jc*l1)*ido] = cc[i + (k + j*l1)*ido] - cc[i - 1 + (k + jc*l1)*ido]; - } - } - } - } else { - for (j=1; j<ipph; j++) { - jc = ip - j; - for (i=2; i<ido; i+=2) { - for (k=0; k<l1; k++) { - ch[i - 1 + (k + j*l1)*ido] = cc[i - 1 + (k + j*l1)*ido] - cc[i + (k + jc*l1)*ido]; - ch[i - 1 + (k + jc*l1)*ido] = cc[i - 1 + (k + j *l1)*ido] + cc[i + (k + jc*l1)*ido]; - ch[i + (k + j*l1)*ido] = cc[i + (k + j*l1)*ido] + cc[i - 1 + (k + jc*l1)*ido]; - ch[i + (k + jc*l1)*ido] = cc[i + (k + j*l1)*ido] - cc[i - 1 + (k + jc*l1)*ido]; - } - } - } - } - for (ik=0; ik<idl1; ik++) cc[ik] = ch[ik]; - for (j=1; j<ip; j++) - for (k=0; k<l1; k++) - cc[(k + j*l1)*ido] = ch[(k + j*l1)*ido]; - if (nbd <= l1) { - is = -ido; - for (j=1; j<ip; j++) { - is += ido; - idij = is-1; - for (i=2; i<ido; i+=2) { - idij += 2; - for (k=0; k<l1; k++) { - cc[i - 1 + (k + j*l1)*ido] = wa[idij - 1]*ch[i - 1 + (k + j*l1)*ido] - wa[idij]* - ch[i + (k + j*l1)*ido]; - cc[i + (k + j*l1)*ido] = wa[idij - 1]*ch[i + (k + j*l1)*ido] + wa[idij]*ch[i - 1 + (k + j*l1)*ido]; - } - } - } - } else { - is = -ido; - for (j=1; j<ip; j++) { - is += ido; - for (k=0; k<l1; k++) { - idij = is - 1; - for (i=2; i<ido; i+=2) { - idij += 2; - cc[i - 1 + (k + j*l1)*ido] = wa[idij-1]*ch[i - 1 + (k + j*l1)*ido] - wa[idij]* - ch[i + (k + j*l1)*ido]; - cc[i + (k + j*l1)*ido] = wa[idij-1]*ch[i + (k + j*l1)*ido] + wa[idij]*ch[i - 1 + (k + j*l1)*ido]; - } - } - } - } - } /* radbg */ - - /* ------------------------------------------------------------ -cfftf1, npy_cfftf, npy_cfftb, cffti1, npy_cffti. Complex FFTs. ---------------------------------------------------------------- */ - -static void cfftf1(int n, Treal c[], Treal ch[], const Treal wa[], const int ifac[MAXFAC+2], int isign) - { - int idot, i; - int k1, l1, l2; - int na, nf, ip, iw, ix2, ix3, ix4, nac, ido, idl1; - Treal *cinput, *coutput; - nf = ifac[1]; - na = 0; - l1 = 1; - iw = 0; - for (k1=2; k1<=nf+1; k1++) { - ip = ifac[k1]; - l2 = ip*l1; - ido = n / l2; - idot = ido + ido; - idl1 = idot*l1; - if (na) { - cinput = ch; - coutput = c; - } else { - cinput = c; - coutput = ch; - } - switch (ip) { - case 4: - ix2 = iw + idot; - ix3 = ix2 + idot; - passf4(idot, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3], isign); - na = !na; - break; - case 2: - passf2(idot, l1, cinput, coutput, &wa[iw], isign); - na = !na; - break; - case 3: - ix2 = iw + idot; - passf3(idot, l1, cinput, coutput, &wa[iw], &wa[ix2], isign); - na = !na; - break; - case 5: - ix2 = iw + idot; - ix3 = ix2 + idot; - ix4 = ix3 + idot; - passf5(idot, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4], isign); - na = !na; - break; - default: - passf(&nac, idot, ip, l1, idl1, cinput, coutput, &wa[iw], isign); - if (nac != 0) na = !na; - } - l1 = l2; - iw += (ip - 1)*idot; - } - if (na == 0) return; - for (i=0; i<2*n; i++) c[i] = ch[i]; - } /* cfftf1 */ - - -NPY_VISIBILITY_HIDDEN void npy_cfftf(int n, Treal c[], Treal wsave[]) - { - int iw1, iw2; - if (n == 1) return; - iw1 = 2*n; - iw2 = iw1 + 2*n; - cfftf1(n, c, wsave, wsave+iw1, (int*)(wsave+iw2), -1); - } /* npy_cfftf */ - - -NPY_VISIBILITY_HIDDEN void npy_cfftb(int n, Treal c[], Treal wsave[]) - { - int iw1, iw2; - if (n == 1) return; - iw1 = 2*n; - iw2 = iw1 + 2*n; - cfftf1(n, c, wsave, wsave+iw1, (int*)(wsave+iw2), +1); - } /* npy_cfftb */ - - -static void factorize(int n, int ifac[MAXFAC+2], const int ntryh[NSPECIAL]) - /* Factorize n in factors in ntryh and rest. On exit, -ifac[0] contains n and ifac[1] contains number of factors, -the factors start from ifac[2]. */ - { - int ntry=3, i, j=0, ib, nf=0, nl=n, nq, nr; -startloop: - if (j < NSPECIAL) - ntry = ntryh[j]; - else - ntry+= 2; - j++; - do { - nq = nl / ntry; - nr = nl - ntry*nq; - if (nr != 0) goto startloop; - nf++; - ifac[nf + 1] = ntry; - nl = nq; - if (ntry == 2 && nf != 1) { - for (i=2; i<=nf; i++) { - ib = nf - i + 2; - ifac[ib + 1] = ifac[ib]; - } - ifac[2] = 2; - } - } while (nl != 1); - ifac[0] = n; - ifac[1] = nf; - } - - -static void cffti1(int n, Treal wa[], int ifac[MAXFAC+2]) - { - int fi, idot, i, j; - int i1, k1, l1, l2; - int ld, ii, nf, ip; - int ido, ipm; - - static const int ntryh[NSPECIAL] = { - 3,4,2,5 }; /* Do not change the order of these. */ - - factorize(n,ifac,ntryh); - nf = ifac[1]; - i = 1; - l1 = 1; - for (k1=1; k1<=nf; k1++) { - ip = ifac[k1+1]; - ld = 0; - l2 = l1*ip; - ido = n / l2; - idot = ido + ido + 2; - ipm = ip - 1; - for (j=1; j<=ipm; j++) { - i1 = i; - wa[i-1] = 1; - wa[i] = 0; - ld += l1; - fi = 0; - for (ii=4; ii<=idot; ii+=2) { - i+= 2; - fi+= 1; - sincos2pi(fi*ld, n, wa+i, wa+i-1); - } - if (ip > 5) { - wa[i1-1] = wa[i-1]; - wa[i1] = wa[i]; - } - } - l1 = l2; - } - } /* cffti1 */ - - -NPY_VISIBILITY_HIDDEN void npy_cffti(int n, Treal wsave[]) - { - int iw1, iw2; - if (n == 1) return; - iw1 = 2*n; - iw2 = iw1 + 2*n; - cffti1(n, wsave+iw1, (int*)(wsave+iw2)); - } /* npy_cffti */ - - /* ------------------------------------------------------------------- -rfftf1, rfftb1, npy_rfftf, npy_rfftb, rffti1, npy_rffti. Treal FFTs. ----------------------------------------------------------------------- */ - -static void rfftf1(int n, Treal c[], Treal ch[], const Treal wa[], const int ifac[MAXFAC+2]) - { - int i; - int k1, l1, l2, na, kh, nf, ip, iw, ix2, ix3, ix4, ido, idl1; - Treal *cinput, *coutput; - nf = ifac[1]; - na = 1; - l2 = n; - iw = n-1; - for (k1 = 1; k1 <= nf; ++k1) { - kh = nf - k1; - ip = ifac[kh + 2]; - l1 = l2 / ip; - ido = n / l2; - idl1 = ido*l1; - iw -= (ip - 1)*ido; - na = !na; - if (na) { - cinput = ch; - coutput = c; - } else { - cinput = c; - coutput = ch; - } - switch (ip) { - case 4: - ix2 = iw + ido; - ix3 = ix2 + ido; - radf4(ido, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3]); - break; - case 2: - radf2(ido, l1, cinput, coutput, &wa[iw]); - break; - case 3: - ix2 = iw + ido; - radf3(ido, l1, cinput, coutput, &wa[iw], &wa[ix2]); - break; - case 5: - ix2 = iw + ido; - ix3 = ix2 + ido; - ix4 = ix3 + ido; - radf5(ido, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4]); - break; - default: - if (ido == 1) - na = !na; - if (na == 0) { - radfg(ido, ip, l1, idl1, c, ch, &wa[iw]); - na = 1; - } else { - radfg(ido, ip, l1, idl1, ch, c, &wa[iw]); - na = 0; - } - } - l2 = l1; - } - if (na == 1) return; - for (i = 0; i < n; i++) c[i] = ch[i]; - } /* rfftf1 */ - - -static void rfftb1(int n, Treal c[], Treal ch[], const Treal wa[], const int ifac[MAXFAC+2]) - { - int i; - int k1, l1, l2, na, nf, ip, iw, ix2, ix3, ix4, ido, idl1; - Treal *cinput, *coutput; - nf = ifac[1]; - na = 0; - l1 = 1; - iw = 0; - for (k1=1; k1<=nf; k1++) { - ip = ifac[k1 + 1]; - l2 = ip*l1; - ido = n / l2; - idl1 = ido*l1; - if (na) { - cinput = ch; - coutput = c; - } else { - cinput = c; - coutput = ch; - } - switch (ip) { - case 4: - ix2 = iw + ido; - ix3 = ix2 + ido; - radb4(ido, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3]); - na = !na; - break; - case 2: - radb2(ido, l1, cinput, coutput, &wa[iw]); - na = !na; - break; - case 3: - ix2 = iw + ido; - radb3(ido, l1, cinput, coutput, &wa[iw], &wa[ix2]); - na = !na; - break; - case 5: - ix2 = iw + ido; - ix3 = ix2 + ido; - ix4 = ix3 + ido; - radb5(ido, l1, cinput, coutput, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4]); - na = !na; - break; - default: - radbg(ido, ip, l1, idl1, cinput, coutput, &wa[iw]); - if (ido == 1) na = !na; - } - l1 = l2; - iw += (ip - 1)*ido; - } - if (na == 0) return; - for (i=0; i<n; i++) c[i] = ch[i]; - } /* rfftb1 */ - - -NPY_VISIBILITY_HIDDEN void npy_rfftf(int n, Treal r[], Treal wsave[]) - { - if (n == 1) return; - rfftf1(n, r, wsave, wsave+n, (int*)(wsave+2*n)); - } /* npy_rfftf */ - - -NPY_VISIBILITY_HIDDEN void npy_rfftb(int n, Treal r[], Treal wsave[]) - { - if (n == 1) return; - rfftb1(n, r, wsave, wsave+n, (int*)(wsave+2*n)); - } /* npy_rfftb */ - - -static void rffti1(int n, Treal wa[], int ifac[MAXFAC+2]) - { - int fi, i, j; - int k1, l1, l2; - int ld, ii, nf, ip, is; - int ido, ipm, nfm1; - static const int ntryh[NSPECIAL] = { - 4,2,3,5 }; /* Do not change the order of these. */ - factorize(n,ifac,ntryh); - nf = ifac[1]; - is = 0; - nfm1 = nf - 1; - l1 = 1; - if (nfm1 == 0) return; - for (k1 = 1; k1 <= nfm1; k1++) { - ip = ifac[k1 + 1]; - ld = 0; - l2 = l1*ip; - ido = n / l2; - ipm = ip - 1; - for (j = 1; j <= ipm; ++j) { - ld += l1; - i = is; - fi = 0; - for (ii = 3; ii <= ido; ii += 2) { - i += 2; - fi += 1; - sincos2pi(fi*ld, n, wa+i-1, wa+i-2); - } - is += ido; - } - l1 = l2; - } - } /* rffti1 */ - - -NPY_VISIBILITY_HIDDEN void npy_rffti(int n, Treal wsave[]) - { - if (n == 1) return; - rffti1(n, wsave+n, (int*)(wsave+2*n)); - } /* npy_rffti */ - -#ifdef __cplusplus -} -#endif diff --git a/numpy/fft/fftpack.h b/numpy/fft/fftpack.h deleted file mode 100644 index 5e8f4631c..000000000 --- a/numpy/fft/fftpack.h +++ /dev/null @@ -1,28 +0,0 @@ -/* - * This file is part of tela the Tensor Language. - * Copyright (c) 1994-1995 Pekka Janhunen - */ - -#ifdef __cplusplus -extern "C" { -#endif - -#define DOUBLE - -#ifdef DOUBLE -#define Treal double -#else -#define Treal float -#endif - -extern NPY_VISIBILITY_HIDDEN void npy_cfftf(int N, Treal data[], const Treal wrk[]); -extern NPY_VISIBILITY_HIDDEN void npy_cfftb(int N, Treal data[], const Treal wrk[]); -extern NPY_VISIBILITY_HIDDEN void npy_cffti(int N, Treal wrk[]); - -extern NPY_VISIBILITY_HIDDEN void npy_rfftf(int N, Treal data[], const Treal wrk[]); -extern NPY_VISIBILITY_HIDDEN void npy_rfftb(int N, Treal data[], const Treal wrk[]); -extern NPY_VISIBILITY_HIDDEN void npy_rffti(int N, Treal wrk[]); - -#ifdef __cplusplus -} -#endif diff --git a/numpy/fft/fftpack_litemodule.c b/numpy/fft/fftpack_litemodule.c deleted file mode 100644 index bd6cfc120..000000000 --- a/numpy/fft/fftpack_litemodule.c +++ /dev/null @@ -1,366 +0,0 @@ -#define NPY_NO_DEPRECATED_API NPY_API_VERSION - -#include "Python.h" -#include "numpy/arrayobject.h" -#include "fftpack.h" - -static PyObject *ErrorObject; - -static const char fftpack_cfftf__doc__[] = ""; - -static PyObject * -fftpack_cfftf(PyObject *NPY_UNUSED(self), PyObject *args) -{ - PyObject *op1, *op2; - PyArrayObject *data; - PyArray_Descr *descr; - double *wsave, *dptr; - npy_intp nsave; - int npts, nrepeats, i; - - if(!PyArg_ParseTuple(args, "OO:cfftf", &op1, &op2)) { - return NULL; - } - data = (PyArrayObject *)PyArray_CopyFromObject(op1, - NPY_CDOUBLE, 1, 0); - if (data == NULL) { - return NULL; - } - descr = PyArray_DescrFromType(NPY_DOUBLE); - if (PyArray_AsCArray(&op2, (void *)&wsave, &nsave, 1, descr) == -1) { - goto fail; - } - if (data == NULL) { - goto fail; - } - - npts = PyArray_DIM(data, PyArray_NDIM(data) - 1); - if (nsave != npts*4 + 15) { - PyErr_SetString(ErrorObject, "invalid work array for fft size"); - goto fail; - } - - nrepeats = PyArray_SIZE(data)/npts; - dptr = (double *)PyArray_DATA(data); - Py_BEGIN_ALLOW_THREADS; - NPY_SIGINT_ON; - for (i = 0; i < nrepeats; i++) { - npy_cfftf(npts, dptr, wsave); - dptr += npts*2; - } - NPY_SIGINT_OFF; - Py_END_ALLOW_THREADS; - PyArray_Free(op2, (char *)wsave); - return (PyObject *)data; - -fail: - PyArray_Free(op2, (char *)wsave); - Py_DECREF(data); - return NULL; -} - -static const char fftpack_cfftb__doc__[] = ""; - -static PyObject * -fftpack_cfftb(PyObject *NPY_UNUSED(self), PyObject *args) -{ - PyObject *op1, *op2; - PyArrayObject *data; - PyArray_Descr *descr; - double *wsave, *dptr; - npy_intp nsave; - int npts, nrepeats, i; - - if(!PyArg_ParseTuple(args, "OO:cfftb", &op1, &op2)) { - return NULL; - } - data = (PyArrayObject *)PyArray_CopyFromObject(op1, - NPY_CDOUBLE, 1, 0); - if (data == NULL) { - return NULL; - } - descr = PyArray_DescrFromType(NPY_DOUBLE); - if (PyArray_AsCArray(&op2, (void *)&wsave, &nsave, 1, descr) == -1) { - goto fail; - } - if (data == NULL) { - goto fail; - } - - npts = PyArray_DIM(data, PyArray_NDIM(data) - 1); - if (nsave != npts*4 + 15) { - PyErr_SetString(ErrorObject, "invalid work array for fft size"); - goto fail; - } - - nrepeats = PyArray_SIZE(data)/npts; - dptr = (double *)PyArray_DATA(data); - Py_BEGIN_ALLOW_THREADS; - NPY_SIGINT_ON; - for (i = 0; i < nrepeats; i++) { - npy_cfftb(npts, dptr, wsave); - dptr += npts*2; - } - NPY_SIGINT_OFF; - Py_END_ALLOW_THREADS; - PyArray_Free(op2, (char *)wsave); - return (PyObject *)data; - -fail: - PyArray_Free(op2, (char *)wsave); - Py_DECREF(data); - return NULL; -} - -static const char fftpack_cffti__doc__[] = ""; - -static PyObject * -fftpack_cffti(PyObject *NPY_UNUSED(self), PyObject *args) -{ - PyArrayObject *op; - npy_intp dim; - long n; - - if (!PyArg_ParseTuple(args, "l:cffti", &n)) { - return NULL; - } - /*Magic size needed by npy_cffti*/ - dim = 4*n + 15; - /*Create a 1 dimensional array of dimensions of type double*/ - op = (PyArrayObject *)PyArray_SimpleNew(1, &dim, NPY_DOUBLE); - if (op == NULL) { - return NULL; - } - - Py_BEGIN_ALLOW_THREADS; - NPY_SIGINT_ON; - npy_cffti(n, (double *)PyArray_DATA((PyArrayObject*)op)); - NPY_SIGINT_OFF; - Py_END_ALLOW_THREADS; - - return (PyObject *)op; -} - -static const char fftpack_rfftf__doc__[] = ""; - -static PyObject * -fftpack_rfftf(PyObject *NPY_UNUSED(self), PyObject *args) -{ - PyObject *op1, *op2; - PyArrayObject *data, *ret; - PyArray_Descr *descr; - double *wsave = NULL, *dptr, *rptr; - npy_intp nsave; - int npts, nrepeats, i, rstep; - - if(!PyArg_ParseTuple(args, "OO:rfftf", &op1, &op2)) { - return NULL; - } - data = (PyArrayObject *)PyArray_ContiguousFromObject(op1, - NPY_DOUBLE, 1, 0); - if (data == NULL) { - return NULL; - } - /* FIXME, direct access changing contents of data->dimensions */ - npts = PyArray_DIM(data, PyArray_NDIM(data) - 1); - PyArray_DIMS(data)[PyArray_NDIM(data) - 1] = npts/2 + 1; - ret = (PyArrayObject *)PyArray_Zeros(PyArray_NDIM(data), - PyArray_DIMS(data), PyArray_DescrFromType(NPY_CDOUBLE), 0); - if (ret == NULL) { - goto fail; - } - PyArray_DIMS(data)[PyArray_NDIM(data) - 1] = npts; - rstep = PyArray_DIM(ret, PyArray_NDIM(ret) - 1)*2; - - descr = PyArray_DescrFromType(NPY_DOUBLE); - if (PyArray_AsCArray(&op2, (void *)&wsave, &nsave, 1, descr) == -1) { - goto fail; - } - if (data == NULL || ret == NULL) { - goto fail; - } - if (nsave != npts*2+15) { - PyErr_SetString(ErrorObject, "invalid work array for fft size"); - goto fail; - } - - nrepeats = PyArray_SIZE(data)/npts; - rptr = (double *)PyArray_DATA(ret); - dptr = (double *)PyArray_DATA(data); - - Py_BEGIN_ALLOW_THREADS; - NPY_SIGINT_ON; - for (i = 0; i < nrepeats; i++) { - memcpy((char *)(rptr+1), dptr, npts*sizeof(double)); - npy_rfftf(npts, rptr+1, wsave); - rptr[0] = rptr[1]; - rptr[1] = 0.0; - rptr += rstep; - dptr += npts; - } - NPY_SIGINT_OFF; - Py_END_ALLOW_THREADS; - PyArray_Free(op2, (char *)wsave); - Py_DECREF(data); - return (PyObject *)ret; - -fail: - PyArray_Free(op2, (char *)wsave); - Py_XDECREF(data); - Py_XDECREF(ret); - return NULL; -} - -static const char fftpack_rfftb__doc__[] = ""; - -static PyObject * -fftpack_rfftb(PyObject *NPY_UNUSED(self), PyObject *args) -{ - PyObject *op1, *op2; - PyArrayObject *data, *ret; - PyArray_Descr *descr; - double *wsave, *dptr, *rptr; - npy_intp nsave; - int npts, nrepeats, i; - - if(!PyArg_ParseTuple(args, "OO:rfftb", &op1, &op2)) { - return NULL; - } - data = (PyArrayObject *)PyArray_ContiguousFromObject(op1, - NPY_CDOUBLE, 1, 0); - if (data == NULL) { - return NULL; - } - npts = PyArray_DIM(data, PyArray_NDIM(data) - 1); - ret = (PyArrayObject *)PyArray_Zeros(PyArray_NDIM(data), PyArray_DIMS(data), - PyArray_DescrFromType(NPY_DOUBLE), 0); - - descr = PyArray_DescrFromType(NPY_DOUBLE); - if (PyArray_AsCArray(&op2, (void *)&wsave, &nsave, 1, descr) == -1) { - goto fail; - } - if (data == NULL || ret == NULL) { - goto fail; - } - if (nsave != npts*2 + 15) { - PyErr_SetString(ErrorObject, "invalid work array for fft size"); - goto fail; - } - - nrepeats = PyArray_SIZE(ret)/npts; - rptr = (double *)PyArray_DATA(ret); - dptr = (double *)PyArray_DATA(data); - - Py_BEGIN_ALLOW_THREADS; - NPY_SIGINT_ON; - for (i = 0; i < nrepeats; i++) { - memcpy((char *)(rptr + 1), (dptr + 2), (npts - 1)*sizeof(double)); - rptr[0] = dptr[0]; - npy_rfftb(npts, rptr, wsave); - rptr += npts; - dptr += npts*2; - } - NPY_SIGINT_OFF; - Py_END_ALLOW_THREADS; - PyArray_Free(op2, (char *)wsave); - Py_DECREF(data); - return (PyObject *)ret; - -fail: - PyArray_Free(op2, (char *)wsave); - Py_XDECREF(data); - Py_XDECREF(ret); - return NULL; -} - -static const char fftpack_rffti__doc__[] = ""; - -static PyObject * -fftpack_rffti(PyObject *NPY_UNUSED(self), PyObject *args) -{ - PyArrayObject *op; - npy_intp dim; - long n; - - if (!PyArg_ParseTuple(args, "l:rffti", &n)) { - return NULL; - } - /*Magic size needed by npy_rffti*/ - dim = 2*n + 15; - /*Create a 1 dimensional array of dimensions of type double*/ - op = (PyArrayObject *)PyArray_SimpleNew(1, &dim, NPY_DOUBLE); - if (op == NULL) { - return NULL; - } - Py_BEGIN_ALLOW_THREADS; - NPY_SIGINT_ON; - npy_rffti(n, (double *)PyArray_DATA((PyArrayObject*)op)); - NPY_SIGINT_OFF; - Py_END_ALLOW_THREADS; - - return (PyObject *)op; -} - - -/* List of methods defined in the module */ - -static struct PyMethodDef fftpack_methods[] = { - {"cfftf", fftpack_cfftf, 1, fftpack_cfftf__doc__}, - {"cfftb", fftpack_cfftb, 1, fftpack_cfftb__doc__}, - {"cffti", fftpack_cffti, 1, fftpack_cffti__doc__}, - {"rfftf", fftpack_rfftf, 1, fftpack_rfftf__doc__}, - {"rfftb", fftpack_rfftb, 1, fftpack_rfftb__doc__}, - {"rffti", fftpack_rffti, 1, fftpack_rffti__doc__}, - {NULL, NULL, 0, NULL} /* sentinel */ -}; - -#if PY_MAJOR_VERSION >= 3 -static struct PyModuleDef moduledef = { - PyModuleDef_HEAD_INIT, - "fftpack_lite", - NULL, - -1, - fftpack_methods, - NULL, - NULL, - NULL, - NULL -}; -#endif - -/* Initialization function for the module */ -#if PY_MAJOR_VERSION >= 3 -#define RETVAL(x) x -PyMODINIT_FUNC PyInit_fftpack_lite(void) -#else -#define RETVAL(x) -PyMODINIT_FUNC -initfftpack_lite(void) -#endif -{ - PyObject *m,*d; -#if PY_MAJOR_VERSION >= 3 - m = PyModule_Create(&moduledef); -#else - static const char fftpack_module_documentation[] = ""; - - m = Py_InitModule4("fftpack_lite", fftpack_methods, - fftpack_module_documentation, - (PyObject*)NULL,PYTHON_API_VERSION); -#endif - if (m == NULL) { - return RETVAL(NULL); - } - - /* Import the array object */ - import_array(); - - /* Add some symbolic constants to the module */ - d = PyModule_GetDict(m); - ErrorObject = PyErr_NewException("fftpack.error", NULL, NULL); - PyDict_SetItemString(d, "error", ErrorObject); - - /* XXXX Add constants here */ - - return RETVAL(m); -} diff --git a/numpy/fft/helper.py b/numpy/fft/helper.py index 9985f6d4c..a920a4ac0 100644 --- a/numpy/fft/helper.py +++ b/numpy/fft/helper.py @@ -4,11 +4,6 @@ Discrete Fourier Transforms - helper.py """ from __future__ import division, absolute_import, print_function -import collections -try: - import threading -except ImportError: - import dummy_threading as threading from numpy.compat import integer_types from numpy.core import integer, empty, arange, asarray, roll from numpy.core.overrides import array_function_dispatch, set_module @@ -227,99 +222,3 @@ def rfftfreq(n, d=1.0): N = n//2 + 1 results = arange(0, N, dtype=int) return results * val - - -class _FFTCache(object): - """ - Cache for the FFT twiddle factors as an LRU (least recently used) cache. - - Parameters - ---------- - max_size_in_mb : int - Maximum memory usage of the cache before items are being evicted. - max_item_count : int - Maximum item count of the cache before items are being evicted. - - Notes - ----- - Items will be evicted if either limit has been reached upon getting and - setting. The maximum memory usages is not strictly the given - ``max_size_in_mb`` but rather - ``max(max_size_in_mb, 1.5 * size_of_largest_item)``. Thus the cache will - never be completely cleared - at least one item will remain and a single - large item can cause the cache to retain several smaller items even if the - given maximum cache size has been exceeded. - """ - def __init__(self, max_size_in_mb, max_item_count): - self._max_size_in_bytes = max_size_in_mb * 1024 ** 2 - self._max_item_count = max_item_count - self._dict = collections.OrderedDict() - self._lock = threading.Lock() - - def put_twiddle_factors(self, n, factors): - """ - Store twiddle factors for an FFT of length n in the cache. - - Putting multiple twiddle factors for a certain n will store it multiple - times. - - Parameters - ---------- - n : int - Data length for the FFT. - factors : ndarray - The actual twiddle values. - """ - with self._lock: - # Pop + later add to move it to the end for LRU behavior. - # Internally everything is stored in a dictionary whose values are - # lists. - try: - value = self._dict.pop(n) - except KeyError: - value = [] - value.append(factors) - self._dict[n] = value - self._prune_cache() - - def pop_twiddle_factors(self, n): - """ - Pop twiddle factors for an FFT of length n from the cache. - - Will return None if the requested twiddle factors are not available in - the cache. - - Parameters - ---------- - n : int - Data length for the FFT. - - Returns - ------- - out : ndarray or None - The retrieved twiddle factors if available, else None. - """ - with self._lock: - if n not in self._dict or not self._dict[n]: - return None - # Pop + later add to move it to the end for LRU behavior. - all_values = self._dict.pop(n) - value = all_values.pop() - # Only put pack if there are still some arrays left in the list. - if all_values: - self._dict[n] = all_values - return value - - def _prune_cache(self): - # Always keep at least one item. - while len(self._dict) > 1 and ( - len(self._dict) > self._max_item_count or self._check_size()): - self._dict.popitem(last=False) - - def _check_size(self): - item_sizes = [sum(_j.nbytes for _j in _i) - for _i in self._dict.values() if _i] - if not item_sizes: - return False - max_size = max(self._max_size_in_bytes, 1.5 * max(item_sizes)) - return sum(item_sizes) > max_size diff --git a/numpy/fft/pocketfft.c b/numpy/fft/pocketfft.c new file mode 100644 index 000000000..10a741b6f --- /dev/null +++ b/numpy/fft/pocketfft.c @@ -0,0 +1,2398 @@ +/* + * This file is part of pocketfft. + * Licensed under a 3-clause BSD style license - see LICENSE.md + */ + +/* + * Main implementation file. + * + * Copyright (C) 2004-2018 Max-Planck-Society + * \author Martin Reinecke + */ + +#include <math.h> +#include <string.h> +#include <stdlib.h> + +#include "npy_config.h" +#define restrict NPY_RESTRICT + +#define RALLOC(type,num) \ + ((type *)malloc((num)*sizeof(type))) +#define DEALLOC(ptr) \ + do { free(ptr); (ptr)=NULL; } while(0) + +#define SWAP(a,b,type) \ + do { type tmp_=(a); (a)=(b); (b)=tmp_; } while(0) + +#ifdef __GNUC__ +#define NOINLINE __attribute__((noinline)) +#define WARN_UNUSED_RESULT __attribute__ ((warn_unused_result)) +#else +#define NOINLINE +#define WARN_UNUSED_RESULT +#endif + +struct cfft_plan_i; +typedef struct cfft_plan_i * cfft_plan; +struct rfft_plan_i; +typedef struct rfft_plan_i * rfft_plan; + +// adapted from https://stackoverflow.com/questions/42792939/ +// CAUTION: this function only works for arguments in the range [-0.25; 0.25]! +static void my_sincosm1pi (double a, double *restrict res) + { + double s = a * a; + /* Approximate cos(pi*x)-1 for x in [-0.25,0.25] */ + double r = -1.0369917389758117e-4; + r = fma (r, s, 1.9294935641298806e-3); + r = fma (r, s, -2.5806887942825395e-2); + r = fma (r, s, 2.3533063028328211e-1); + r = fma (r, s, -1.3352627688538006e+0); + r = fma (r, s, 4.0587121264167623e+0); + r = fma (r, s, -4.9348022005446790e+0); + double c = r*s; + /* Approximate sin(pi*x) for x in [-0.25,0.25] */ + r = 4.6151442520157035e-4; + r = fma (r, s, -7.3700183130883555e-3); + r = fma (r, s, 8.2145868949323936e-2); + r = fma (r, s, -5.9926452893214921e-1); + r = fma (r, s, 2.5501640398732688e+0); + r = fma (r, s, -5.1677127800499516e+0); + s = s * a; + r = r * s; + s = fma (a, 3.1415926535897931e+0, r); + res[0] = c; + res[1] = s; + } + +NOINLINE static void calc_first_octant(size_t den, double * restrict res) + { + size_t n = (den+4)>>3; + if (n==0) return; + res[0]=1.; res[1]=0.; + if (n==1) return; + size_t l1=(size_t)sqrt(n); + for (size_t i=1; i<l1; ++i) + my_sincosm1pi((2.*i)/den,&res[2*i]); + size_t start=l1; + while(start<n) + { + double cs[2]; + my_sincosm1pi((2.*start)/den,cs); + res[2*start] = cs[0]+1.; + res[2*start+1] = cs[1]; + size_t end = l1; + if (start+end>n) end = n-start; + for (size_t i=1; i<end; ++i) + { + double csx[2]={res[2*i], res[2*i+1]}; + res[2*(start+i)] = ((cs[0]*csx[0] - cs[1]*csx[1] + cs[0]) + csx[0]) + 1.; + res[2*(start+i)+1] = (cs[0]*csx[1] + cs[1]*csx[0]) + cs[1] + csx[1]; + } + start += l1; + } + for (size_t i=1; i<l1; ++i) + res[2*i] += 1.; + } + +NOINLINE static void calc_first_quadrant(size_t n, double * restrict res) + { + double * restrict p = res+n; + calc_first_octant(n<<1, p); + size_t ndone=(n+2)>>2; + size_t i=0, idx1=0, idx2=2*ndone-2; + for (; i+1<ndone; i+=2, idx1+=2, idx2-=2) + { + res[idx1] = p[2*i]; + res[idx1+1] = p[2*i+1]; + res[idx2] = p[2*i+3]; + res[idx2+1] = p[2*i+2]; + } + if (i!=ndone) + { + res[idx1 ] = p[2*i]; + res[idx1+1] = p[2*i+1]; + } + } + +NOINLINE static void calc_first_half(size_t n, double * restrict res) + { + int ndone=(n+1)>>1; + double * p = res+n-1; + calc_first_octant(n<<2, p); + int i4=0, in=n, i=0; + for (; i4<=in-i4; ++i, i4+=4) // octant 0 + { + res[2*i] = p[2*i4]; res[2*i+1] = p[2*i4+1]; + } + for (; i4-in <= 0; ++i, i4+=4) // octant 1 + { + int xm = in-i4; + res[2*i] = p[2*xm+1]; res[2*i+1] = p[2*xm]; + } + for (; i4<=3*in-i4; ++i, i4+=4) // octant 2 + { + int xm = i4-in; + res[2*i] = -p[2*xm+1]; res[2*i+1] = p[2*xm]; + } + for (; i<ndone; ++i, i4+=4) // octant 3 + { + int xm = 2*in-i4; + res[2*i] = -p[2*xm]; res[2*i+1] = p[2*xm+1]; + } + } + +NOINLINE static void fill_first_quadrant(size_t n, double * restrict res) + { + const double hsqt2 = 0.707106781186547524400844362104849; + size_t quart = n>>2; + if ((n&7)==0) + res[quart] = res[quart+1] = hsqt2; + for (size_t i=2, j=2*quart-2; i<quart; i+=2, j-=2) + { + res[j ] = res[i+1]; + res[j+1] = res[i ]; + } + } + +NOINLINE static void fill_first_half(size_t n, double * restrict res) + { + size_t half = n>>1; + if ((n&3)==0) + for (size_t i=0; i<half; i+=2) + { + res[i+half] = -res[i+1]; + res[i+half+1] = res[i ]; + } + else + for (size_t i=2, j=2*half-2; i<half; i+=2, j-=2) + { + res[j ] = -res[i ]; + res[j+1] = res[i+1]; + } + } + +NOINLINE static void fill_second_half(size_t n, double * restrict res) + { + if ((n&1)==0) + for (size_t i=0; i<n; ++i) + res[i+n] = -res[i]; + else + for (size_t i=2, j=2*n-2; i<n; i+=2, j-=2) + { + res[j ] = res[i ]; + res[j+1] = -res[i+1]; + } + } + +NOINLINE static void sincos_2pibyn_half(size_t n, double * restrict res) + { + if ((n&3)==0) + { + calc_first_octant(n, res); + fill_first_quadrant(n, res); + fill_first_half(n, res); + } + else if ((n&1)==0) + { + calc_first_quadrant(n, res); + fill_first_half(n, res); + } + else + calc_first_half(n, res); + } + +NOINLINE static void sincos_2pibyn(size_t n, double * restrict res) + { + sincos_2pibyn_half(n, res); + fill_second_half(n, res); + } + +NOINLINE static size_t largest_prime_factor (size_t n) + { + size_t res=1; + size_t tmp; + while (((tmp=(n>>1))<<1)==n) + { res=2; n=tmp; } + + size_t limit=(size_t)sqrt(n+0.01); + for (size_t x=3; x<=limit; x+=2) + while (((tmp=(n/x))*x)==n) + { + res=x; + n=tmp; + limit=(size_t)sqrt(n+0.01); + } + if (n>1) res=n; + + return res; + } + +NOINLINE static double cost_guess (size_t n) + { + const double lfp=1.1; // penalty for non-hardcoded larger factors + size_t ni=n; + double result=0.; + size_t tmp; + while (((tmp=(n>>1))<<1)==n) + { result+=2; n=tmp; } + + size_t limit=(size_t)sqrt(n+0.01); + for (size_t x=3; x<=limit; x+=2) + while ((tmp=(n/x))*x==n) + { + result+= (x<=5) ? x : lfp*x; // penalize larger prime factors + n=tmp; + limit=(size_t)sqrt(n+0.01); + } + if (n>1) result+=(n<=5) ? n : lfp*n; + + return result*ni; + } + +/* returns the smallest composite of 2, 3, 5, 7 and 11 which is >= n */ +NOINLINE static size_t good_size(size_t n) + { + if (n<=6) return n; + + size_t bestfac=2*n; + for (size_t f2=1; f2<bestfac; f2*=2) + for (size_t f23=f2; f23<bestfac; f23*=3) + for (size_t f235=f23; f235<bestfac; f235*=5) + for (size_t f2357=f235; f2357<bestfac; f2357*=7) + for (size_t f235711=f2357; f235711<bestfac; f235711*=11) + if (f235711>=n) bestfac=f235711; + return bestfac; + } + +typedef struct cmplx { + double r,i; +} cmplx; + +#define NFCT 25 +typedef struct cfftp_fctdata + { + size_t fct; + cmplx *tw, *tws; + } cfftp_fctdata; + +typedef struct cfftp_plan_i + { + size_t length, nfct; + cmplx *mem; + cfftp_fctdata fct[NFCT]; + } cfftp_plan_i; +typedef struct cfftp_plan_i * cfftp_plan; + +#define PMC(a,b,c,d) { a.r=c.r+d.r; a.i=c.i+d.i; b.r=c.r-d.r; b.i=c.i-d.i; } +#define ADDC(a,b,c) { a.r=b.r+c.r; a.i=b.i+c.i; } +#define SCALEC(a,b) { a.r*=b; a.i*=b; } +#define ROT90(a) { double tmp_=a.r; a.r=-a.i; a.i=tmp_; } +#define ROTM90(a) { double tmp_=-a.r; a.r=a.i; a.i=tmp_; } +#define CH(a,b,c) ch[(a)+ido*((b)+l1*(c))] +#define CC(a,b,c) cc[(a)+ido*((b)+cdim*(c))] +#define WA(x,i) wa[(i)-1+(x)*(ido-1)] +/* a = b*c */ +#define A_EQ_B_MUL_C(a,b,c) { a.r=b.r*c.r-b.i*c.i; a.i=b.r*c.i+b.i*c.r; } +/* a = conj(b)*c*/ +#define A_EQ_CB_MUL_C(a,b,c) { a.r=b.r*c.r+b.i*c.i; a.i=b.r*c.i-b.i*c.r; } + +#define PMSIGNC(a,b,c,d) { a.r=c.r+sign*d.r; a.i=c.i+sign*d.i; b.r=c.r-sign*d.r; b.i=c.i-sign*d.i; } +/* a = b*c */ +#define MULPMSIGNC(a,b,c) { a.r=b.r*c.r-sign*b.i*c.i; a.i=b.r*c.i+sign*b.i*c.r; } +/* a *= b */ +#define MULPMSIGNCEQ(a,b) { double xtmp=a.r; a.r=b.r*a.r-sign*b.i*a.i; a.i=b.r*a.i+sign*b.i*xtmp; } + +NOINLINE static void pass2b (size_t ido, size_t l1, const cmplx * restrict cc, + cmplx * restrict ch, const cmplx * restrict wa) + { + const size_t cdim=2; + + if (ido==1) + for (size_t k=0; k<l1; ++k) + PMC (CH(0,k,0),CH(0,k,1),CC(0,0,k),CC(0,1,k)) + else + for (size_t k=0; k<l1; ++k) + { + PMC (CH(0,k,0),CH(0,k,1),CC(0,0,k),CC(0,1,k)) + for (size_t i=1; i<ido; ++i) + { + cmplx t; + PMC (CH(i,k,0),t,CC(i,0,k),CC(i,1,k)) + A_EQ_B_MUL_C (CH(i,k,1),WA(0,i),t) + } + } + } + +NOINLINE static void pass2f (size_t ido, size_t l1, const cmplx * restrict cc, + cmplx * restrict ch, const cmplx * restrict wa) + { + const size_t cdim=2; + + if (ido==1) + for (size_t k=0; k<l1; ++k) + PMC (CH(0,k,0),CH(0,k,1),CC(0,0,k),CC(0,1,k)) + else + for (size_t k=0; k<l1; ++k) + { + PMC (CH(0,k,0),CH(0,k,1),CC(0,0,k),CC(0,1,k)) + for (size_t i=1; i<ido; ++i) + { + cmplx t; + PMC (CH(i,k,0),t,CC(i,0,k),CC(i,1,k)) + A_EQ_CB_MUL_C (CH(i,k,1),WA(0,i),t) + } + } + } + +#define PREP3(idx) \ + cmplx t0 = CC(idx,0,k), t1, t2; \ + PMC (t1,t2,CC(idx,1,k),CC(idx,2,k)) \ + CH(idx,k,0).r=t0.r+t1.r; \ + CH(idx,k,0).i=t0.i+t1.i; +#define PARTSTEP3a(u1,u2,twr,twi) \ + { \ + cmplx ca,cb; \ + ca.r=t0.r+twr*t1.r; \ + ca.i=t0.i+twr*t1.i; \ + cb.i=twi*t2.r; \ + cb.r=-(twi*t2.i); \ + PMC(CH(0,k,u1),CH(0,k,u2),ca,cb) \ + } + +#define PARTSTEP3b(u1,u2,twr,twi) \ + { \ + cmplx ca,cb,da,db; \ + ca.r=t0.r+twr*t1.r; \ + ca.i=t0.i+twr*t1.i; \ + cb.i=twi*t2.r; \ + cb.r=-(twi*t2.i); \ + PMC(da,db,ca,cb) \ + A_EQ_B_MUL_C (CH(i,k,u1),WA(u1-1,i),da) \ + A_EQ_B_MUL_C (CH(i,k,u2),WA(u2-1,i),db) \ + } +NOINLINE static void pass3b (size_t ido, size_t l1, const cmplx * restrict cc, + cmplx * restrict ch, const cmplx * restrict wa) + { + const size_t cdim=3; + const double tw1r=-0.5, tw1i= 0.86602540378443864676; + + if (ido==1) + for (size_t k=0; k<l1; ++k) + { + PREP3(0) + PARTSTEP3a(1,2,tw1r,tw1i) + } + else + for (size_t k=0; k<l1; ++k) + { + { + PREP3(0) + PARTSTEP3a(1,2,tw1r,tw1i) + } + for (size_t i=1; i<ido; ++i) + { + PREP3(i) + PARTSTEP3b(1,2,tw1r,tw1i) + } + } + } +#define PARTSTEP3f(u1,u2,twr,twi) \ + { \ + cmplx ca,cb,da,db; \ + ca.r=t0.r+twr*t1.r; \ + ca.i=t0.i+twr*t1.i; \ + cb.i=twi*t2.r; \ + cb.r=-(twi*t2.i); \ + PMC(da,db,ca,cb) \ + A_EQ_CB_MUL_C (CH(i,k,u1),WA(u1-1,i),da) \ + A_EQ_CB_MUL_C (CH(i,k,u2),WA(u2-1,i),db) \ + } +NOINLINE static void pass3f (size_t ido, size_t l1, const cmplx * restrict cc, + cmplx * restrict ch, const cmplx * restrict wa) + { + const size_t cdim=3; + const double tw1r=-0.5, tw1i= -0.86602540378443864676; + + if (ido==1) + for (size_t k=0; k<l1; ++k) + { + PREP3(0) + PARTSTEP3a(1,2,tw1r,tw1i) + } + else + for (size_t k=0; k<l1; ++k) + { + { + PREP3(0) + PARTSTEP3a(1,2,tw1r,tw1i) + } + for (size_t i=1; i<ido; ++i) + { + PREP3(i) + PARTSTEP3f(1,2,tw1r,tw1i) + } + } + } + +NOINLINE static void pass4b (size_t ido, size_t l1, const cmplx * restrict cc, + cmplx * restrict ch, const cmplx * restrict wa) + { + const size_t cdim=4; + + if (ido==1) + for (size_t k=0; k<l1; ++k) + { + cmplx t1, t2, t3, t4; + PMC(t2,t1,CC(0,0,k),CC(0,2,k)) + PMC(t3,t4,CC(0,1,k),CC(0,3,k)) + ROT90(t4) + PMC(CH(0,k,0),CH(0,k,2),t2,t3) + PMC(CH(0,k,1),CH(0,k,3),t1,t4) + } + else + for (size_t k=0; k<l1; ++k) + { + { + cmplx t1, t2, t3, t4; + PMC(t2,t1,CC(0,0,k),CC(0,2,k)) + PMC(t3,t4,CC(0,1,k),CC(0,3,k)) + ROT90(t4) + PMC(CH(0,k,0),CH(0,k,2),t2,t3) + PMC(CH(0,k,1),CH(0,k,3),t1,t4) + } + for (size_t i=1; i<ido; ++i) + { + cmplx c2, c3, c4, t1, t2, t3, t4; + cmplx cc0=CC(i,0,k), cc1=CC(i,1,k),cc2=CC(i,2,k),cc3=CC(i,3,k); + PMC(t2,t1,cc0,cc2) + PMC(t3,t4,cc1,cc3) + ROT90(t4) + cmplx wa0=WA(0,i), wa1=WA(1,i),wa2=WA(2,i); + PMC(CH(i,k,0),c3,t2,t3) + PMC(c2,c4,t1,t4) + A_EQ_B_MUL_C (CH(i,k,1),wa0,c2) + A_EQ_B_MUL_C (CH(i,k,2),wa1,c3) + A_EQ_B_MUL_C (CH(i,k,3),wa2,c4) + } + } + } +NOINLINE static void pass4f (size_t ido, size_t l1, const cmplx * restrict cc, + cmplx * restrict ch, const cmplx * restrict wa) + { + const size_t cdim=4; + + if (ido==1) + for (size_t k=0; k<l1; ++k) + { + cmplx t1, t2, t3, t4; + PMC(t2,t1,CC(0,0,k),CC(0,2,k)) + PMC(t3,t4,CC(0,1,k),CC(0,3,k)) + ROTM90(t4) + PMC(CH(0,k,0),CH(0,k,2),t2,t3) + PMC(CH(0,k,1),CH(0,k,3),t1,t4) + } + else + for (size_t k=0; k<l1; ++k) + { + { + cmplx t1, t2, t3, t4; + PMC(t2,t1,CC(0,0,k),CC(0,2,k)) + PMC(t3,t4,CC(0,1,k),CC(0,3,k)) + ROTM90(t4) + PMC(CH(0,k,0),CH(0,k,2),t2,t3) + PMC (CH(0,k,1),CH(0,k,3),t1,t4) + } + for (size_t i=1; i<ido; ++i) + { + cmplx c2, c3, c4, t1, t2, t3, t4; + cmplx cc0=CC(i,0,k), cc1=CC(i,1,k),cc2=CC(i,2,k),cc3=CC(i,3,k); + PMC(t2,t1,cc0,cc2) + PMC(t3,t4,cc1,cc3) + ROTM90(t4) + cmplx wa0=WA(0,i), wa1=WA(1,i),wa2=WA(2,i); + PMC(CH(i,k,0),c3,t2,t3) + PMC(c2,c4,t1,t4) + A_EQ_CB_MUL_C (CH(i,k,1),wa0,c2) + A_EQ_CB_MUL_C (CH(i,k,2),wa1,c3) + A_EQ_CB_MUL_C (CH(i,k,3),wa2,c4) + } + } + } + +#define PREP5(idx) \ + cmplx t0 = CC(idx,0,k), t1, t2, t3, t4; \ + PMC (t1,t4,CC(idx,1,k),CC(idx,4,k)) \ + PMC (t2,t3,CC(idx,2,k),CC(idx,3,k)) \ + CH(idx,k,0).r=t0.r+t1.r+t2.r; \ + CH(idx,k,0).i=t0.i+t1.i+t2.i; + +#define PARTSTEP5a(u1,u2,twar,twbr,twai,twbi) \ + { \ + cmplx ca,cb; \ + ca.r=t0.r+twar*t1.r+twbr*t2.r; \ + ca.i=t0.i+twar*t1.i+twbr*t2.i; \ + cb.i=twai*t4.r twbi*t3.r; \ + cb.r=-(twai*t4.i twbi*t3.i); \ + PMC(CH(0,k,u1),CH(0,k,u2),ca,cb) \ + } + +#define PARTSTEP5b(u1,u2,twar,twbr,twai,twbi) \ + { \ + cmplx ca,cb,da,db; \ + ca.r=t0.r+twar*t1.r+twbr*t2.r; \ + ca.i=t0.i+twar*t1.i+twbr*t2.i; \ + cb.i=twai*t4.r twbi*t3.r; \ + cb.r=-(twai*t4.i twbi*t3.i); \ + PMC(da,db,ca,cb) \ + A_EQ_B_MUL_C (CH(i,k,u1),WA(u1-1,i),da) \ + A_EQ_B_MUL_C (CH(i,k,u2),WA(u2-1,i),db) \ + } +NOINLINE static void pass5b (size_t ido, size_t l1, const cmplx * restrict cc, + cmplx * restrict ch, const cmplx * restrict wa) + { + const size_t cdim=5; + const double tw1r= 0.3090169943749474241, + tw1i= 0.95105651629515357212, + tw2r= -0.8090169943749474241, + tw2i= 0.58778525229247312917; + + if (ido==1) + for (size_t k=0; k<l1; ++k) + { + PREP5(0) + PARTSTEP5a(1,4,tw1r,tw2r,+tw1i,+tw2i) + PARTSTEP5a(2,3,tw2r,tw1r,+tw2i,-tw1i) + } + else + for (size_t k=0; k<l1; ++k) + { + { + PREP5(0) + PARTSTEP5a(1,4,tw1r,tw2r,+tw1i,+tw2i) + PARTSTEP5a(2,3,tw2r,tw1r,+tw2i,-tw1i) + } + for (size_t i=1; i<ido; ++i) + { + PREP5(i) + PARTSTEP5b(1,4,tw1r,tw2r,+tw1i,+tw2i) + PARTSTEP5b(2,3,tw2r,tw1r,+tw2i,-tw1i) + } + } + } +#define PARTSTEP5f(u1,u2,twar,twbr,twai,twbi) \ + { \ + cmplx ca,cb,da,db; \ + ca.r=t0.r+twar*t1.r+twbr*t2.r; \ + ca.i=t0.i+twar*t1.i+twbr*t2.i; \ + cb.i=twai*t4.r twbi*t3.r; \ + cb.r=-(twai*t4.i twbi*t3.i); \ + PMC(da,db,ca,cb) \ + A_EQ_CB_MUL_C (CH(i,k,u1),WA(u1-1,i),da) \ + A_EQ_CB_MUL_C (CH(i,k,u2),WA(u2-1,i),db) \ + } +NOINLINE static void pass5f (size_t ido, size_t l1, const cmplx * restrict cc, + cmplx * restrict ch, const cmplx * restrict wa) + { + const size_t cdim=5; + const double tw1r= 0.3090169943749474241, + tw1i= -0.95105651629515357212, + tw2r= -0.8090169943749474241, + tw2i= -0.58778525229247312917; + + if (ido==1) + for (size_t k=0; k<l1; ++k) + { + PREP5(0) + PARTSTEP5a(1,4,tw1r,tw2r,+tw1i,+tw2i) + PARTSTEP5a(2,3,tw2r,tw1r,+tw2i,-tw1i) + } + else + for (size_t k=0; k<l1; ++k) + { + { + PREP5(0) + PARTSTEP5a(1,4,tw1r,tw2r,+tw1i,+tw2i) + PARTSTEP5a(2,3,tw2r,tw1r,+tw2i,-tw1i) + } + for (size_t i=1; i<ido; ++i) + { + PREP5(i) + PARTSTEP5f(1,4,tw1r,tw2r,+tw1i,+tw2i) + PARTSTEP5f(2,3,tw2r,tw1r,+tw2i,-tw1i) + } + } + } + +#define PREP7(idx) \ + cmplx t1 = CC(idx,0,k), t2, t3, t4, t5, t6, t7; \ + PMC (t2,t7,CC(idx,1,k),CC(idx,6,k)) \ + PMC (t3,t6,CC(idx,2,k),CC(idx,5,k)) \ + PMC (t4,t5,CC(idx,3,k),CC(idx,4,k)) \ + CH(idx,k,0).r=t1.r+t2.r+t3.r+t4.r; \ + CH(idx,k,0).i=t1.i+t2.i+t3.i+t4.i; + +#define PARTSTEP7a0(u1,u2,x1,x2,x3,y1,y2,y3,out1,out2) \ + { \ + cmplx ca,cb; \ + ca.r=t1.r+x1*t2.r+x2*t3.r+x3*t4.r; \ + ca.i=t1.i+x1*t2.i+x2*t3.i+x3*t4.i; \ + cb.i=y1*t7.r y2*t6.r y3*t5.r; \ + cb.r=-(y1*t7.i y2*t6.i y3*t5.i); \ + PMC(out1,out2,ca,cb) \ + } +#define PARTSTEP7a(u1,u2,x1,x2,x3,y1,y2,y3) \ + PARTSTEP7a0(u1,u2,x1,x2,x3,y1,y2,y3,CH(0,k,u1),CH(0,k,u2)) +#define PARTSTEP7(u1,u2,x1,x2,x3,y1,y2,y3) \ + { \ + cmplx da,db; \ + PARTSTEP7a0(u1,u2,x1,x2,x3,y1,y2,y3,da,db) \ + MULPMSIGNC (CH(i,k,u1),WA(u1-1,i),da) \ + MULPMSIGNC (CH(i,k,u2),WA(u2-1,i),db) \ + } + +NOINLINE static void pass7(size_t ido, size_t l1, const cmplx * restrict cc, + cmplx * restrict ch, const cmplx * restrict wa, const int sign) + { + const size_t cdim=7; + const double tw1r= 0.623489801858733530525, + tw1i= sign * 0.7818314824680298087084, + tw2r= -0.222520933956314404289, + tw2i= sign * 0.9749279121818236070181, + tw3r= -0.9009688679024191262361, + tw3i= sign * 0.4338837391175581204758; + + if (ido==1) + for (size_t k=0; k<l1; ++k) + { + PREP7(0) + PARTSTEP7a(1,6,tw1r,tw2r,tw3r,+tw1i,+tw2i,+tw3i) + PARTSTEP7a(2,5,tw2r,tw3r,tw1r,+tw2i,-tw3i,-tw1i) + PARTSTEP7a(3,4,tw3r,tw1r,tw2r,+tw3i,-tw1i,+tw2i) + } + else + for (size_t k=0; k<l1; ++k) + { + { + PREP7(0) + PARTSTEP7a(1,6,tw1r,tw2r,tw3r,+tw1i,+tw2i,+tw3i) + PARTSTEP7a(2,5,tw2r,tw3r,tw1r,+tw2i,-tw3i,-tw1i) + PARTSTEP7a(3,4,tw3r,tw1r,tw2r,+tw3i,-tw1i,+tw2i) + } + for (size_t i=1; i<ido; ++i) + { + PREP7(i) + PARTSTEP7(1,6,tw1r,tw2r,tw3r,+tw1i,+tw2i,+tw3i) + PARTSTEP7(2,5,tw2r,tw3r,tw1r,+tw2i,-tw3i,-tw1i) + PARTSTEP7(3,4,tw3r,tw1r,tw2r,+tw3i,-tw1i,+tw2i) + } + } + } + +#define PREP11(idx) \ + cmplx t1 = CC(idx,0,k), t2, t3, t4, t5, t6, t7, t8, t9, t10, t11; \ + PMC (t2,t11,CC(idx,1,k),CC(idx,10,k)) \ + PMC (t3,t10,CC(idx,2,k),CC(idx, 9,k)) \ + PMC (t4,t9 ,CC(idx,3,k),CC(idx, 8,k)) \ + PMC (t5,t8 ,CC(idx,4,k),CC(idx, 7,k)) \ + PMC (t6,t7 ,CC(idx,5,k),CC(idx, 6,k)) \ + CH(idx,k,0).r=t1.r+t2.r+t3.r+t4.r+t5.r+t6.r; \ + CH(idx,k,0).i=t1.i+t2.i+t3.i+t4.i+t5.i+t6.i; + +#define PARTSTEP11a0(u1,u2,x1,x2,x3,x4,x5,y1,y2,y3,y4,y5,out1,out2) \ + { \ + cmplx ca,cb; \ + ca.r=t1.r+x1*t2.r+x2*t3.r+x3*t4.r+x4*t5.r+x5*t6.r; \ + ca.i=t1.i+x1*t2.i+x2*t3.i+x3*t4.i+x4*t5.i+x5*t6.i; \ + cb.i=y1*t11.r y2*t10.r y3*t9.r y4*t8.r y5*t7.r; \ + cb.r=-(y1*t11.i y2*t10.i y3*t9.i y4*t8.i y5*t7.i ); \ + PMC(out1,out2,ca,cb) \ + } +#define PARTSTEP11a(u1,u2,x1,x2,x3,x4,x5,y1,y2,y3,y4,y5) \ + PARTSTEP11a0(u1,u2,x1,x2,x3,x4,x5,y1,y2,y3,y4,y5,CH(0,k,u1),CH(0,k,u2)) +#define PARTSTEP11(u1,u2,x1,x2,x3,x4,x5,y1,y2,y3,y4,y5) \ + { \ + cmplx da,db; \ + PARTSTEP11a0(u1,u2,x1,x2,x3,x4,x5,y1,y2,y3,y4,y5,da,db) \ + MULPMSIGNC (CH(i,k,u1),WA(u1-1,i),da) \ + MULPMSIGNC (CH(i,k,u2),WA(u2-1,i),db) \ + } + +NOINLINE static void pass11 (size_t ido, size_t l1, const cmplx * restrict cc, + cmplx * restrict ch, const cmplx * restrict wa, const int sign) + { + const size_t cdim=11; + const double tw1r = 0.8412535328311811688618, + tw1i = sign * 0.5406408174555975821076, + tw2r = 0.4154150130018864255293, + tw2i = sign * 0.9096319953545183714117, + tw3r = -0.1423148382732851404438, + tw3i = sign * 0.9898214418809327323761, + tw4r = -0.6548607339452850640569, + tw4i = sign * 0.755749574354258283774, + tw5r = -0.9594929736144973898904, + tw5i = sign * 0.2817325568414296977114; + + if (ido==1) + for (size_t k=0; k<l1; ++k) + { + PREP11(0) + PARTSTEP11a(1,10,tw1r,tw2r,tw3r,tw4r,tw5r,+tw1i,+tw2i,+tw3i,+tw4i,+tw5i) + PARTSTEP11a(2, 9,tw2r,tw4r,tw5r,tw3r,tw1r,+tw2i,+tw4i,-tw5i,-tw3i,-tw1i) + PARTSTEP11a(3, 8,tw3r,tw5r,tw2r,tw1r,tw4r,+tw3i,-tw5i,-tw2i,+tw1i,+tw4i) + PARTSTEP11a(4, 7,tw4r,tw3r,tw1r,tw5r,tw2r,+tw4i,-tw3i,+tw1i,+tw5i,-tw2i) + PARTSTEP11a(5, 6,tw5r,tw1r,tw4r,tw2r,tw3r,+tw5i,-tw1i,+tw4i,-tw2i,+tw3i) + } + else + for (size_t k=0; k<l1; ++k) + { + { + PREP11(0) + PARTSTEP11a(1,10,tw1r,tw2r,tw3r,tw4r,tw5r,+tw1i,+tw2i,+tw3i,+tw4i,+tw5i) + PARTSTEP11a(2, 9,tw2r,tw4r,tw5r,tw3r,tw1r,+tw2i,+tw4i,-tw5i,-tw3i,-tw1i) + PARTSTEP11a(3, 8,tw3r,tw5r,tw2r,tw1r,tw4r,+tw3i,-tw5i,-tw2i,+tw1i,+tw4i) + PARTSTEP11a(4, 7,tw4r,tw3r,tw1r,tw5r,tw2r,+tw4i,-tw3i,+tw1i,+tw5i,-tw2i) + PARTSTEP11a(5, 6,tw5r,tw1r,tw4r,tw2r,tw3r,+tw5i,-tw1i,+tw4i,-tw2i,+tw3i) + } + for (size_t i=1; i<ido; ++i) + { + PREP11(i) + PARTSTEP11(1,10,tw1r,tw2r,tw3r,tw4r,tw5r,+tw1i,+tw2i,+tw3i,+tw4i,+tw5i) + PARTSTEP11(2, 9,tw2r,tw4r,tw5r,tw3r,tw1r,+tw2i,+tw4i,-tw5i,-tw3i,-tw1i) + PARTSTEP11(3, 8,tw3r,tw5r,tw2r,tw1r,tw4r,+tw3i,-tw5i,-tw2i,+tw1i,+tw4i) + PARTSTEP11(4, 7,tw4r,tw3r,tw1r,tw5r,tw2r,+tw4i,-tw3i,+tw1i,+tw5i,-tw2i) + PARTSTEP11(5, 6,tw5r,tw1r,tw4r,tw2r,tw3r,+tw5i,-tw1i,+tw4i,-tw2i,+tw3i) + } + } + } + +#define CX(a,b,c) cc[(a)+ido*((b)+l1*(c))] +#define CX2(a,b) cc[(a)+idl1*(b)] +#define CH2(a,b) ch[(a)+idl1*(b)] + +NOINLINE static int passg (size_t ido, size_t ip, size_t l1, + cmplx * restrict cc, cmplx * restrict ch, const cmplx * restrict wa, + const cmplx * restrict csarr, const int sign) + { + const size_t cdim=ip; + size_t ipph = (ip+1)/2; + size_t idl1 = ido*l1; + + cmplx * restrict wal=RALLOC(cmplx,ip); + if (!wal) return -1; + wal[0]=(cmplx){1.,0.}; + for (size_t i=1; i<ip; ++i) + wal[i]=(cmplx){csarr[i].r,sign*csarr[i].i}; + + for (size_t k=0; k<l1; ++k) + for (size_t i=0; i<ido; ++i) + CH(i,k,0) = CC(i,0,k); + for (size_t j=1, jc=ip-1; j<ipph; ++j, --jc) + for (size_t k=0; k<l1; ++k) + for (size_t i=0; i<ido; ++i) + PMC(CH(i,k,j),CH(i,k,jc),CC(i,j,k),CC(i,jc,k)) + for (size_t k=0; k<l1; ++k) + for (size_t i=0; i<ido; ++i) + { + cmplx tmp = CH(i,k,0); + for (size_t j=1; j<ipph; ++j) + ADDC(tmp,tmp,CH(i,k,j)) + CX(i,k,0) = tmp; + } + for (size_t l=1, lc=ip-1; l<ipph; ++l, --lc) + { + // j=0 + for (size_t ik=0; ik<idl1; ++ik) + { + CX2(ik,l).r = CH2(ik,0).r+wal[l].r*CH2(ik,1).r+wal[2*l].r*CH2(ik,2).r; + CX2(ik,l).i = CH2(ik,0).i+wal[l].r*CH2(ik,1).i+wal[2*l].r*CH2(ik,2).i; + CX2(ik,lc).r=-wal[l].i*CH2(ik,ip-1).i-wal[2*l].i*CH2(ik,ip-2).i; + CX2(ik,lc).i=wal[l].i*CH2(ik,ip-1).r+wal[2*l].i*CH2(ik,ip-2).r; + } + + size_t iwal=2*l; + size_t j=3, jc=ip-3; + for (; j<ipph-1; j+=2, jc-=2) + { + iwal+=l; if (iwal>ip) iwal-=ip; + cmplx xwal=wal[iwal]; + iwal+=l; if (iwal>ip) iwal-=ip; + cmplx xwal2=wal[iwal]; + for (size_t ik=0; ik<idl1; ++ik) + { + CX2(ik,l).r += CH2(ik,j).r*xwal.r+CH2(ik,j+1).r*xwal2.r; + CX2(ik,l).i += CH2(ik,j).i*xwal.r+CH2(ik,j+1).i*xwal2.r; + CX2(ik,lc).r -= CH2(ik,jc).i*xwal.i+CH2(ik,jc-1).i*xwal2.i; + CX2(ik,lc).i += CH2(ik,jc).r*xwal.i+CH2(ik,jc-1).r*xwal2.i; + } + } + for (; j<ipph; ++j, --jc) + { + iwal+=l; if (iwal>ip) iwal-=ip; + cmplx xwal=wal[iwal]; + for (size_t ik=0; ik<idl1; ++ik) + { + CX2(ik,l).r += CH2(ik,j).r*xwal.r; + CX2(ik,l).i += CH2(ik,j).i*xwal.r; + CX2(ik,lc).r -= CH2(ik,jc).i*xwal.i; + CX2(ik,lc).i += CH2(ik,jc).r*xwal.i; + } + } + } + DEALLOC(wal); + + // shuffling and twiddling + if (ido==1) + for (size_t j=1, jc=ip-1; j<ipph; ++j, --jc) + for (size_t ik=0; ik<idl1; ++ik) + { + cmplx t1=CX2(ik,j), t2=CX2(ik,jc); + PMC(CX2(ik,j),CX2(ik,jc),t1,t2) + } + else + { + for (size_t j=1, jc=ip-1; j<ipph; ++j,--jc) + for (size_t k=0; k<l1; ++k) + { + cmplx t1=CX(0,k,j), t2=CX(0,k,jc); + PMC(CX(0,k,j),CX(0,k,jc),t1,t2) + for (size_t i=1; i<ido; ++i) + { + cmplx x1, x2; + PMC(x1,x2,CX(i,k,j),CX(i,k,jc)) + size_t idij=(j-1)*(ido-1)+i-1; + MULPMSIGNC (CX(i,k,j),wa[idij],x1) + idij=(jc-1)*(ido-1)+i-1; + MULPMSIGNC (CX(i,k,jc),wa[idij],x2) + } + } + } + return 0; + } + +#undef CH2 +#undef CX2 +#undef CX + +NOINLINE WARN_UNUSED_RESULT static int pass_all(cfftp_plan plan, cmplx c[], double fct, + const int sign) + { + if (plan->length==1) return 0; + size_t len=plan->length; + size_t l1=1, nf=plan->nfct; + cmplx *ch = RALLOC(cmplx, len); + if (!ch) return -1; + cmplx *p1=c, *p2=ch; + + for(size_t k1=0; k1<nf; k1++) + { + size_t ip=plan->fct[k1].fct; + size_t l2=ip*l1; + size_t ido = len/l2; + if (ip==4) + sign>0 ? pass4b (ido, l1, p1, p2, plan->fct[k1].tw) + : pass4f (ido, l1, p1, p2, plan->fct[k1].tw); + else if(ip==2) + sign>0 ? pass2b (ido, l1, p1, p2, plan->fct[k1].tw) + : pass2f (ido, l1, p1, p2, plan->fct[k1].tw); + else if(ip==3) + sign>0 ? pass3b (ido, l1, p1, p2, plan->fct[k1].tw) + : pass3f (ido, l1, p1, p2, plan->fct[k1].tw); + else if(ip==5) + sign>0 ? pass5b (ido, l1, p1, p2, plan->fct[k1].tw) + : pass5f (ido, l1, p1, p2, plan->fct[k1].tw); + else if(ip==7) pass7 (ido, l1, p1, p2, plan->fct[k1].tw, sign); + else if(ip==11) pass11(ido, l1, p1, p2, plan->fct[k1].tw, sign); + else + { + if (passg(ido, ip, l1, p1, p2, plan->fct[k1].tw, plan->fct[k1].tws, sign)) + { DEALLOC(ch); return -1; } + SWAP(p1,p2,cmplx *); + } + SWAP(p1,p2,cmplx *); + l1=l2; + } + if (p1!=c) + { + if (fct!=1.) + for (size_t i=0; i<len; ++i) + { + c[i].r = ch[i].r*fct; + c[i].i = ch[i].i*fct; + } + else + memcpy (c,p1,len*sizeof(cmplx)); + } + else + if (fct!=1.) + for (size_t i=0; i<len; ++i) + { + c[i].r *= fct; + c[i].i *= fct; + } + DEALLOC(ch); + return 0; + } + +#undef PMSIGNC +#undef A_EQ_B_MUL_C +#undef A_EQ_CB_MUL_C +#undef MULPMSIGNC +#undef MULPMSIGNCEQ + +#undef WA +#undef CC +#undef CH +#undef ROT90 +#undef SCALEC +#undef ADDC +#undef PMC + +NOINLINE WARN_UNUSED_RESULT +static int cfftp_forward(cfftp_plan plan, double c[], double fct) + { return pass_all(plan,(cmplx *)c, fct, -1); } + +NOINLINE WARN_UNUSED_RESULT +static int cfftp_backward(cfftp_plan plan, double c[], double fct) + { return pass_all(plan,(cmplx *)c, fct, 1); } + +NOINLINE WARN_UNUSED_RESULT +static int cfftp_factorize (cfftp_plan plan) + { + size_t length=plan->length; + size_t nfct=0; + while ((length%4)==0) + { if (nfct>=NFCT) return -1; plan->fct[nfct++].fct=4; length>>=2; } + if ((length%2)==0) + { + length>>=1; + // factor 2 should be at the front of the factor list + if (nfct>=NFCT) return -1; + plan->fct[nfct++].fct=2; + SWAP(plan->fct[0].fct, plan->fct[nfct-1].fct,size_t); + } + size_t maxl=(size_t)(sqrt((double)length))+1; + for (size_t divisor=3; (length>1)&&(divisor<maxl); divisor+=2) + if ((length%divisor)==0) + { + while ((length%divisor)==0) + { + if (nfct>=NFCT) return -1; + plan->fct[nfct++].fct=divisor; + length/=divisor; + } + maxl=(size_t)(sqrt((double)length))+1; + } + if (length>1) plan->fct[nfct++].fct=length; + plan->nfct=nfct; + return 0; + } + +NOINLINE static size_t cfftp_twsize (cfftp_plan plan) + { + size_t twsize=0, l1=1; + for (size_t k=0; k<plan->nfct; ++k) + { + size_t ip=plan->fct[k].fct, ido= plan->length/(l1*ip); + twsize+=(ip-1)*(ido-1); + if (ip>11) + twsize+=ip; + l1*=ip; + } + return twsize; + } + +NOINLINE WARN_UNUSED_RESULT static int cfftp_comp_twiddle (cfftp_plan plan) + { + size_t length=plan->length; + double *twid = RALLOC(double, 2*length); + if (!twid) return -1; + sincos_2pibyn(length, twid); + size_t l1=1; + size_t memofs=0; + for (size_t k=0; k<plan->nfct; ++k) + { + size_t ip=plan->fct[k].fct, ido= length/(l1*ip); + plan->fct[k].tw=plan->mem+memofs; + memofs+=(ip-1)*(ido-1); + for (size_t j=1; j<ip; ++j) + for (size_t i=1; i<ido; ++i) + { + plan->fct[k].tw[(j-1)*(ido-1)+i-1].r = twid[2*j*l1*i]; + plan->fct[k].tw[(j-1)*(ido-1)+i-1].i = twid[2*j*l1*i+1]; + } + if (ip>11) + { + plan->fct[k].tws=plan->mem+memofs; + memofs+=ip; + for (size_t j=0; j<ip; ++j) + { + plan->fct[k].tws[j].r = twid[2*j*l1*ido]; + plan->fct[k].tws[j].i = twid[2*j*l1*ido+1]; + } + } + l1*=ip; + } + DEALLOC(twid); + return 0; + } + +static cfftp_plan make_cfftp_plan (size_t length) + { + if (length==0) return NULL; + cfftp_plan plan = RALLOC(cfftp_plan_i,1); + if (!plan) return NULL; + plan->length=length; + plan->nfct=0; + for (size_t i=0; i<NFCT; ++i) + plan->fct[i]=(cfftp_fctdata){0,0,0}; + plan->mem=0; + if (length==1) return plan; + if (cfftp_factorize(plan)!=0) { DEALLOC(plan); return NULL; } + size_t tws=cfftp_twsize(plan); + plan->mem=RALLOC(cmplx,tws); + if (!plan->mem) { DEALLOC(plan); return NULL; } + if (cfftp_comp_twiddle(plan)!=0) + { DEALLOC(plan->mem); DEALLOC(plan); return NULL; } + return plan; + } + +static void destroy_cfftp_plan (cfftp_plan plan) + { + DEALLOC(plan->mem); + DEALLOC(plan); + } + +typedef struct rfftp_fctdata + { + size_t fct; + double *tw, *tws; + } rfftp_fctdata; + +typedef struct rfftp_plan_i + { + size_t length, nfct; + double *mem; + rfftp_fctdata fct[NFCT]; + } rfftp_plan_i; +typedef struct rfftp_plan_i * rfftp_plan; + +#define WA(x,i) wa[(i)+(x)*(ido-1)] +#define PM(a,b,c,d) { a=c+d; b=c-d; } +/* (a+ib) = conj(c+id) * (e+if) */ +#define MULPM(a,b,c,d,e,f) { a=c*e+d*f; b=c*f-d*e; } + +#define CC(a,b,c) cc[(a)+ido*((b)+l1*(c))] +#define CH(a,b,c) ch[(a)+ido*((b)+cdim*(c))] + +NOINLINE static void radf2 (size_t ido, size_t l1, const double * restrict cc, + double * restrict ch, const double * restrict wa) + { + const size_t cdim=2; + + for (size_t k=0; k<l1; k++) + PM (CH(0,0,k),CH(ido-1,1,k),CC(0,k,0),CC(0,k,1)) + if ((ido&1)==0) + for (size_t k=0; k<l1; k++) + { + CH( 0,1,k) = -CC(ido-1,k,1); + CH(ido-1,0,k) = CC(ido-1,k,0); + } + if (ido<=2) return; + for (size_t k=0; k<l1; k++) + for (size_t i=2; i<ido; i+=2) + { + size_t ic=ido-i; + double tr2, ti2; + MULPM (tr2,ti2,WA(0,i-2),WA(0,i-1),CC(i-1,k,1),CC(i,k,1)) + PM (CH(i-1,0,k),CH(ic-1,1,k),CC(i-1,k,0),tr2) + PM (CH(i ,0,k),CH(ic ,1,k),ti2,CC(i ,k,0)) + } + } + +NOINLINE static void radf3(size_t ido, size_t l1, const double * restrict cc, + double * restrict ch, const double * restrict wa) + { + const size_t cdim=3; + static const double taur=-0.5, taui=0.86602540378443864676; + + for (size_t k=0; k<l1; k++) + { + double cr2=CC(0,k,1)+CC(0,k,2); + CH(0,0,k) = CC(0,k,0)+cr2; + CH(0,2,k) = taui*(CC(0,k,2)-CC(0,k,1)); + CH(ido-1,1,k) = CC(0,k,0)+taur*cr2; + } + if (ido==1) return; + for (size_t k=0; k<l1; k++) + for (size_t i=2; i<ido; i+=2) + { + size_t ic=ido-i; + double di2, di3, dr2, dr3; + MULPM (dr2,di2,WA(0,i-2),WA(0,i-1),CC(i-1,k,1),CC(i,k,1)) // d2=conj(WA0)*CC1 + MULPM (dr3,di3,WA(1,i-2),WA(1,i-1),CC(i-1,k,2),CC(i,k,2)) // d3=conj(WA1)*CC2 + double cr2=dr2+dr3; // c add + double ci2=di2+di3; + CH(i-1,0,k) = CC(i-1,k,0)+cr2; // c add + CH(i ,0,k) = CC(i ,k,0)+ci2; + double tr2 = CC(i-1,k,0)+taur*cr2; // c add + double ti2 = CC(i ,k,0)+taur*ci2; + double tr3 = taui*(di2-di3); // t3 = taui*i*(d3-d2)? + double ti3 = taui*(dr3-dr2); + PM(CH(i-1,2,k),CH(ic-1,1,k),tr2,tr3) // PM(i) = t2+t3 + PM(CH(i ,2,k),CH(ic ,1,k),ti3,ti2) // PM(ic) = conj(t2-t3) + } + } + +NOINLINE static void radf4(size_t ido, size_t l1, const double * restrict cc, + double * restrict ch, const double * restrict wa) + { + const size_t cdim=4; + static const double hsqt2=0.70710678118654752440; + + for (size_t k=0; k<l1; k++) + { + double tr1,tr2; + PM (tr1,CH(0,2,k),CC(0,k,3),CC(0,k,1)) + PM (tr2,CH(ido-1,1,k),CC(0,k,0),CC(0,k,2)) + PM (CH(0,0,k),CH(ido-1,3,k),tr2,tr1) + } + if ((ido&1)==0) + for (size_t k=0; k<l1; k++) + { + double ti1=-hsqt2*(CC(ido-1,k,1)+CC(ido-1,k,3)); + double tr1= hsqt2*(CC(ido-1,k,1)-CC(ido-1,k,3)); + PM (CH(ido-1,0,k),CH(ido-1,2,k),CC(ido-1,k,0),tr1) + PM (CH( 0,3,k),CH( 0,1,k),ti1,CC(ido-1,k,2)) + } + if (ido<=2) return; + for (size_t k=0; k<l1; k++) + for (size_t i=2; i<ido; i+=2) + { + size_t ic=ido-i; + double ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4; + MULPM(cr2,ci2,WA(0,i-2),WA(0,i-1),CC(i-1,k,1),CC(i,k,1)) + MULPM(cr3,ci3,WA(1,i-2),WA(1,i-1),CC(i-1,k,2),CC(i,k,2)) + MULPM(cr4,ci4,WA(2,i-2),WA(2,i-1),CC(i-1,k,3),CC(i,k,3)) + PM(tr1,tr4,cr4,cr2) + PM(ti1,ti4,ci2,ci4) + PM(tr2,tr3,CC(i-1,k,0),cr3) + PM(ti2,ti3,CC(i ,k,0),ci3) + PM(CH(i-1,0,k),CH(ic-1,3,k),tr2,tr1) + PM(CH(i ,0,k),CH(ic ,3,k),ti1,ti2) + PM(CH(i-1,2,k),CH(ic-1,1,k),tr3,ti4) + PM(CH(i ,2,k),CH(ic ,1,k),tr4,ti3) + } + } + +NOINLINE static void radf5(size_t ido, size_t l1, const double * restrict cc, + double * restrict ch, const double * restrict wa) + { + const size_t cdim=5; + static const double tr11= 0.3090169943749474241, ti11=0.95105651629515357212, + tr12=-0.8090169943749474241, ti12=0.58778525229247312917; + + for (size_t k=0; k<l1; k++) + { + double cr2, cr3, ci4, ci5; + PM (cr2,ci5,CC(0,k,4),CC(0,k,1)) + PM (cr3,ci4,CC(0,k,3),CC(0,k,2)) + CH(0,0,k)=CC(0,k,0)+cr2+cr3; + CH(ido-1,1,k)=CC(0,k,0)+tr11*cr2+tr12*cr3; + CH(0,2,k)=ti11*ci5+ti12*ci4; + CH(ido-1,3,k)=CC(0,k,0)+tr12*cr2+tr11*cr3; + CH(0,4,k)=ti12*ci5-ti11*ci4; + } + if (ido==1) return; + for (size_t k=0; k<l1;++k) + for (size_t i=2; i<ido; i+=2) + { + double ci2, di2, ci4, ci5, di3, di4, di5, ci3, cr2, cr3, dr2, dr3, + dr4, dr5, cr5, cr4, ti2, ti3, ti5, ti4, tr2, tr3, tr4, tr5; + size_t ic=ido-i; + MULPM (dr2,di2,WA(0,i-2),WA(0,i-1),CC(i-1,k,1),CC(i,k,1)) + MULPM (dr3,di3,WA(1,i-2),WA(1,i-1),CC(i-1,k,2),CC(i,k,2)) + MULPM (dr4,di4,WA(2,i-2),WA(2,i-1),CC(i-1,k,3),CC(i,k,3)) + MULPM (dr5,di5,WA(3,i-2),WA(3,i-1),CC(i-1,k,4),CC(i,k,4)) + PM(cr2,ci5,dr5,dr2) + PM(ci2,cr5,di2,di5) + PM(cr3,ci4,dr4,dr3) + PM(ci3,cr4,di3,di4) + CH(i-1,0,k)=CC(i-1,k,0)+cr2+cr3; + CH(i ,0,k)=CC(i ,k,0)+ci2+ci3; + tr2=CC(i-1,k,0)+tr11*cr2+tr12*cr3; + ti2=CC(i ,k,0)+tr11*ci2+tr12*ci3; + tr3=CC(i-1,k,0)+tr12*cr2+tr11*cr3; + ti3=CC(i ,k,0)+tr12*ci2+tr11*ci3; + MULPM(tr5,tr4,cr5,cr4,ti11,ti12) + MULPM(ti5,ti4,ci5,ci4,ti11,ti12) + PM(CH(i-1,2,k),CH(ic-1,1,k),tr2,tr5) + PM(CH(i ,2,k),CH(ic ,1,k),ti5,ti2) + PM(CH(i-1,4,k),CH(ic-1,3,k),tr3,tr4) + PM(CH(i ,4,k),CH(ic ,3,k),ti4,ti3) + } + } + +#undef CC +#undef CH +#define C1(a,b,c) cc[(a)+ido*((b)+l1*(c))] +#define C2(a,b) cc[(a)+idl1*(b)] +#define CH2(a,b) ch[(a)+idl1*(b)] +#define CC(a,b,c) cc[(a)+ido*((b)+cdim*(c))] +#define CH(a,b,c) ch[(a)+ido*((b)+l1*(c))] +NOINLINE static void radfg(size_t ido, size_t ip, size_t l1, + double * restrict cc, double * restrict ch, const double * restrict wa, + const double * restrict csarr) + { + const size_t cdim=ip; + size_t ipph=(ip+1)/2; + size_t idl1 = ido*l1; + + if (ido>1) + { + for (size_t j=1, jc=ip-1; j<ipph; ++j,--jc) // 114 + { + size_t is=(j-1)*(ido-1), + is2=(jc-1)*(ido-1); + for (size_t k=0; k<l1; ++k) // 113 + { + size_t idij=is; + size_t idij2=is2; + for (size_t i=1; i<=ido-2; i+=2) // 112 + { + double t1=C1(i,k,j ), t2=C1(i+1,k,j ), + t3=C1(i,k,jc), t4=C1(i+1,k,jc); + double x1=wa[idij]*t1 + wa[idij+1]*t2, + x2=wa[idij]*t2 - wa[idij+1]*t1, + x3=wa[idij2]*t3 + wa[idij2+1]*t4, + x4=wa[idij2]*t4 - wa[idij2+1]*t3; + C1(i ,k,j ) = x1+x3; + C1(i ,k,jc) = x2-x4; + C1(i+1,k,j ) = x2+x4; + C1(i+1,k,jc) = x3-x1; + idij+=2; + idij2+=2; + } + } + } + } + + for (size_t j=1, jc=ip-1; j<ipph; ++j,--jc) // 123 + for (size_t k=0; k<l1; ++k) // 122 + { + double t1=C1(0,k,j), t2=C1(0,k,jc); + C1(0,k,j ) = t1+t2; + C1(0,k,jc) = t2-t1; + } + +//everything in C +//memset(ch,0,ip*l1*ido*sizeof(double)); + + for (size_t l=1,lc=ip-1; l<ipph; ++l,--lc) // 127 + { + for (size_t ik=0; ik<idl1; ++ik) // 124 + { + CH2(ik,l ) = C2(ik,0)+csarr[2*l]*C2(ik,1)+csarr[4*l]*C2(ik,2); + CH2(ik,lc) = csarr[2*l+1]*C2(ik,ip-1)+csarr[4*l+1]*C2(ik,ip-2); + } + size_t iang = 2*l; + size_t j=3, jc=ip-3; + for (; j<ipph-3; j+=4,jc-=4) // 126 + { + iang+=l; if (iang>=ip) iang-=ip; + double ar1=csarr[2*iang], ai1=csarr[2*iang+1]; + iang+=l; if (iang>=ip) iang-=ip; + double ar2=csarr[2*iang], ai2=csarr[2*iang+1]; + iang+=l; if (iang>=ip) iang-=ip; + double ar3=csarr[2*iang], ai3=csarr[2*iang+1]; + iang+=l; if (iang>=ip) iang-=ip; + double ar4=csarr[2*iang], ai4=csarr[2*iang+1]; + for (size_t ik=0; ik<idl1; ++ik) // 125 + { + CH2(ik,l ) += ar1*C2(ik,j )+ar2*C2(ik,j +1) + +ar3*C2(ik,j +2)+ar4*C2(ik,j +3); + CH2(ik,lc) += ai1*C2(ik,jc)+ai2*C2(ik,jc-1) + +ai3*C2(ik,jc-2)+ai4*C2(ik,jc-3); + } + } + for (; j<ipph-1; j+=2,jc-=2) // 126 + { + iang+=l; if (iang>=ip) iang-=ip; + double ar1=csarr[2*iang], ai1=csarr[2*iang+1]; + iang+=l; if (iang>=ip) iang-=ip; + double ar2=csarr[2*iang], ai2=csarr[2*iang+1]; + for (size_t ik=0; ik<idl1; ++ik) // 125 + { + CH2(ik,l ) += ar1*C2(ik,j )+ar2*C2(ik,j +1); + CH2(ik,lc) += ai1*C2(ik,jc)+ai2*C2(ik,jc-1); + } + } + for (; j<ipph; ++j,--jc) // 126 + { + iang+=l; if (iang>=ip) iang-=ip; + double ar=csarr[2*iang], ai=csarr[2*iang+1]; + for (size_t ik=0; ik<idl1; ++ik) // 125 + { + CH2(ik,l ) += ar*C2(ik,j ); + CH2(ik,lc) += ai*C2(ik,jc); + } + } + } + for (size_t ik=0; ik<idl1; ++ik) // 101 + CH2(ik,0) = C2(ik,0); + for (size_t j=1; j<ipph; ++j) // 129 + for (size_t ik=0; ik<idl1; ++ik) // 128 + CH2(ik,0) += C2(ik,j); + +// everything in CH at this point! +//memset(cc,0,ip*l1*ido*sizeof(double)); + + for (size_t k=0; k<l1; ++k) // 131 + for (size_t i=0; i<ido; ++i) // 130 + CC(i,0,k) = CH(i,k,0); + + for (size_t j=1, jc=ip-1; j<ipph; ++j,--jc) // 137 + { + size_t j2=2*j-1; + for (size_t k=0; k<l1; ++k) // 136 + { + CC(ido-1,j2,k) = CH(0,k,j); + CC(0,j2+1,k) = CH(0,k,jc); + } + } + + if (ido==1) return; + + for (size_t j=1, jc=ip-1; j<ipph; ++j,--jc) // 140 + { + size_t j2=2*j-1; + for(size_t k=0; k<l1; ++k) // 139 + for(size_t i=1, ic=ido-i-2; i<=ido-2; i+=2, ic-=2) // 138 + { + CC(i ,j2+1,k) = CH(i ,k,j )+CH(i ,k,jc); + CC(ic ,j2 ,k) = CH(i ,k,j )-CH(i ,k,jc); + CC(i+1 ,j2+1,k) = CH(i+1,k,j )+CH(i+1,k,jc); + CC(ic+1,j2 ,k) = CH(i+1,k,jc)-CH(i+1,k,j ); + } + } + } +#undef C1 +#undef C2 +#undef CH2 + +#undef CH +#undef CC +#define CH(a,b,c) ch[(a)+ido*((b)+l1*(c))] +#define CC(a,b,c) cc[(a)+ido*((b)+cdim*(c))] + +NOINLINE static void radb2(size_t ido, size_t l1, const double * restrict cc, + double * restrict ch, const double * restrict wa) + { + const size_t cdim=2; + + for (size_t k=0; k<l1; k++) + PM (CH(0,k,0),CH(0,k,1),CC(0,0,k),CC(ido-1,1,k)) + if ((ido&1)==0) + for (size_t k=0; k<l1; k++) + { + CH(ido-1,k,0) = 2.*CC(ido-1,0,k); + CH(ido-1,k,1) =-2.*CC(0 ,1,k); + } + if (ido<=2) return; + for (size_t k=0; k<l1;++k) + for (size_t i=2; i<ido; i+=2) + { + size_t ic=ido-i; + double ti2, tr2; + PM (CH(i-1,k,0),tr2,CC(i-1,0,k),CC(ic-1,1,k)) + PM (ti2,CH(i ,k,0),CC(i ,0,k),CC(ic ,1,k)) + MULPM (CH(i,k,1),CH(i-1,k,1),WA(0,i-2),WA(0,i-1),ti2,tr2) + } + } + +NOINLINE static void radb3(size_t ido, size_t l1, const double * restrict cc, + double * restrict ch, const double * restrict wa) + { + const size_t cdim=3; + static const double taur=-0.5, taui=0.86602540378443864676; + + for (size_t k=0; k<l1; k++) + { + double tr2=2.*CC(ido-1,1,k); + double cr2=CC(0,0,k)+taur*tr2; + CH(0,k,0)=CC(0,0,k)+tr2; + double ci3=2.*taui*CC(0,2,k); + PM (CH(0,k,2),CH(0,k,1),cr2,ci3); + } + if (ido==1) return; + for (size_t k=0; k<l1; k++) + for (size_t i=2; i<ido; i+=2) + { + size_t ic=ido-i; + double tr2=CC(i-1,2,k)+CC(ic-1,1,k); // t2=CC(I) + conj(CC(ic)) + double ti2=CC(i ,2,k)-CC(ic ,1,k); + double cr2=CC(i-1,0,k)+taur*tr2; // c2=CC +taur*t2 + double ci2=CC(i ,0,k)+taur*ti2; + CH(i-1,k,0)=CC(i-1,0,k)+tr2; // CH=CC+t2 + CH(i ,k,0)=CC(i ,0,k)+ti2; + double cr3=taui*(CC(i-1,2,k)-CC(ic-1,1,k));// c3=taui*(CC(i)-conj(CC(ic))) + double ci3=taui*(CC(i ,2,k)+CC(ic ,1,k)); + double di2, di3, dr2, dr3; + PM(dr3,dr2,cr2,ci3) // d2= (cr2-ci3, ci2+cr3) = c2+i*c3 + PM(di2,di3,ci2,cr3) // d3= (cr2+ci3, ci2-cr3) = c2-i*c3 + MULPM(CH(i,k,1),CH(i-1,k,1),WA(0,i-2),WA(0,i-1),di2,dr2) // ch = WA*d2 + MULPM(CH(i,k,2),CH(i-1,k,2),WA(1,i-2),WA(1,i-1),di3,dr3) + } + } + +NOINLINE static void radb4(size_t ido, size_t l1, const double * restrict cc, + double * restrict ch, const double * restrict wa) + { + const size_t cdim=4; + static const double sqrt2=1.41421356237309504880; + + for (size_t k=0; k<l1; k++) + { + double tr1, tr2; + PM (tr2,tr1,CC(0,0,k),CC(ido-1,3,k)) + double tr3=2.*CC(ido-1,1,k); + double tr4=2.*CC(0,2,k); + PM (CH(0,k,0),CH(0,k,2),tr2,tr3) + PM (CH(0,k,3),CH(0,k,1),tr1,tr4) + } + if ((ido&1)==0) + for (size_t k=0; k<l1; k++) + { + double tr1,tr2,ti1,ti2; + PM (ti1,ti2,CC(0 ,3,k),CC(0 ,1,k)) + PM (tr2,tr1,CC(ido-1,0,k),CC(ido-1,2,k)) + CH(ido-1,k,0)=tr2+tr2; + CH(ido-1,k,1)=sqrt2*(tr1-ti1); + CH(ido-1,k,2)=ti2+ti2; + CH(ido-1,k,3)=-sqrt2*(tr1+ti1); + } + if (ido<=2) return; + for (size_t k=0; k<l1;++k) + for (size_t i=2; i<ido; i+=2) + { + double ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4; + size_t ic=ido-i; + PM (tr2,tr1,CC(i-1,0,k),CC(ic-1,3,k)) + PM (ti1,ti2,CC(i ,0,k),CC(ic ,3,k)) + PM (tr4,ti3,CC(i ,2,k),CC(ic ,1,k)) + PM (tr3,ti4,CC(i-1,2,k),CC(ic-1,1,k)) + PM (CH(i-1,k,0),cr3,tr2,tr3) + PM (CH(i ,k,0),ci3,ti2,ti3) + PM (cr4,cr2,tr1,tr4) + PM (ci2,ci4,ti1,ti4) + MULPM (CH(i,k,1),CH(i-1,k,1),WA(0,i-2),WA(0,i-1),ci2,cr2) + MULPM (CH(i,k,2),CH(i-1,k,2),WA(1,i-2),WA(1,i-1),ci3,cr3) + MULPM (CH(i,k,3),CH(i-1,k,3),WA(2,i-2),WA(2,i-1),ci4,cr4) + } + } + +NOINLINE static void radb5(size_t ido, size_t l1, const double * restrict cc, + double * restrict ch, const double * restrict wa) + { + const size_t cdim=5; + static const double tr11= 0.3090169943749474241, ti11=0.95105651629515357212, + tr12=-0.8090169943749474241, ti12=0.58778525229247312917; + + for (size_t k=0; k<l1; k++) + { + double ti5=CC(0,2,k)+CC(0,2,k); + double ti4=CC(0,4,k)+CC(0,4,k); + double tr2=CC(ido-1,1,k)+CC(ido-1,1,k); + double tr3=CC(ido-1,3,k)+CC(ido-1,3,k); + CH(0,k,0)=CC(0,0,k)+tr2+tr3; + double cr2=CC(0,0,k)+tr11*tr2+tr12*tr3; + double cr3=CC(0,0,k)+tr12*tr2+tr11*tr3; + double ci4, ci5; + MULPM(ci5,ci4,ti5,ti4,ti11,ti12) + PM(CH(0,k,4),CH(0,k,1),cr2,ci5) + PM(CH(0,k,3),CH(0,k,2),cr3,ci4) + } + if (ido==1) return; + for (size_t k=0; k<l1;++k) + for (size_t i=2; i<ido; i+=2) + { + size_t ic=ido-i; + double tr2, tr3, tr4, tr5, ti2, ti3, ti4, ti5; + PM(tr2,tr5,CC(i-1,2,k),CC(ic-1,1,k)) + PM(ti5,ti2,CC(i ,2,k),CC(ic ,1,k)) + PM(tr3,tr4,CC(i-1,4,k),CC(ic-1,3,k)) + PM(ti4,ti3,CC(i ,4,k),CC(ic ,3,k)) + CH(i-1,k,0)=CC(i-1,0,k)+tr2+tr3; + CH(i ,k,0)=CC(i ,0,k)+ti2+ti3; + double cr2=CC(i-1,0,k)+tr11*tr2+tr12*tr3; + double ci2=CC(i ,0,k)+tr11*ti2+tr12*ti3; + double cr3=CC(i-1,0,k)+tr12*tr2+tr11*tr3; + double ci3=CC(i ,0,k)+tr12*ti2+tr11*ti3; + double ci4, ci5, cr5, cr4; + MULPM(cr5,cr4,tr5,tr4,ti11,ti12) + MULPM(ci5,ci4,ti5,ti4,ti11,ti12) + double dr2, dr3, dr4, dr5, di2, di3, di4, di5; + PM(dr4,dr3,cr3,ci4) + PM(di3,di4,ci3,cr4) + PM(dr5,dr2,cr2,ci5) + PM(di2,di5,ci2,cr5) + MULPM(CH(i,k,1),CH(i-1,k,1),WA(0,i-2),WA(0,i-1),di2,dr2) + MULPM(CH(i,k,2),CH(i-1,k,2),WA(1,i-2),WA(1,i-1),di3,dr3) + MULPM(CH(i,k,3),CH(i-1,k,3),WA(2,i-2),WA(2,i-1),di4,dr4) + MULPM(CH(i,k,4),CH(i-1,k,4),WA(3,i-2),WA(3,i-1),di5,dr5) + } + } + +#undef CC +#undef CH +#define CC(a,b,c) cc[(a)+ido*((b)+cdim*(c))] +#define CH(a,b,c) ch[(a)+ido*((b)+l1*(c))] +#define C1(a,b,c) cc[(a)+ido*((b)+l1*(c))] +#define C2(a,b) cc[(a)+idl1*(b)] +#define CH2(a,b) ch[(a)+idl1*(b)] + +NOINLINE static void radbg(size_t ido, size_t ip, size_t l1, + double * restrict cc, double * restrict ch, const double * restrict wa, + const double * restrict csarr) + { + const size_t cdim=ip; + size_t ipph=(ip+1)/ 2; + size_t idl1 = ido*l1; + + for (size_t k=0; k<l1; ++k) // 102 + for (size_t i=0; i<ido; ++i) // 101 + CH(i,k,0) = CC(i,0,k); + for (size_t j=1, jc=ip-1; j<ipph; ++j, --jc) // 108 + { + size_t j2=2*j-1; + for (size_t k=0; k<l1; ++k) + { + CH(0,k,j ) = 2*CC(ido-1,j2,k); + CH(0,k,jc) = 2*CC(0,j2+1,k); + } + } + + if (ido!=1) + { + for (size_t j=1, jc=ip-1; j<ipph; ++j,--jc) // 111 + { + size_t j2=2*j-1; + for (size_t k=0; k<l1; ++k) + for (size_t i=1, ic=ido-i-2; i<=ido-2; i+=2, ic-=2) // 109 + { + CH(i ,k,j ) = CC(i ,j2+1,k)+CC(ic ,j2,k); + CH(i ,k,jc) = CC(i ,j2+1,k)-CC(ic ,j2,k); + CH(i+1,k,j ) = CC(i+1,j2+1,k)-CC(ic+1,j2,k); + CH(i+1,k,jc) = CC(i+1,j2+1,k)+CC(ic+1,j2,k); + } + } + } + for (size_t l=1,lc=ip-1; l<ipph; ++l,--lc) + { + for (size_t ik=0; ik<idl1; ++ik) + { + C2(ik,l ) = CH2(ik,0)+csarr[2*l]*CH2(ik,1)+csarr[4*l]*CH2(ik,2); + C2(ik,lc) = csarr[2*l+1]*CH2(ik,ip-1)+csarr[4*l+1]*CH2(ik,ip-2); + } + size_t iang=2*l; + size_t j=3,jc=ip-3; + for(; j<ipph-3; j+=4,jc-=4) + { + iang+=l; if(iang>ip) iang-=ip; + double ar1=csarr[2*iang], ai1=csarr[2*iang+1]; + iang+=l; if(iang>ip) iang-=ip; + double ar2=csarr[2*iang], ai2=csarr[2*iang+1]; + iang+=l; if(iang>ip) iang-=ip; + double ar3=csarr[2*iang], ai3=csarr[2*iang+1]; + iang+=l; if(iang>ip) iang-=ip; + double ar4=csarr[2*iang], ai4=csarr[2*iang+1]; + for (size_t ik=0; ik<idl1; ++ik) + { + C2(ik,l ) += ar1*CH2(ik,j )+ar2*CH2(ik,j +1) + +ar3*CH2(ik,j +2)+ar4*CH2(ik,j +3); + C2(ik,lc) += ai1*CH2(ik,jc)+ai2*CH2(ik,jc-1) + +ai3*CH2(ik,jc-2)+ai4*CH2(ik,jc-3); + } + } + for(; j<ipph-1; j+=2,jc-=2) + { + iang+=l; if(iang>ip) iang-=ip; + double ar1=csarr[2*iang], ai1=csarr[2*iang+1]; + iang+=l; if(iang>ip) iang-=ip; + double ar2=csarr[2*iang], ai2=csarr[2*iang+1]; + for (size_t ik=0; ik<idl1; ++ik) + { + C2(ik,l ) += ar1*CH2(ik,j )+ar2*CH2(ik,j +1); + C2(ik,lc) += ai1*CH2(ik,jc)+ai2*CH2(ik,jc-1); + } + } + for(; j<ipph; ++j,--jc) + { + iang+=l; if(iang>ip) iang-=ip; + double war=csarr[2*iang], wai=csarr[2*iang+1]; + for (size_t ik=0; ik<idl1; ++ik) + { + C2(ik,l ) += war*CH2(ik,j ); + C2(ik,lc) += wai*CH2(ik,jc); + } + } + } + for (size_t j=1; j<ipph; ++j) + for (size_t ik=0; ik<idl1; ++ik) + CH2(ik,0) += CH2(ik,j); + for (size_t j=1, jc=ip-1; j<ipph; ++j,--jc) // 124 + for (size_t k=0; k<l1; ++k) + { + CH(0,k,j ) = C1(0,k,j)-C1(0,k,jc); + CH(0,k,jc) = C1(0,k,j)+C1(0,k,jc); + } + + if (ido==1) return; + + for (size_t j=1, jc=ip-1; j<ipph; ++j, --jc) // 127 + for (size_t k=0; k<l1; ++k) + for (size_t i=1; i<=ido-2; i+=2) + { + CH(i ,k,j ) = C1(i ,k,j)-C1(i+1,k,jc); + CH(i ,k,jc) = C1(i ,k,j)+C1(i+1,k,jc); + CH(i+1,k,j ) = C1(i+1,k,j)+C1(i ,k,jc); + CH(i+1,k,jc) = C1(i+1,k,j)-C1(i ,k,jc); + } + +// All in CH + + for (size_t j=1; j<ip; ++j) + { + size_t is = (j-1)*(ido-1); + for (size_t k=0; k<l1; ++k) + { + size_t idij = is; + for (size_t i=1; i<=ido-2; i+=2) + { + double t1=CH(i,k,j), t2=CH(i+1,k,j); + CH(i ,k,j) = wa[idij]*t1-wa[idij+1]*t2; + CH(i+1,k,j) = wa[idij]*t2+wa[idij+1]*t1; + idij+=2; + } + } + } + } +#undef C1 +#undef C2 +#undef CH2 + +#undef CC +#undef CH +#undef PM +#undef MULPM +#undef WA + +static void copy_and_norm(double *c, double *p1, size_t n, double fct) + { + if (p1!=c) + { + if (fct!=1.) + for (size_t i=0; i<n; ++i) + c[i] = fct*p1[i]; + else + memcpy (c,p1,n*sizeof(double)); + } + else + if (fct!=1.) + for (size_t i=0; i<n; ++i) + c[i] *= fct; + } + +WARN_UNUSED_RESULT +static int rfftp_forward(rfftp_plan plan, double c[], double fct) + { + if (plan->length==1) return 0; + size_t n=plan->length; + size_t l1=n, nf=plan->nfct; + double *ch = RALLOC(double, n); + if (!ch) return -1; + double *p1=c, *p2=ch; + + for(size_t k1=0; k1<nf;++k1) + { + size_t k=nf-k1-1; + size_t ip=plan->fct[k].fct; + size_t ido=n / l1; + l1 /= ip; + if(ip==4) + radf4(ido, l1, p1, p2, plan->fct[k].tw); + else if(ip==2) + radf2(ido, l1, p1, p2, plan->fct[k].tw); + else if(ip==3) + radf3(ido, l1, p1, p2, plan->fct[k].tw); + else if(ip==5) + radf5(ido, l1, p1, p2, plan->fct[k].tw); + else + { + radfg(ido, ip, l1, p1, p2, plan->fct[k].tw, plan->fct[k].tws); + SWAP (p1,p2,double *); + } + SWAP (p1,p2,double *); + } + copy_and_norm(c,p1,n,fct); + DEALLOC(ch); + return 0; + } + +WARN_UNUSED_RESULT +static int rfftp_backward(rfftp_plan plan, double c[], double fct) + { + if (plan->length==1) return 0; + size_t n=plan->length; + size_t l1=1, nf=plan->nfct; + double *ch = RALLOC(double, n); + if (!ch) return -1; + double *p1=c, *p2=ch; + + for(size_t k=0; k<nf; k++) + { + size_t ip = plan->fct[k].fct, + ido= n/(ip*l1); + if(ip==4) + radb4(ido, l1, p1, p2, plan->fct[k].tw); + else if(ip==2) + radb2(ido, l1, p1, p2, plan->fct[k].tw); + else if(ip==3) + radb3(ido, l1, p1, p2, plan->fct[k].tw); + else if(ip==5) + radb5(ido, l1, p1, p2, plan->fct[k].tw); + else + radbg(ido, ip, l1, p1, p2, plan->fct[k].tw, plan->fct[k].tws); + SWAP (p1,p2,double *); + l1*=ip; + } + copy_and_norm(c,p1,n,fct); + DEALLOC(ch); + return 0; + } + +WARN_UNUSED_RESULT +static int rfftp_factorize (rfftp_plan plan) + { + size_t length=plan->length; + size_t nfct=0; + while ((length%4)==0) + { if (nfct>=NFCT) return -1; plan->fct[nfct++].fct=4; length>>=2; } + if ((length%2)==0) + { + length>>=1; + // factor 2 should be at the front of the factor list + if (nfct>=NFCT) return -1; + plan->fct[nfct++].fct=2; + SWAP(plan->fct[0].fct, plan->fct[nfct-1].fct,size_t); + } + size_t maxl=(size_t)(sqrt((double)length))+1; + for (size_t divisor=3; (length>1)&&(divisor<maxl); divisor+=2) + if ((length%divisor)==0) + { + while ((length%divisor)==0) + { + if (nfct>=NFCT) return -1; + plan->fct[nfct++].fct=divisor; + length/=divisor; + } + maxl=(size_t)(sqrt((double)length))+1; + } + if (length>1) plan->fct[nfct++].fct=length; + plan->nfct=nfct; + return 0; + } + +static size_t rfftp_twsize(rfftp_plan plan) + { + size_t twsize=0, l1=1; + for (size_t k=0; k<plan->nfct; ++k) + { + size_t ip=plan->fct[k].fct, ido= plan->length/(l1*ip); + twsize+=(ip-1)*(ido-1); + if (ip>5) twsize+=2*ip; + l1*=ip; + } + return twsize; + return 0; + } + +WARN_UNUSED_RESULT NOINLINE static int rfftp_comp_twiddle (rfftp_plan plan) + { + size_t length=plan->length; + double *twid = RALLOC(double, 2*length); + if (!twid) return -1; + sincos_2pibyn_half(length, twid); + size_t l1=1; + double *ptr=plan->mem; + for (size_t k=0; k<plan->nfct; ++k) + { + size_t ip=plan->fct[k].fct, ido=length/(l1*ip); + if (k<plan->nfct-1) // last factor doesn't need twiddles + { + plan->fct[k].tw=ptr; ptr+=(ip-1)*(ido-1); + for (size_t j=1; j<ip; ++j) + for (size_t i=1; i<=(ido-1)/2; ++i) + { + plan->fct[k].tw[(j-1)*(ido-1)+2*i-2] = twid[2*j*l1*i]; + plan->fct[k].tw[(j-1)*(ido-1)+2*i-1] = twid[2*j*l1*i+1]; + } + } + if (ip>5) // special factors required by *g functions + { + plan->fct[k].tws=ptr; ptr+=2*ip; + plan->fct[k].tws[0] = 1.; + plan->fct[k].tws[1] = 0.; + for (size_t i=1; i<=(ip>>1); ++i) + { + plan->fct[k].tws[2*i ] = twid[2*i*(length/ip)]; + plan->fct[k].tws[2*i+1] = twid[2*i*(length/ip)+1]; + plan->fct[k].tws[2*(ip-i) ] = twid[2*i*(length/ip)]; + plan->fct[k].tws[2*(ip-i)+1] = -twid[2*i*(length/ip)+1]; + } + } + l1*=ip; + } + DEALLOC(twid); + return 0; + } + +NOINLINE static rfftp_plan make_rfftp_plan (size_t length) + { + if (length==0) return NULL; + rfftp_plan plan = RALLOC(rfftp_plan_i,1); + if (!plan) return NULL; + plan->length=length; + plan->nfct=0; + plan->mem=NULL; + for (size_t i=0; i<NFCT; ++i) + plan->fct[i]=(rfftp_fctdata){0,0,0}; + if (length==1) return plan; + if (rfftp_factorize(plan)!=0) { DEALLOC(plan); return NULL; } + size_t tws=rfftp_twsize(plan); + plan->mem=RALLOC(double,tws); + if (!plan->mem) { DEALLOC(plan); return NULL; } + if (rfftp_comp_twiddle(plan)!=0) + { DEALLOC(plan->mem); DEALLOC(plan); return NULL; } + return plan; + } + +NOINLINE static void destroy_rfftp_plan (rfftp_plan plan) + { + DEALLOC(plan->mem); + DEALLOC(plan); + } + +typedef struct fftblue_plan_i + { + size_t n, n2; + cfftp_plan plan; + double *mem; + double *bk, *bkf; + } fftblue_plan_i; +typedef struct fftblue_plan_i * fftblue_plan; + +NOINLINE static fftblue_plan make_fftblue_plan (size_t length) + { + fftblue_plan plan = RALLOC(fftblue_plan_i,1); + if (!plan) return NULL; + plan->n = length; + plan->n2 = good_size(plan->n*2-1); + plan->mem = RALLOC(double, 2*plan->n+2*plan->n2); + if (!plan->mem) { DEALLOC(plan); return NULL; } + plan->bk = plan->mem; + plan->bkf = plan->bk+2*plan->n; + +/* initialize b_k */ + double *tmp = RALLOC(double,4*plan->n); + if (!tmp) { DEALLOC(plan->mem); DEALLOC(plan); return NULL; } + sincos_2pibyn(2*plan->n,tmp); + plan->bk[0] = 1; + plan->bk[1] = 0; + + size_t coeff=0; + for (size_t m=1; m<plan->n; ++m) + { + coeff+=2*m-1; + if (coeff>=2*plan->n) coeff-=2*plan->n; + plan->bk[2*m ] = tmp[2*coeff ]; + plan->bk[2*m+1] = tmp[2*coeff+1]; + } + + /* initialize the zero-padded, Fourier transformed b_k. Add normalisation. */ + double xn2 = 1./plan->n2; + plan->bkf[0] = plan->bk[0]*xn2; + plan->bkf[1] = plan->bk[1]*xn2; + for (size_t m=2; m<2*plan->n; m+=2) + { + plan->bkf[m] = plan->bkf[2*plan->n2-m] = plan->bk[m] *xn2; + plan->bkf[m+1] = plan->bkf[2*plan->n2-m+1] = plan->bk[m+1] *xn2; + } + for (size_t m=2*plan->n;m<=(2*plan->n2-2*plan->n+1);++m) + plan->bkf[m]=0.; + plan->plan=make_cfftp_plan(plan->n2); + if (!plan->plan) + { DEALLOC(tmp); DEALLOC(plan->mem); DEALLOC(plan); return NULL; } + if (cfftp_forward(plan->plan,plan->bkf,1.)!=0) + { DEALLOC(tmp); DEALLOC(plan->mem); DEALLOC(plan); return NULL; } + DEALLOC(tmp); + + return plan; + } + +NOINLINE static void destroy_fftblue_plan (fftblue_plan plan) + { + DEALLOC(plan->mem); + destroy_cfftp_plan(plan->plan); + DEALLOC(plan); + } + +NOINLINE WARN_UNUSED_RESULT +static int fftblue_fft(fftblue_plan plan, double c[], int isign, double fct) + { + size_t n=plan->n; + size_t n2=plan->n2; + double *bk = plan->bk; + double *bkf = plan->bkf; + double *akf = RALLOC(double, 2*n2); + if (!akf) return -1; + +/* initialize a_k and FFT it */ + if (isign>0) + for (size_t m=0; m<2*n; m+=2) + { + akf[m] = c[m]*bk[m] - c[m+1]*bk[m+1]; + akf[m+1] = c[m]*bk[m+1] + c[m+1]*bk[m]; + } + else + for (size_t m=0; m<2*n; m+=2) + { + akf[m] = c[m]*bk[m] + c[m+1]*bk[m+1]; + akf[m+1] =-c[m]*bk[m+1] + c[m+1]*bk[m]; + } + for (size_t m=2*n; m<2*n2; ++m) + akf[m]=0; + + if (cfftp_forward (plan->plan,akf,fct)!=0) + { DEALLOC(akf); return -1; } + +/* do the convolution */ + if (isign>0) + for (size_t m=0; m<2*n2; m+=2) + { + double im = -akf[m]*bkf[m+1] + akf[m+1]*bkf[m]; + akf[m ] = akf[m]*bkf[m] + akf[m+1]*bkf[m+1]; + akf[m+1] = im; + } + else + for (size_t m=0; m<2*n2; m+=2) + { + double im = akf[m]*bkf[m+1] + akf[m+1]*bkf[m]; + akf[m ] = akf[m]*bkf[m] - akf[m+1]*bkf[m+1]; + akf[m+1] = im; + } + +/* inverse FFT */ + if (cfftp_backward (plan->plan,akf,1.)!=0) + { DEALLOC(akf); return -1; } + +/* multiply by b_k */ + if (isign>0) + for (size_t m=0; m<2*n; m+=2) + { + c[m] = bk[m] *akf[m] - bk[m+1]*akf[m+1]; + c[m+1] = bk[m+1]*akf[m] + bk[m] *akf[m+1]; + } + else + for (size_t m=0; m<2*n; m+=2) + { + c[m] = bk[m] *akf[m] + bk[m+1]*akf[m+1]; + c[m+1] =-bk[m+1]*akf[m] + bk[m] *akf[m+1]; + } + DEALLOC(akf); + return 0; + } + +WARN_UNUSED_RESULT +static int cfftblue_backward(fftblue_plan plan, double c[], double fct) + { return fftblue_fft(plan,c,1,fct); } + +WARN_UNUSED_RESULT +static int cfftblue_forward(fftblue_plan plan, double c[], double fct) + { return fftblue_fft(plan,c,-1,fct); } + +WARN_UNUSED_RESULT +static int rfftblue_backward(fftblue_plan plan, double c[], double fct) + { + size_t n=plan->n; + double *tmp = RALLOC(double,2*n); + if (!tmp) return -1; + tmp[0]=c[0]; + tmp[1]=0.; + memcpy (tmp+2,c+1, (n-1)*sizeof(double)); + if ((n&1)==0) tmp[n+1]=0.; + for (size_t m=2; m<n; m+=2) + { + tmp[2*n-m]=tmp[m]; + tmp[2*n-m+1]=-tmp[m+1]; + } + if (fftblue_fft(plan,tmp,1,fct)!=0) + { DEALLOC(tmp); return -1; } + for (size_t m=0; m<n; ++m) + c[m] = tmp[2*m]; + DEALLOC(tmp); + return 0; + } + +WARN_UNUSED_RESULT +static int rfftblue_forward(fftblue_plan plan, double c[], double fct) + { + size_t n=plan->n; + double *tmp = RALLOC(double,2*n); + if (!tmp) return -1; + for (size_t m=0; m<n; ++m) + { + tmp[2*m] = c[m]; + tmp[2*m+1] = 0.; + } + if (fftblue_fft(plan,tmp,-1,fct)!=0) + { DEALLOC(tmp); return -1; } + c[0] = tmp[0]; + memcpy (c+1, tmp+2, (n-1)*sizeof(double)); + DEALLOC(tmp); + return 0; + } + +typedef struct cfft_plan_i + { + cfftp_plan packplan; + fftblue_plan blueplan; + } cfft_plan_i; + +static cfft_plan make_cfft_plan (size_t length) + { + if (length==0) return NULL; + cfft_plan plan = RALLOC(cfft_plan_i,1); + if (!plan) return NULL; + plan->blueplan=0; + plan->packplan=0; + if ((length<50) || (largest_prime_factor(length)<=sqrt(length))) + { + plan->packplan=make_cfftp_plan(length); + if (!plan->packplan) { DEALLOC(plan); return NULL; } + return plan; + } + double comp1 = cost_guess(length); + double comp2 = 2*cost_guess(good_size(2*length-1)); + comp2*=1.5; /* fudge factor that appears to give good overall performance */ + if (comp2<comp1) // use Bluestein + { + plan->blueplan=make_fftblue_plan(length); + if (!plan->blueplan) { DEALLOC(plan); return NULL; } + } + else + { + plan->packplan=make_cfftp_plan(length); + if (!plan->packplan) { DEALLOC(plan); return NULL; } + } + return plan; + } + +static void destroy_cfft_plan (cfft_plan plan) + { + if (plan->blueplan) + destroy_fftblue_plan(plan->blueplan); + if (plan->packplan) + destroy_cfftp_plan(plan->packplan); + DEALLOC(plan); + } + +WARN_UNUSED_RESULT static int cfft_backward(cfft_plan plan, double c[], double fct) + { + if (plan->packplan) + return cfftp_backward(plan->packplan,c,fct); + // if (plan->blueplan) + return cfftblue_backward(plan->blueplan,c,fct); + } + +WARN_UNUSED_RESULT static int cfft_forward(cfft_plan plan, double c[], double fct) + { + if (plan->packplan) + return cfftp_forward(plan->packplan,c,fct); + // if (plan->blueplan) + return cfftblue_forward(plan->blueplan,c,fct); + } + +typedef struct rfft_plan_i + { + rfftp_plan packplan; + fftblue_plan blueplan; + } rfft_plan_i; + +static rfft_plan make_rfft_plan (size_t length) + { + if (length==0) return NULL; + rfft_plan plan = RALLOC(rfft_plan_i,1); + if (!plan) return NULL; + plan->blueplan=0; + plan->packplan=0; + if ((length<50) || (largest_prime_factor(length)<=sqrt(length))) + { + plan->packplan=make_rfftp_plan(length); + if (!plan->packplan) { DEALLOC(plan); return NULL; } + return plan; + } + double comp1 = 0.5*cost_guess(length); + double comp2 = 2*cost_guess(good_size(2*length-1)); + comp2*=1.5; /* fudge factor that appears to give good overall performance */ + if (comp2<comp1) // use Bluestein + { + plan->blueplan=make_fftblue_plan(length); + if (!plan->blueplan) { DEALLOC(plan); return NULL; } + } + else + { + plan->packplan=make_rfftp_plan(length); + if (!plan->packplan) { DEALLOC(plan); return NULL; } + } + return plan; + } + +static void destroy_rfft_plan (rfft_plan plan) + { + if (plan->blueplan) + destroy_fftblue_plan(plan->blueplan); + if (plan->packplan) + destroy_rfftp_plan(plan->packplan); + DEALLOC(plan); + } + +WARN_UNUSED_RESULT static int rfft_backward(rfft_plan plan, double c[], double fct) + { + if (plan->packplan) + return rfftp_backward(plan->packplan,c,fct); + else // if (plan->blueplan) + return rfftblue_backward(plan->blueplan,c,fct); + } + +WARN_UNUSED_RESULT static int rfft_forward(rfft_plan plan, double c[], double fct) + { + if (plan->packplan) + return rfftp_forward(plan->packplan,c,fct); + else // if (plan->blueplan) + return rfftblue_forward(plan->blueplan,c,fct); + } + +#define NPY_NO_DEPRECATED_API NPY_API_VERSION + +#include "Python.h" +#include "numpy/arrayobject.h" + +static PyObject * +execute_complex(PyObject *a1, int is_forward, double fct) +{ + PyArrayObject *data = (PyArrayObject *)PyArray_CopyFromObject(a1, NPY_CDOUBLE, 1, 0); + if (!data) return NULL; + + int npts = PyArray_DIM(data, PyArray_NDIM(data) - 1); + cfft_plan plan=NULL; + + int nrepeats = PyArray_SIZE(data)/npts; + double *dptr = (double *)PyArray_DATA(data); + int fail=0; + Py_BEGIN_ALLOW_THREADS; + NPY_SIGINT_ON; + plan = make_cfft_plan(npts); + if (!plan) fail=1; + if (!fail) + for (int i = 0; i < nrepeats; i++) { + int res = is_forward ? + cfft_forward(plan, dptr, fct) : cfft_backward(plan, dptr, fct); + if (res!=0) { fail=1; break; } + dptr += npts*2; + } + if (plan) destroy_cfft_plan(plan); + NPY_SIGINT_OFF; + Py_END_ALLOW_THREADS; + if (fail) { + Py_XDECREF(data); + return PyErr_NoMemory(); + } + return (PyObject *)data; +} + +static PyObject * +execute_real_forward(PyObject *a1, double fct) +{ + rfft_plan plan=NULL; + int fail = 0; + PyArrayObject *data = (PyArrayObject *)PyArray_ContiguousFromObject(a1, + NPY_DOUBLE, 1, 0); + if (!data) return NULL; + + int ndim = PyArray_NDIM(data); + const npy_intp *odim = PyArray_DIMS(data); + int npts = odim[ndim - 1]; + npy_intp *tdim=(npy_intp *)malloc(ndim*sizeof(npy_intp)); + if (!tdim) + { Py_XDECREF(data); return NULL; } + for (int d=0; d<ndim-1; ++d) + tdim[d] = odim[d]; + tdim[ndim-1] = npts/2 + 1; + PyArrayObject *ret = (PyArrayObject *)PyArray_Empty(ndim, + tdim, PyArray_DescrFromType(NPY_CDOUBLE), 0); + free(tdim); + if (!ret) fail=1; + if (!fail) { + int rstep = PyArray_DIM(ret, PyArray_NDIM(ret) - 1)*2; + + int nrepeats = PyArray_SIZE(data)/npts; + double *rptr = (double *)PyArray_DATA(ret), + *dptr = (double *)PyArray_DATA(data); + + Py_BEGIN_ALLOW_THREADS; + NPY_SIGINT_ON; + plan = make_rfft_plan(npts); + if (!plan) fail=1; + if (!fail) + for (int i = 0; i < nrepeats; i++) { + rptr[rstep-1] = 0.0; + memcpy((char *)(rptr+1), dptr, npts*sizeof(double)); + if (rfft_forward(plan, rptr+1, fct)!=0) {fail=1; break;} + rptr[0] = rptr[1]; + rptr[1] = 0.0; + rptr += rstep; + dptr += npts; + } + if (plan) destroy_rfft_plan(plan); + NPY_SIGINT_OFF; + Py_END_ALLOW_THREADS; + } + if (fail) { + Py_XDECREF(data); + Py_XDECREF(ret); + return PyErr_NoMemory(); + } + Py_DECREF(data); + return (PyObject *)ret; +} +static PyObject * +execute_real_backward(PyObject *a1, double fct) +{ + rfft_plan plan=NULL; + PyArrayObject *data = (PyArrayObject *)PyArray_ContiguousFromObject(a1, + NPY_CDOUBLE, 1, 0); + if (!data) return NULL; + int npts = PyArray_DIM(data, PyArray_NDIM(data) - 1); + PyArrayObject *ret = (PyArrayObject *)PyArray_Empty(PyArray_NDIM(data), + PyArray_DIMS(data), PyArray_DescrFromType(NPY_DOUBLE), 0); + int fail = 0; + if (!ret) fail=1; + if (!fail) { + int nrepeats = PyArray_SIZE(ret)/npts; + double *rptr = (double *)PyArray_DATA(ret), + *dptr = (double *)PyArray_DATA(data); + + Py_BEGIN_ALLOW_THREADS; + NPY_SIGINT_ON; + plan = make_rfft_plan(npts); + if (!plan) fail=1; + if (!fail) { + for (int i = 0; i < nrepeats; i++) { + memcpy((char *)(rptr + 1), (dptr + 2), (npts - 1)*sizeof(double)); + rptr[0] = dptr[0]; + if (rfft_backward(plan, rptr, fct)!=0) {fail=1; break;} + rptr += npts; + dptr += npts*2; + } + } + if (plan) destroy_rfft_plan(plan); + NPY_SIGINT_OFF; + Py_END_ALLOW_THREADS; + } + if (fail) { + Py_XDECREF(data); + Py_XDECREF(ret); + return PyErr_NoMemory(); + } + Py_DECREF(data); + return (PyObject *)ret; +} + +static PyObject * +execute_real(PyObject *a1, int is_forward, double fct) +{ + return is_forward ? execute_real_forward(a1, fct) + : execute_real_backward(a1, fct); +} + +static const char execute__doc__[] = ""; + +static PyObject * +execute(PyObject *NPY_UNUSED(self), PyObject *args) +{ + PyObject *a1; + int is_real, is_forward; + double fct; + + if(!PyArg_ParseTuple(args, "Oiid:execute", &a1, &is_real, &is_forward, &fct)) { + return NULL; + } + + return is_real ? execute_real(a1, is_forward, fct) + : execute_complex(a1, is_forward, fct); +} + +/* List of methods defined in the module */ + +static struct PyMethodDef methods[] = { + {"execute", execute, 1, execute__doc__}, + {NULL, NULL, 0, NULL} /* sentinel */ +}; + +#if PY_MAJOR_VERSION >= 3 +static struct PyModuleDef moduledef = { + PyModuleDef_HEAD_INIT, + "pocketfft_internal", + NULL, + -1, + methods, + NULL, + NULL, + NULL, + NULL +}; +#endif + +/* Initialization function for the module */ +#if PY_MAJOR_VERSION >= 3 +#define RETVAL(x) x +PyMODINIT_FUNC PyInit_pocketfft_internal(void) +#else +#define RETVAL(x) +PyMODINIT_FUNC +initpocketfft_internal(void) +#endif +{ + PyObject *m; +#if PY_MAJOR_VERSION >= 3 + m = PyModule_Create(&moduledef); +#else + static const char module_documentation[] = ""; + + m = Py_InitModule4("pocketfft_internal", methods, + module_documentation, + (PyObject*)NULL,PYTHON_API_VERSION); +#endif + if (m == NULL) { + return RETVAL(NULL); + } + + /* Import the array object */ + import_array(); + + /* XXXX Add constants here */ + + return RETVAL(m); +} diff --git a/numpy/fft/fftpack.py b/numpy/fft/pocketfft.py index d0df6fb48..794d13937 100644 --- a/numpy/fft/fftpack.py +++ b/numpy/fft/pocketfft.py @@ -26,9 +26,6 @@ n = n-dimensional transform (Note: 2D routines are just nD routines with different default behavior.) -The underlying code for these functions is an f2c-translated and modified -version of the FFTPACK routines. - """ from __future__ import division, absolute_import, print_function @@ -37,26 +34,18 @@ __all__ = ['fft', 'ifft', 'rfft', 'irfft', 'hfft', 'ihfft', 'rfftn', import functools -from numpy.core import (array, asarray, zeros, swapaxes, shape, conjugate, - take, sqrt) +from numpy.core import asarray, zeros, swapaxes, conjugate, take, sqrt +from . import pocketfft_internal as pfi from numpy.core.multiarray import normalize_axis_index from numpy.core import overrides -from . import fftpack_lite as fftpack -from .helper import _FFTCache - -_fft_cache = _FFTCache(max_size_in_mb=100, max_item_count=32) -_real_fft_cache = _FFTCache(max_size_in_mb=100, max_item_count=32) array_function_dispatch = functools.partial( overrides.array_function_dispatch, module='numpy.fft') -def _raw_fft(a, n=None, axis=-1, init_function=fftpack.cffti, - work_function=fftpack.cfftf, fft_cache=_fft_cache): - a = asarray(a) +def _raw_fft(a, n, axis, is_real, is_forward, fct): axis = normalize_axis_index(axis, a.ndim) - if n is None: n = a.shape[axis] @@ -64,15 +53,6 @@ def _raw_fft(a, n=None, axis=-1, init_function=fftpack.cffti, raise ValueError("Invalid number of FFT data points (%d) specified." % n) - # We have to ensure that only a single thread can access a wsave array - # at any given time. Thus we remove it from the cache and insert it - # again after it has been used. Multiple threads might create multiple - # copies of the wsave array. This is intentional and a limitation of - # the current C code. - wsave = fft_cache.pop_twiddle_factors(n) - if wsave is None: - wsave = init_function(n) - if a.shape[axis] != n: s = list(a.shape) if s[axis] > n: @@ -87,25 +67,22 @@ def _raw_fft(a, n=None, axis=-1, init_function=fftpack.cffti, z[tuple(index)] = a a = z - if axis != a.ndim - 1: + if axis == a.ndim-1: + r = pfi.execute(a, is_real, is_forward, fct) + else: a = swapaxes(a, axis, -1) - r = work_function(a, wsave) - if axis != a.ndim - 1: + r = pfi.execute(a, is_real, is_forward, fct) r = swapaxes(r, axis, -1) - - # As soon as we put wsave back into the cache, another thread could pick it - # up and start using it, so we must not do this until after we're - # completely done using it ourselves. - fft_cache.put_twiddle_factors(n, wsave) - return r def _unitary(norm): - if norm not in (None, "ortho"): - raise ValueError("Invalid norm value %s, should be None or \"ortho\"." - % norm) - return norm is not None + if norm is None: + return False + if norm=="ortho": + return True + raise ValueError("Invalid norm value %s, should be None or \"ortho\"." + % norm) def _fft_dispatcher(a, n=None, axis=None, norm=None): @@ -198,12 +175,13 @@ def fft(a, n=None, axis=-1, norm=None): """ - a = asarray(a).astype(complex, copy=False) + a = asarray(a) if n is None: n = a.shape[axis] - output = _raw_fft(a, n, axis, fftpack.cffti, fftpack.cfftf, _fft_cache) - if _unitary(norm): - output *= 1 / sqrt(n) + fct = 1 + if norm is not None and _unitary(norm): + fct = 1 / sqrt(n) + output = _raw_fft(a, n, axis, False, True, fct) return output @@ -294,13 +272,14 @@ def ifft(a, n=None, axis=-1, norm=None): >>> plt.show() """ - # The copy may be required for multithreading. - a = array(a, copy=True, dtype=complex) + a = asarray(a) if n is None: n = a.shape[axis] - unitary = _unitary(norm) - output = _raw_fft(a, n, axis, fftpack.cffti, fftpack.cfftb, _fft_cache) - return output * (1 / (sqrt(n) if unitary else n)) + fct = 1/n + if norm is not None and _unitary(norm): + fct = 1/sqrt(n) + output = _raw_fft(a, n, axis, False, False, fct) + return output @@ -383,14 +362,13 @@ def rfft(a, n=None, axis=-1, norm=None): exploited to compute only the non-negative frequency terms. """ - # The copy may be required for multithreading. - a = array(a, copy=True, dtype=float) - output = _raw_fft(a, n, axis, fftpack.rffti, fftpack.rfftf, - _real_fft_cache) - if _unitary(norm): + a = asarray(a) + fct = 1 + if norm is not None and _unitary(norm): if n is None: n = a.shape[axis] - output *= 1 / sqrt(n) + fct = 1/sqrt(n) + output = _raw_fft(a, n, axis, True, True, fct) return output @@ -475,14 +453,14 @@ def irfft(a, n=None, axis=-1, norm=None): specified, and the output array is purely real. """ - # The copy may be required for multithreading. - a = array(a, copy=True, dtype=complex) + a = asarray(a) if n is None: n = (a.shape[axis] - 1) * 2 - unitary = _unitary(norm) - output = _raw_fft(a, n, axis, fftpack.rffti, fftpack.rfftb, - _real_fft_cache) - return output * (1 / (sqrt(n) if unitary else n)) + fct = 1/n + if norm is not None and _unitary(norm): + fct = 1/sqrt(n) + output = _raw_fft(a, n, axis, True, False, fct) + return output @array_function_dispatch(_fft_dispatcher) @@ -560,8 +538,7 @@ def hfft(a, n=None, axis=-1, norm=None): [ 2., -2.]]) """ - # The copy may be required for multithreading. - a = array(a, copy=True, dtype=complex) + a = asarray(a) if n is None: n = (a.shape[axis] - 1) * 2 unitary = _unitary(norm) @@ -621,8 +598,7 @@ def ihfft(a, n=None, axis=-1, norm=None): array([ 1.-0.j, 2.-0.j, 3.-0.j, 4.-0.j]) # may vary """ - # The copy may be required for multithreading. - a = array(a, copy=True, dtype=float) + a = asarray(a) if n is None: n = a.shape[axis] unitary = _unitary(norm) @@ -1122,8 +1098,7 @@ def rfftn(a, s=None, axes=None, norm=None): [0.+0.j, 0.+0.j]]]) """ - # The copy may be required for multithreading. - a = array(a, copy=True, dtype=float) + a = asarray(a) s, axes = _cook_nd_args(a, s, axes) a = rfft(a, s[-1], axes[-1], norm) for ii in range(len(axes)-1): @@ -1255,8 +1230,7 @@ def irfftn(a, s=None, axes=None, norm=None): [1., 1.]]]) """ - # The copy may be required for multithreading. - a = array(a, copy=True, dtype=complex) + a = asarray(a) s, axes = _cook_nd_args(a, s, axes, invreal=1) for ii in range(len(axes)-1): a = ifft(a, s[ii], axes[ii], norm) diff --git a/numpy/fft/setup.py b/numpy/fft/setup.py index cd99a82d7..6c3548b65 100644 --- a/numpy/fft/setup.py +++ b/numpy/fft/setup.py @@ -7,9 +7,9 @@ def configuration(parent_package='',top_path=None): config.add_data_dir('tests') - # Configure fftpack_lite - config.add_extension('fftpack_lite', - sources=['fftpack_litemodule.c', 'fftpack.c'] + # Configure pocketfft_internal + config.add_extension('pocketfft_internal', + sources=['pocketfft.c'] ) return config diff --git a/numpy/fft/tests/test_helper.py b/numpy/fft/tests/test_helper.py index 8d315fa02..6613c8002 100644 --- a/numpy/fft/tests/test_helper.py +++ b/numpy/fft/tests/test_helper.py @@ -7,7 +7,6 @@ from __future__ import division, absolute_import, print_function import numpy as np from numpy.testing import assert_array_almost_equal, assert_equal from numpy import fft, pi -from numpy.fft.helper import _FFTCache class TestFFTShift(object): @@ -168,81 +167,3 @@ class TestIRFFTN(object): # Should not raise error fft.irfftn(a, axes=axes) - - -class TestFFTCache(object): - - def test_basic_behaviour(self): - c = _FFTCache(max_size_in_mb=1, max_item_count=4) - - # Put - c.put_twiddle_factors(1, np.ones(2, dtype=np.float32)) - c.put_twiddle_factors(2, np.zeros(2, dtype=np.float32)) - - # Get - assert_array_almost_equal(c.pop_twiddle_factors(1), - np.ones(2, dtype=np.float32)) - assert_array_almost_equal(c.pop_twiddle_factors(2), - np.zeros(2, dtype=np.float32)) - - # Nothing should be left. - assert_equal(len(c._dict), 0) - - # Now put everything in twice so it can be retrieved once and each will - # still have one item left. - for _ in range(2): - c.put_twiddle_factors(1, np.ones(2, dtype=np.float32)) - c.put_twiddle_factors(2, np.zeros(2, dtype=np.float32)) - assert_array_almost_equal(c.pop_twiddle_factors(1), - np.ones(2, dtype=np.float32)) - assert_array_almost_equal(c.pop_twiddle_factors(2), - np.zeros(2, dtype=np.float32)) - assert_equal(len(c._dict), 2) - - def test_automatic_pruning(self): - # That's around 2600 single precision samples. - c = _FFTCache(max_size_in_mb=0.01, max_item_count=4) - - c.put_twiddle_factors(1, np.ones(200, dtype=np.float32)) - c.put_twiddle_factors(2, np.ones(200, dtype=np.float32)) - assert_equal(list(c._dict.keys()), [1, 2]) - - # This is larger than the limit but should still be kept. - c.put_twiddle_factors(3, np.ones(3000, dtype=np.float32)) - assert_equal(list(c._dict.keys()), [1, 2, 3]) - # Add one more. - c.put_twiddle_factors(4, np.ones(3000, dtype=np.float32)) - # The other three should no longer exist. - assert_equal(list(c._dict.keys()), [4]) - - # Now test the max item count pruning. - c = _FFTCache(max_size_in_mb=0.01, max_item_count=2) - c.put_twiddle_factors(2, np.empty(2)) - c.put_twiddle_factors(1, np.empty(2)) - # Can still be accessed. - assert_equal(list(c._dict.keys()), [2, 1]) - - c.put_twiddle_factors(3, np.empty(2)) - # 1 and 3 can still be accessed - c[2] has been touched least recently - # and is thus evicted. - assert_equal(list(c._dict.keys()), [1, 3]) - - # One last test. We will add a single large item that is slightly - # bigger then the cache size. Some small items can still be added. - c = _FFTCache(max_size_in_mb=0.01, max_item_count=5) - c.put_twiddle_factors(1, np.ones(3000, dtype=np.float32)) - c.put_twiddle_factors(2, np.ones(2, dtype=np.float32)) - c.put_twiddle_factors(3, np.ones(2, dtype=np.float32)) - c.put_twiddle_factors(4, np.ones(2, dtype=np.float32)) - assert_equal(list(c._dict.keys()), [1, 2, 3, 4]) - - # One more big item. This time it is 6 smaller ones but they are - # counted as one big item. - for _ in range(6): - c.put_twiddle_factors(5, np.ones(500, dtype=np.float32)) - # '1' no longer in the cache. Rest still in the cache. - assert_equal(list(c._dict.keys()), [2, 3, 4, 5]) - - # Another big item - should now be the only item in the cache. - c.put_twiddle_factors(6, np.ones(4000, dtype=np.float32)) - assert_equal(list(c._dict.keys()), [6]) diff --git a/numpy/fft/tests/test_fftpack.py b/numpy/fft/tests/test_pocketfft.py index 8d6cd8407..0552f6afd 100644 --- a/numpy/fft/tests/test_fftpack.py +++ b/numpy/fft/tests/test_pocketfft.py @@ -28,6 +28,16 @@ class TestFFTShift(object): class TestFFT1D(object): + def test_identity(self): + maxlen = 512 + x = random(maxlen) + 1j*random(maxlen) + xr = random(maxlen) + for i in range(1,maxlen): + assert_array_almost_equal(np.fft.ifft(np.fft.fft(x[0:i])), x[0:i], + decimal=12) + assert_array_almost_equal(np.fft.irfft(np.fft.rfft(xr[0:i]),i), + xr[0:i], decimal=12) + def test_fft(self): x = random(30) + 1j*random(30) assert_array_almost_equal(fft1(x), np.fft.fft(x)) |