#include <complex.h>
#include <math.h>
typedef double _Complex c;

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));}

#include <stdlib.h>
double rnd(){static int f = 0;
double u2(){return 0x1.p-30*random()-1;}
  static double m;
  if(f) {f = 0; return m;}
  while (1){double x = u2(), y = u2(), s = x*x + y*y; 
    if(s<1) {double r = sqrt(- 2*log(s)/s);
    f = 1; m = x*r; return y*r;}}}

#include <assert.h>
#define M 41227
#define p 0
#define D 1000
#define CS (1<<12)
#define rp 25
#define s 0.16
// Signal at receiver:
float R[rp*(CS-D)+D];
int main(){
assert(D < CS);
float sk[D]; {int j=D; while(j--) sk[j] = (random()&16)?1:-1;}
c cc(int j) {c su = 0; int k = D; while(k--) su += sk[k]*R[k+j]; return su;}
{const int ro = M; assert (0<=ro && ro+D <= rp*(CS-D));
 {int j=rp*(CS-D)+D; while(j--) R[j] = 1?rnd():0;}
 {int j=D; while(j--) R[j+ro] += s*sk[j];}} // add signal

 {c sgc[CS]; {int j=D; while(j--) sgc[j] = sk[j];}
  {int j; for(j=D; j<CS; ++j) sgc[j] = 0;}
  fft(CS, sgc);
if(0) {// receive logic:
  c ws[CS];
  int k=rp; while(k--) {{int j=CS; while(j--) ws[j] = R[k*(CS-D)+j];}
  if(0) {int i=CS; while(i--) {pc(ws[i]); printf("x %d %d\n", i, k);}}
  fft(CS, ws);
  {int j=CS; while(j--) ws[j] = sgc[j]*~ws[j]/CS;} fft(CS, ws);
{int i=CS; while(i--)
  if(k==p && 0 || creal(ws[i])>120) {
     pc(ws[i]); printf(" k=%d i=%d t=%d\n", k, i, k*(CS-D)+i);}}}
if(1) {void px(int i){pc(cc(i)); printf(" sp %d\n", i);}
        px(M), px(1); px(0);}}
else{ // compute C
    {int j=CS; while(j--) sgc[j] *= ~sgc[j]/CS;} fft(CS, sgc);
    {double mx = -100; int wh;
      {int j=CS; while(j--) {double x = creal(sgc[j]);
      if(x>mx && j>0) {mx=x; wh=j;}
      if(fabs(x) > 0.01) {pc(sgc[j]); printf(" %d\n", j);}}}
      printf("Max C = %10.6f at %d\n", mx, wh);}}}
return 0;}