/*  sfft.h  by K.Tsuru  */
/***************************************************
Provides fast Fourie transform(FFT).

Aug 24, 2000
After version 2.10 Ooura's FFT program is used and
"fftMaxArraySize" is enlarged.
System    fftType       table   USES_SIN_TABLE_FFT | USES_LONG_DBL_FFT | fftSqrt |  sin, cos
 32 bits   double            O          1               |        0                  |  sqrt    |    X
                                  X          0               |        0                 |  sqrt    |   sin, cos
 64 bits   double            O          1               |        0                 |  sqrt    |    X
                                  X          0               |        0                 |  sqrt    |   sin, cos
             long double      O          1               |        1                 |  sqrtl  |    X
                                  X          0               |        1                  |  sqrtl  |  sinl, cosl
---------------------------------------
long double ldpi =        3.1415926535 89793 23846; // give 21 digits
cout << ldpi << endl; // 3.1415926535 89793 11600 //16 digits agree only.  The same as double.
Then "long double" cannot be used FFT multiplication GCC at the present time.

****************************************************/
#ifndef S_FFT_H
#define S_FFT_H

#ifndef  DEF_CONST_H
#include "defconst.h" // for M_PI_4 
#endif
#ifndef SN_CONST_H
#include "snconst.h"
#endif  // SN_CONST_H
//type of real number using in FFT
//If you use long double or quadruple precision, change here.
//specify the function of sqrt(), modf(), etc.
//Speed depends on the performance of FPU.
#include <math.h>
#define fftM_PI_4 M_PI_4 // =pi/4 Has 20 digits enough for long double
#define USES_SIN_TABLE_FFT 0
#define USES_LONG_DBL_FFT 0

#if USES_LONG_DBL_FFT  // =uses long double============================
const bool Uses_long_double_FFT = true;
const bool Uses_double_FFT = false;
typedef long double fftType;
#define fftSqrt sqrtl
fftType fftModfl(const fftType x, fftType *ip); // modfl() my version // for FFTVerify()
#define fftFmod fmodl //
#define fftFabs fabsl //
#define fftLdexp ldexpl //
#define fftSin sinl
#define fftCos cosl
#define fftType_MANT_DIG LDBL_MANT_DIG // = 64
const int FFT_MAX_SIZE_BITS = 27; // bits of maximum size of FFT work area, fftMaxSize=2^27 = 134217728.
#else // end of USES_LONG_DBL_FFT
// below === uses double =====================
const bool Uses_long_double_FFT = false;
const bool Uses_double_FFT = true;
typedef double fftType;
#define fftSqrt sqrt
#define fftModf modf
#define fftFmod fmod //
#define fftFabs fabs
#define fftLdexp ldexp //
#define fftSin sin
#define fftCos cos
#define fftType_MANT_DIG DBL_MANT_DIG // = 53
const int FFT_MAX_SIZE_BITS = 23;   // 20 --> 22 --> 23 --> 24(NG); // bits of maximum size of FFT work area fftMaxSize=2^23 = 8388608.
#endif // end of === uses double =====================

const uint fftMaxArraySize = 1u << FFT_MAX_SIZE_BITS;

#if USES_SIN_TABLE_FFT
const bool UsesFFTSineTable = true;
long FFTSineTableSize();
long FFTSineTableUsedMemory();
int MakeSinTable(int maxsizebit);    // give FFT_MAX_SIZE_BITS
fftType fftSinR(long j, long w);        // = sin (pi/4)*(j/w)
fftType fftCosR(long j, long w);    // = cos (pi/4)*(j/w)
#else
const bool UsesFFTSineTable = false;
fftType fftSinR(long n, long d);
fftType fftCosR(long n, long d);
#endif // USES_SIN_TABLE_FFT

// return the information of work area on static memory.
uint FFTSize();     // figure of N
uint FFTsinSize();  // table size of sine.
uint FFTWorkSize(); // size of work area for sine/cosine table and data
                        // To get memory size please multiply "sizeof(fftType)"
                        // which will be order of 10 MB.
/*************************** To use Ooura's rdft(). ***************************/
struct FFTWorkArea{
    int size, ipsz;
    int* ip;    // for bit reverse table
    fftType* w; // for sine/cosine table
//constructor
    FFTWorkArea():size(0),ipsz(0),ip(NULL),w(NULL){ }
    void Allocate(int N);   // 217
    void MemFree();     // 217
    ~FFTWorkArea(){}
};
void rdft(int n, int isgn, fftType *a, int *ip, fftType *w);
/*******************************************************************************/
#endif      // S_FFT_H