#include #include #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; }