mardi 20 juillet 2021

Gaussian Random Number Generator using ran2.c and boxmuller algorithm

I am trying to generate normally distributed random numbers using the 'ran2.c' and Box-Muller Method. I have written the code and saved as 'gau_rit.c' . It is showing this warning message,

In file included from gau_rit.c:5:
ran2.c:16:19: note: expected ‘long int *’ but argument is of type ‘long int’

Also for the output,

Segmentation fault (core dumped)

I am new in C language, so I need help in this regard. I understand that something is happening with the variable declaration. But failed to figure out.

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#define PI         3.14159265   // The value of pi
#include "ran2.c"
double norm(double mean, double sd);  // Returns a normal random variable
//double ran2(int seed);   
int main()
{
int i=45;//counter 
float x,s;
double mean = 0.0, sd = 1.0;
static long seed1=123456789;
static long seed2=987654321;
double norm(double mean, double sd){
    double  u, r, theta;
    u = 0.0;
    while(u==0.0)
        u=ran2(seed1);
    r = sqrt(-2.0 * log(u));
    theta = 0.0;
    while(theta==0.0)
        theta = 2.0 * PI * ran2(seed2);
    x = r * cos(theta); //for arbitrary mean and sd, we need to change this x
    return(x);
}

for(i=1;i<=45;i++)
{
  s=norm(mean,sd);
  printf("%f\n",s);
 }
}

And the ran2.c is here

#define IM1 2147483563
#define IM2 2147483399
#define AM (1.0/IM1)
#define IMM1 (IM1-1)
#define IA1 40014
#define IA2 40692
#define IQ1 53668
#define IQ2 52774
#define IR1 12211
#define IR2 3791
#define NTAB 32
#define NDIV (1+IMM1/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)

double ran2(long *idum)
{
  int j;
  long k;
  static long idum2=123456789;
  static long iy=0;
  static long iv[NTAB];
  double temp;

  if (*idum <= 0) {
    if (-(*idum) < 1) *idum=1;
    else *idum = -(*idum);
    idum2=(*idum);
    for (j=NTAB+7;j>=0;j--) {
        k=(*idum)/IQ1;
        *idum=IA1*(*idum-k*IQ1)-k*IR1;
        if (*idum < 0) *idum += IM1;
        if (j < NTAB) iv[j] = *idum;
    }
    iy=iv[0];
}
k=(*idum)/IQ1;
*idum=IA1*(*idum-k*IQ1)-k*IR1;
if (*idum < 0) *idum += IM1;
k=idum2/IQ2;
idum2=IA2*(idum2-k*IQ2)-k*IR2;
if (idum2 < 0) idum2 += IM2;
j=iy/NDIV;
iy=iv[j]-idum2;
iv[j] = *idum;
if (iy < 1) iy += IMM1;
if ((temp=AM*iy) > RNMX) return RNMX;
else return temp;
}
#undef IM1
#undef IM2
#undef AM
#undef IMM1
#undef IA1
#undef IA2
#undef IQ1
#undef IQ2
#undef IR1
#undef IR2
#undef NTAB
#undef NDIV
#undef EPS
#undef RNMX
/* (C) Copr. 1986-92 Numerical Recipes Software #.3. 



Aucun commentaire:

Enregistrer un commentaire