#include <stdlib.h>
#include <math.h>
#include <stdio.h>
// Each invocation of rnd returns a new random normal deviate.
typedef double R;
R sq(R x){return x*x;}
R scale = 0x1.p-30;
R rnd(){static int f = 0;
  static R m;
  if(f) {f = 0; return m;}
  while (1){R x = scale*random()-1, y = scale*random()-1,
    s = x*x +y*y; 
    if(s<1) {R r = sqrt(- 2*log(s)/s);
    f = 1; m = x*r; return y*r;}}}

typedef double _Complex C;
typedef float _Complex c;
C cmp(R x, R y){C z; __real__ z = x; __imag__ z = y; return z;}
R m2(C p){return sq(__real__ p) + sq(__imag__ p);}
void pc(c a){printf("%14.7e + %14.7ei\n", __real__ a, __imag__ a);}
C it(C w){C r; __real__ r = - __imag__ w;
   __imag__ r = __real__ w; return r;}
C Sqrt(C x){
R A = __real__ x, B = __imag__ x;
if(A < 0) return it(Sqrt(-x));
{R r = sqrt((A + sqrt(A*A + B*B))/2);
return r + it(B/(2*r));}}

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 = Sqrt(b);}}}

#define P (1<<20)
R const tp = 2.*3.1415926535897932;
c d[P];
void inir(){for(int j=0; j<P; ++j) d[j] = cmp(rnd(), rnd());}
int main(){
if(0){inir();
fft(P, d);
{R s = 0; C t=0;
  for(int j=0; j<P; ++j) {s += m2(d[j]); t += d[j];}
  printf("s = %e rms=%e av=%e_%ei\n", s, sqrt(s)/P, __real__ t/P, __imag__ t/P);}}
if(0){R s=0, m=0; for(int j = 0; j<7000; ++j) {R d = rnd(); m += d; s += sq(d);}
   printf("m=%e s=%e\n", m, s);}
if(1){for(int j=0; j<P; ++j) d[j] = rnd() + 0.01*cos(326447*tp*j/P);
   fft(P, d);
   for(int j=0; j<P; ++j) {R n = m2(d[j]);
      if (n > 11200000) {printf("j=%7d n=%e ", j, n); pc(d[j]);}}}
return 0;}