High precision Discrete Fourier Transform in c

Question

I'm trying to do a high precision discrete fourier transform on a signal. To examine the precision, I use a gaussian function as the signal, because the fourier transform is also a gaussian function.

The code is as follows

#include <stdio.h>
#include <tgmath.h>

typedef long double complex cplx;
long double PI;
void _fft(cplx buf[], cplx out[], int n, int step)
{
if (step < n) {
    _fft(out, buf, n, step * 2);
    _fft(out + step, buf + step, n, step * 2);

    for (int i = 0; i < n; i += 2 * step) {
        cplx t = cexpl(- I * PI * i / n) * out[i + step];
        buf[i / 2]     = out[i] + t;
        buf[(i + n)/2] = out[i] - t;
    }
}
}

void fft(cplx buf[], int n)
{
cplx out[n];
for (int i = 0; i < n; i++) out[i] = buf[i];

_fft(buf, out, n, 1);
}


int main()
{
const int nPoints = pow(2, 12);
PI = atan2(1, 1) * 4;
cplx dt = 1e-3;
cplx dOmega = 1 / (dt * nPoints);
long double T[nPoints];
long double DOmega[nPoints];

cplx At[nPoints];
cplx tau = 28.4e-3;
for (int i = 0; i < nPoints; ++i)
{
    T[i] = dt*(i-nPoints/2);
    DOmega[i] = dOmega * (i - nPoints / 2);
    At[i] = cexpl(-T[i]*T[i]/2/(tau*tau));
}

fft(At, nPoints);
FILE* fw;
fw = fopen("fft_01.txt", "w+");

for (int i = 0; i < nPoints; ++i)
{
    fprintf(fw, "%.15Le, %.15Le\n", DOmega[i], fabs(At[i]) );
}

return 0;
}

When I change the defined type from float complex to double complex, the result did show a improvement in precision. However, When I change the data type from double complex to long double complex, there's no improvement while the numbers are far away from the minimal limit.

I think that the two wings of the peak should go down like the results by Mathematica, but here it seems stuck around 10^(-35) and can not further get down.

When I further apply a iDFT to the result of DFT, it's more apparent that the long double precision has been lost. The blue one is the original signal, with high precision. The green one is the result after iDFT(DFT(x)), which should be coincide with the blue one, but the minimal numbers are much larger.

Anyone can tell me where is the problem, thanks in advance.

PS: I am on Mac OS 10.13.1, using gcc as the compiler. I have checked the numerical limits on my platform

storage size for float: 4
minimum float positive value: 1.175494e-38
Maximum float positive value: 3.402823e+38
precision value: 6

storage size for double: 8
minimum double positive value: 2.225074e-308
Maximum double positive value: 1.797693e+308
precision value: 15

storage size for long double: 16
minimum long double positive value: 3.362103e-4932
Maximum long double positive value: 1.189731e+4932
precision value: 18

Answer 1

In addition to the various comments that point out that most built-in functions such as atan or exp operate on double precision variables, there is also the issue that long double on typical x86-type hardware is not actually much more accurate (if at all) than regular double. In particular, long double is not an arbitrary precision data type such as the ones that Mathematica uses -- Mathematica lets you choose the accuracy of variables however you want it, and it will compute things to this accuracy, but this is not the case for C's long double data type that has a fixed accuracy just like double.