// http://en.wikipedia.org/wiki/Fast_Hadamard_transform
#include <stdio.h>
typedef double r;
void hm(int m, r * x){int b; for(b=m; b>1; b>>=1)
{int c=b>>1; int k;
 if(0) printf("bb %d %d\n", b, c);
 for(k=0; k<m; k+= b) {int j; for(j=0; j<c; ++j) {
  r t = x[j+k+c];
  x[j+k+c] = x[j+k] - t;
  if(0) printf("sw j=%d k=%d c=%d\n", j, k, c);
  x[j+k] += t;}}}}

#include <math.h>
r f(r x){return x+3;}
#define s (1<<22)
// The shift amount above can be an integer from 0 thru at least 26.
r w[s];
void p(r a) {int j; if(0)for(j=0; j<s; ++j) {if(!(j%8)) printf("\n%4d ", j); printf("%e ", a*w[j]);}}
int main(){
{int j=s; while(j--) w[j]=f(j);}
p(1); hm(s, w); p(1./sqrt(s)); hm(s, w); p(1./s);
{int j=s; while(j--) if(w[j]!=f(j)*s) printf("gox %d %13.12f\n", j, w[j]);}
return 0;}