#include <stdio.h>
#include <assert.h>

#define MAXBIT 30

void sobseq(float *x, float *y) {

  // Returns two quasi-random numbers for a 2-dimensional Sobel
  // sequence. Adapted from Numerical Recipies in C.

  int j, k, l;
  unsigned long i, im,ipp;

  // The following variables are "static" since we want their
  // values to remain stored after the function returns. These
  // values represent the state of the quasi-random number generator.

  static float fac;
  static int init=0;
  static unsigned long in, ix[3], *iu[2*MAXBIT+1];
  static unsigned long mdeg[3]={0,1,2};
  static unsigned long ip[3]={0,0,1};
  static unsigned long iv[2*MAXBIT+1]=
     {0,1,1,1,1,1,1,3,1,3,3,1,1,5,7,7,3,3,5,15,11,5,15,13,9};

  // Initialise the generator the first time the function is called.

  if (!init) {
    init = 1;
    for (j = 1, k = 0; j <= MAXBIT; j++, k += 2) iu[j] = &iv[k];
    for (k = 1; k <= 2; k++) {
      for (j = 1; j <= mdeg[k]; j++) iu[j][k] <<= (MAXBIT-j);
      for (j = mdeg[k]+1; j <= MAXBIT; j++) {
	ipp = ip[k];
	i = iu[j-mdeg[k]][k];
	i ^= (i >> mdeg[k]);
	for (l = mdeg[k]-1; l >= 1; l--) {
	  if (ipp & 1) i ^= iu[j-l][k];
	  ipp >>= 1;
	}
	iu[j][k] = i;
      }
    }
    fac = 1.0/(1L << MAXBIT);
    in = 0;
  }

  // Now calculate the next pair of numbers in the 2-dimensional
  // Sobel sequence.

  im = in;
  for (j = 1; j <= MAXBIT; j++) {
    if (!(im & 1)) break;
    im >>= 1;
  }
  assert(j <= MAXBIT);
  im = (j-1)*2;
  ix[1] ^= iv[im+1];
  ix[2] ^= iv[im+2];
  *x = ix[1] * fac;
  *y = ix[2] * fac;
  in++;
}


int main(void) {

  // Test the Sobel generator.

  int i;
  float x, y;

  for (i = 0; i < 1000; i++) {
    sobseq(&x, &y);
    printf("%f %f\n", x, y);
  }
  return 0;
}