#include <stdio.h>
#include <math.h>

#define NMAX 8192
#define NMAXSQRT 64

void cdft(int, int, double *, int *, double *);

int main()
{
    int n, ip[NMAXSQRT + 2];
    double a[NMAX + 1], w[NMAX * 5 / 4];

    n = 1024;		  /* # point of DFT */
    if (2 * n > NMAX) {
        printf("2 * %d > %d (NMAX)\n", n, NMAX);
        return(1);
    }
    printf("n %d\n", n);

    ip[0] = 0;

    /* cdft(), complex fast dft */
    int i;
    double* p;
    p = a;
    for (i = 0; i < n; i += 1) {
        *p++ = cos(3 * 2*M_PI*i/n) / 3
             + cos(5 * 2*M_PI*i/n) / 5
             + cos(7 * 2*M_PI*i/n) / 7;
        *p++ = 0;
    }

    cdft(2*n, 1, a, ip, w);

    p = a;
    for (i = 0; i < n / 2; i += 1, p += 2) {    /* n/2 for amplitude spectrum */
        printf("%4d: %16f %16f %16f\n",
            i, *(p+0), *(p+1),
            2 * sqrt(*(p+0) * *(p+0) + *(p+1) * *(p+1)) / n);
    }
    printf("\n");

    return 0;
}