// Illustrate circular convolution thru fft.
#include <complex.h>
#include <math.h>
typedef double _Complex c;
// In this code ~ (which is not -) is the complex conjugate operator.
// http://gcc.gnu.org/onlinedocs/gcc-4.8.2/gcc/Complex.html#Complex

void fft(const int n, c a[n]){
  {int j=0, i; for(i=0; i<n-1; ++i) {
    if(i<j) {c t=a[j]; a[j]=a[i]; a[i]=t;}
    {int m=n/2; while(m<=j) {j -= m; m /= 2;} j += m;}}} 
  {int j=1, m; c b=-1; while(j<n) {c w = 1;
    for(m=0; m<j; ++m) {int i;
      for(i=m; i<n; i+=j+j)
        {c t = w*a[i+j]; a[i+j]=a[i]-t; a[i] += t;}
      w *= b;}
    j *= 2; b = csqrt(b);}}}

#include <stdio.h>
void pc(c a){printf("%16.12f + %16.12fi", creal(a),  cimag(a));}

#define D (1<<5)
int main(){c a[D], b[D], c[D];
  {int j=D; while(j--) {a[j] = (j==5 || j==29); b[j] = (j==2 || j==11);}}
  fft(D, a); fft(D, b); {int j=D; while(j--) c[j] = ~(a[j]*b[j])/D;}
  fft(D, c);
  {int j; for(j=0; j<D; ++j) if (cabs(c[j]) > 1.e-8) {pc(c[j]); printf(" %d\n", j);}}
  return 0;}
// c[k] = sum[i+j=k mod D] a[i]*b[j] ; convolution

/*
  1.000000000000 +  -0.000000000000i 7
  1.000000000000 +  -0.000000000000i 8
  1.000000000000 +  -0.000000000000i 16
  1.000000000000 +  -0.000000000000i 31
*/

// Correlation:
// Substituting (~a[j])*b[j] for ~(a[j]*b[j]) is for c[k] = sum[i-j=k mod D] a[i]*b[j]
/*
  1.000000000000 +   0.000000000000i 3
  1.000000000000 +   0.000000000000i 18
  1.000000000000 +   0.000000000000i 26
  1.000000000000 +   0.000000000000i 27
*/