I have 2 pieces of code called GENERO1.c and gauss2.c which will generate two random numbers in Polar marsaglia method. However, something is not working properly because when I run another code where I call the function gauss2() , its stopping in that function like being in an infinite loop. Could anyone tell me why is not working?
Bellow is the code?
Please I would welcome any help, because I really cant find something.
I compile the program with:
gcc -o check -lm gauss2.c GENER01.c check.c
and run by: ./check
GENER01.c
Code:/* NOTE works for 32 bit !!!!!!!!!! use "-m32" */ #include <stdio.h> #include<math.h> #define znew (z=36969*(z&65535)+(z>>16)) #define wnew (w=18000*(w&65535)+(w>>16)) #define MWC ((znew<<16)+wnew ) #define SHR3 (jsr^=(jsr<<17), jsr^=(jsr>>13), jsr^=(jsr<<5)) #define CONG (jcong=69069*jcong+1234567) #define FIB ((b=a+b),(a=b-a)) #define KISS ((MWC^CONG)+SHR3) #define LFIB4 (c++,t[c]=t[c]+t[UC(c+58)]+t[UC(c+119)]+t[UC(c+178)]) #define SWB (c++,bro=(x<y),t[c]=(x=t[UC(c+34)])-(y=t[UC(c+19)]+bro)) #define UNI (KISS*2.328306e-10) #define VNI ((long) KISS)*4.656613e-10 #define UC (unsigned char) /*a cast operation*/ typedef unsigned long UL; /* Global static variables: */ static UL z=362436069, w=521288629, jsr=123456789, jcong=380116160; static UL a=224466889, b=7584631, t[256]; /* Use random seeds to reset z,w,jsr,jcong,a,b, and the table t[256]*/ /*static UL x=0,y=0,bro; static unsigned char c=0;*/ /* Example procedure to set the table, using KISS: */ void settable(UL i1,UL i2,UL i3,UL i4,UL i5, UL i6) { int i; z=i1;w=i2,jsr=i3; jcong=i4; a=i5; b=i6; for(i=0;i<256;i=i+1) t[i]=KISS; } int Threepoint(){ return(KISS%3); } int Eightpoint(){ return(KISS%8); } int Fourpoint(){ return(KISS%4); } double GENER02(double lambda){ /* lambda > 0 */ /* function is a translation from P.E. Kloeden E. Platen and H. Schurz Numerical Solution of SDE Trough Computer experiments, Springer 1994. page 11. */ double xu; double Uniform(); xu=Uniform(); if(xu>0.0) xu=-log(xu)/lambda; else xu=0.0; return(xu); } double Uniform(){ return((double) UNI); } double GENER01(double sqrtDT){ /* *WT=-sqrtDT+2*(KISS%2)*sqrtDT; */ if(KISS%2) return(-sqrtDT); else return(sqrtDT); } double GENER03(double sqrt3DT){ long I; double Y; /* generate a random variable on set {1,2,3,4,5,6} with uniform prob. */ I = KISS%6+1; /* six point random variable 1, 2,...,6. Each element prob. 1/6 */ if(I>2) /* prob. 2/3 */ Y=0.0; else{ if(I==1) /* prob. 1/6 */ Y=-sqrt3DT; else /* prob. 1/6 */ Y=sqrt3DT; } /* printf("%g %d\n",Y,I); */ /* printf("%g\n",Y); */ return( Y ); }
gauss2.c
check.cCode:#include<math.h> /* Polar Marsaglia method for gaussian dist. pseudo-random no.generator*/ void gauss2(double *z1,double *z2){ double y, v1, v2,r, fac; double Uniform(); double temp; r=100.0; while(r>1.0 || r<=0.0){ y=Uniform(); /* rand. no [0:1] */ v1 = 2.0*y - 1.0; /* rand. no [-1:1] */ y=Uniform(); v2 = 2.0*y - 1.0; /* rand. no [-1:1] */ r = v1 * v1 + v2 * v2; } fac = (sqrt( -2.0*log(r) / r)); *z1 = v1 * fac; *z2 = v2 * fac; return; }
Code:#include <stdlib.h> #include <stdio.h> #include <time.h> #include <math.h> typedef unsigned long UL; #define pi 3.14159 double Uniform(); void gauss2(double *z1,double *z2); main() { int i; int j=100; int new; double z1, z2; UL i1, i2,i3,i4,i5,i6; /* to initialize Uniform (AND gauss2) */ i1=time(NULL); i2= i1+23;i3= i1+243;i4= i1+26;i5= i2+28;i6= i1+27; settable(i1,i2,i3,i4,i5,i6); for(i=0;i<j;i++) { printf("OK\n"); gauss2(&z1,&z2); printf("%lf",z1); } }
