/* Gaussian random number sequence generator by Walsh transform (WT) */
/* uniform disribution ---> Guassian distribution                    */
/*                 g[] ---> g[] (inplace processing)                 */
/*                                                                   */
/* (c) 2025 cepstrum.co.jp                                           */

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

#define WALSH_N   1048576          // 1048576=2^20

float g[WALSH_N];                  // Walsh transform in/out
float g2[WALSH_N];                 // work space for Walsh transform
float tmp[WALSH_N];                // save input sequence

#define SEED_VAL  0x5a5aa5a5       // xorshift seed value
#define OUT_FNAME "wt_out.txt"     // output data filne name

//================================================================
// integer log2()

unsigned ilog2(unsigned x) {
  unsigned i, j;

  i=x;
  j=0;
  do {
    j=j+1;
    i=i/2;
  } while (i!=1);
  return j;
}

//================================================================
// fast Walsh transform routine
// converted from FORTRAN code in the old book published more than 40 years ago
// floor() is meaningless in C program

void walsh(unsigned length, int flag) {
  unsigned i1, i2, i3, j, k, n1, n2, n3;
  int      pn;

  for (i1=0; i1<ilog2(length); i1=i1+1) {
    i2=pow(2, i1+1);
    for (j=0; j<i2; j=j+1) {
      i3=length/powf(2, i1+1);
      for (k=0; k<i3; k=k+1) {
        n1=k*powf(2, i1+1)+j;
        n2=2*k*powf(2, i1)+floor(j/2);
        n3=(2*k+1)*powf(2, i1)+floor(j/2);
        pn=pow(-1, floor((j+1)/2));
        g2[n1]=g[n2]+pn*g[n3];
        //printf("%6u %6u %6u %6i\n", n1, n2, n3, pn);   // for test
      }
    }
    for (j=0; j<length; j=j+1) g[j]=g2[j];
  }
  if (0<flag) {
    for (j=0; j<length; j=j+1) g[j]=g[j]/(float)length;
  }
}

//================================================================
// xorshift random number generator (32bit)

uint32_t uint_xorshift32bit(uint32_t seed) {
  static uint32_t y=0x0f0fa5a5;

  if (seed!=0) {
    y=seed;
    return seed;
  }

  y^=(y<<5);
  y^=(y>>13);
  y^=(y<<6);
  return y;
}

//================================================================

int main() {
  FILE     *fp;
  unsigned i;
  float    offset, maxval;

  uint_xorshift32bit(SEED_VAL);     // set xorshift seed

  // generate uniform distribution signal
  for (i=0; i<WALSH_N; i=i+1) g[i]=uint_xorshift32bit(0);
  // remove DC offset
  offset=0.0;
  for (i=0; i<WALSH_N; i=i+1) offset=offset+g[i];
  for (i=0; i<WALSH_N; i=i+1) g[i]=g[i]-offset/(float)WALSH_N;
  // normalize [-1.0, 1.0]
  maxval=0.0;
  for (i=0; i<WALSH_N; i=i+1) {
    if (maxval<fabs(g[i])) maxval=fabs(g[i]);
  }
  for (i=0; i<WALSH_N; i=i+1) g[i]=g[i]/maxval;

  // save source data
  for (i=0; i<WALSH_N; i=i+1) tmp[i]=g[i];

  walsh(WALSH_N, -1);                // Walsh transform

  // output Gaussian random sequence
  fp=fopen(OUT_FNAME, "w");
  for (i=0; i<WALSH_N; i=i+1) fprintf(fp, "%e %e\n", tmp[i], g[i]);
  fclose(fp);

}