/*
 *  libbigint is a big integer arithmetic library
 *  Copyright (C) 2001 Ray Strode
 *
 *  This library is free software; you can redistribute it and/or
 *  modify it under the terms of the GNU Lesser General Public
 *  License as published by the Free Software Foundation; either
 *  version 2.1 of the License, or (at your option) any later version.
 *
 *  This library is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 *  Lesser General Public License for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public
 *  License along with this library; if not, write to the Free Software
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 */

/*
 * This is an implementation of the Discrete Fourier Transformation
 * (iterative FFT).  It is incomplete, and not currently used by libbigint.
 * The algorithm comes from the CLRS.
 */

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

#define lg(n) ((int) (log10(n) / .30102999566398))
#define delta .00000001

typedef struct {
    double real;
    double imaginary;
} complex;

complex complex_mul(complex A, complex B) {
    complex C;

    C.real = (A.real * B.real) - (A.imaginary * B.imaginary);
    C.imaginary = (A.real * B.imaginary) + (A.imaginary * B.real);

    return C;
}

complex complex_add(complex A, complex B) {
    complex C;

    C.real = A.real + B.real;
    C.imaginary = A.imaginary + B.imaginary;

    return C;
}

complex complex_sub(complex A, complex B) {
    complex C;

    C.real = A.real - B.real;
    C.imaginary = A.imaginary - B.imaginary;

    return C;
}

complex complex_make_nice(complex A) {

    complex B;

    if (abs(((long) A.real) - A.real) < delta) {
        B.real = 0;
    } else {
        B.real = A.real;
    }

    if (abs(((long) A.imaginary) - A.imaginary) < delta) {
        B.imaginary=0;
    } else {
        B.imaginary = A.imaginary;
    }

    return B;
}

int rev(int n, int nbits) {

    int i = 0, r = 0;

    for (i = 0; i < nbits; i++) {

        /*
         * Grabs a bit from n at the ith position and places it in the mirror
         * of it's position in r.
         */
        r += (((n & (1 << i)) != 0) << (nbits - i - 2));
    }

    return r;
}

void bit_reverse_copy(complex *a, int n, complex **A) {

    int k;
    int r;
    int nbits;

    if (a == NULL || n == 0 || A == NULL || a == *A) return;

    nbits = lg(n) + 1;

    *A = (complex *) malloc(n * sizeof(complex));

    for (k = 0; k < n; k++) {
        int r = rev(k, nbits);
        (*A)[r] = a[k];
    }

}

void dft(complex *a, int n, complex **output) {

    complex *A;
    int s, m, j, k, o;
    complex omega, omega_m, t, u;

    bit_reverse_copy(a, n, &A);

    for (s = 1; s <= lg(n); s++) {

        m = 1 << s;

        /*
         * omega_m = e^(2pi/m) (because e^ui = cos u + i(sin u))
         */
        omega_m.real = cos((2 * M_PI)/m);
        omega_m.imaginary = sin((2 * M_PI)/m);
        omega_m = complex_make_nice(omega_m);

        for (k = 0; k < n; k += m) {

            omega.real = 1;
            omega.imaginary = 0;

            for (j = 0; j < (m / 2); j++) {

                t = complex_mul(omega, A[k + j + (m/2)]);

                u = A[k + j];

                A[k + j] = complex_add(u, t);

                A[k + j + (m/2)] = complex_sub(u, t);

                /*
                 * FIXME: Recomputing cosine and sine each time is inefficient.
                 *        But the obvious optimization doesn't work.
                 */
#if 0
                omega = complex_mul(omega, omega_m);
#endif
                omega.real = cos(2 * (j + 1) * M_PI/(m));
                omega.imaginary = sin(2 * (j + 1) * M_PI/(m));
            }
        }
    }
    *output = A;
}

complex *inverse_dft(complex *a, int n) {

    /*
     * FIXME: implement the inverse
     */
}

main() {
    int i;
    complex A[8];
    complex *b, *c;

    A[0].real = 4;
    A[0].imaginary = 0;
    A[1].real = 0;
    A[1].imaginary = 0;
    A[2].real = 3;
    A[2].imaginary = 0;
    A[3].real = 6;
    A[3].imaginary = 0;
    A[4].real = 2;
    A[4].imaginary = 0;
    A[5].real = 9;
    A[5].imaginary = 0;
    A[6].real = 6;
    A[6].imaginary = 0;
    A[7].real = 5;
    A[7].imaginary = 0;

    dft(A, 8, &b);

    for (i = 0; i < 8; i++) printf("%g + %gi\n", b[i].real, b[i].imaginary);

    return 0;
}