0

i am working with monte carlo simulations and i need to understand this function in C language. I hope someone knows how to explain the modus operandi. Thx. // function that gen random numbers

        #include <stdlib.h>
        #define MBIG 1000000000
        #define MSEED 161803398
        #define MZ 0
        #define FAC (1.0/MBIG)

        double ran3(long *idum)
        {
          static int inext,inextp;
          static long ma[56];
          static int iff=0;
          long mj,mk;
          int i,ii,k;

          if (*idum < 0 || iff == 0){
              iff=1;
              mj = labs(MSEED-labs(*idum));
              mj %= MBIG;
              ma[55]=mj;
              mk=1;
              for (i=1;i<=54;i++){
                  ii=(21*i) % 55;
                  ma[ii]=mk;
                  mk=mj-mk;
                  if (mk < MZ) mk += MBIG;
                  mj =ma[ii];
              }
              for (k=1;k<=4;k++)
                  for (i=1;i<=55;i++) {
                      ma[i] -= ma[1+(i+30) % 55];
                      if (ma[i] < MZ) ma[i] += MBIG;
                  }
              inext=0;
              inext=31;
              *idum=1;
          }
          if (++inext == 56) inext=1;
          if (++inextp == 56) inextp=1;
          mj =ma[inext]-ma[inextp];
          if (mj < MZ) mj += MBIG;
          ma[inext] =mj;
          return mj*FAC;
        }
John Paul
  • 15
  • 1
  • 5
    I'd definitely expect something random to come out of this jumble. – Quentin Feb 08 '15 at 01:50
  • 2
    That's actually Knuth's Subtractive RNG. See e.g. http://stackoverflow.com/a/23691787/214671 or http://rosettacode.org/wiki/Subtractive_generator – Matteo Italia Feb 08 '15 at 01:52

1 Answers1

1

That function is from Numerical Recipes in C, Chapter 7, Random Numbers. The best explanation you will find is in there.

  • Alas, there's no explanation of the algorithm in "Numerical Recipes". There is in one in "The Art of Computer Programming" by Don Knuth. – Wintermute Feb 08 '15 at 01:58