FW: bad rand generator: distribution/prevalences/periodicity

[afanofosc]

There is no API function called srandEx.

Look at the source code and you can see what exactly these values do.

Read wikipedia and other internet articles about the park-miller/lehmer LCG to find out what seeds are appropriate.

I would use the C FAQ recommendation for getting subsets of the random number range but from what I have read the values I am using with the RNG are "modulus compatible".

Since Seed is a signed long any value greater than its maximum positive value (LONG_MAX) will wrap negative and any negative number will cause the seed to be set back to the previous seed.

John Hansen

[afanofosc]

edit: Doc deleted his last post for some reason...

You call srand to seed the RNG and you call rand to get a random number. I have already said that on this thread. There is no timer function so you can look at the source code I posted a link to to see what it does with seed values that are equal to 0 to less than zero. Just run all your same tests with srand(1) or (2) or whichever seed you used in your other tests and %16 or %100 on rand() and see what happens.

John Hansen

[HaWe]

ooops, sorry, wanted to write new....
but anyway, the reply says it all..
well, I have no doubt that the new functon works like documented - opposite to the old rand function.

[linusa]

I posted the quotation because I thought those who following this thread would enjoy von Neumann's joke.

Thanks, I really did.

[timpattinson]

I made a little program that displays targets on the screen in Random() places.
I ran it for 5mins and i looked like this:

nxt.gif (3.9 KiB) Viewed 8555 times

Image drawing routine:

Code: Select all

#define DrawTarget(x){ \
if(x.sleep) {}         \
else{                  \
CircleOut(x.X,x.Y,1);  \
CircleOut(x.X,x.Y,3);  \
CircleOut(x.X,x.Y,5);  \
    }                  \
}
#define WakeTarget(x){ \
x.sleep = false;               \
x.X = Random(99);       \
x.Y = Random(63);       \
}
bool sleep = false;
while(1){
    if(sleep)
    {
    sleep = false;
    SleepTarget(MyTarget);
    }
    else{sleep=1;}
    DrawTarget(MyTarget);
    Wait(1000);
    WakeTarget(MyTarget);
    }

Full Code:

Code: Select all

struct target{
    byte X;
    byte Y;
    bool sleep;
    unsigned long wakeTick;
};

#define DrawTarget(x){ \
if(x.sleep) {}         \
else{                  \
CircleOut(x.X,x.Y,1);  \
CircleOut(x.X,x.Y,3);  \
CircleOut(x.X,x.Y,5);  \
    }                  \
}

#define SleepTarget(x){                         \
x.sleep = true;                                 \
x.wakeTick = CurrentTick() + 2000;   \
}

#define WakeTarget(x){ \
x.sleep = false;               \
x.X = Random(99);       \
x.Y = Random(63);       \
}



task main(){
target MyTarget;
MyTarget.X = 30;
MyTarget.Y = 30;
MyTarget.sleep = 0;
bool slp = false;
while(1){
    if(slp)
    {
    slp = false;
    SleepTarget(MyTarget);
    }
    else{slp=1;}
    DrawTarget(MyTarget);
    Wait(1000);
    WakeTarget(MyTarget);
    }
}

[HaWe]

thank you Tim,
very impressive, this again proves: 1 picture says more than 1000 words (1 word= 4 byte)^^

After 5 Minutes the Mersenne Twister output looks like this:

TimGraphicMrand.jpg (28.67 KiB) Viewed 8544 times

(maybe be better compare number of loops instead of minutes in future)

Code: Select all

//-------------------------------------------------------
// Kernighan & Ritchie LCG (linear congruence generator)
//-------------------------------------------------------

unsigned long _RAND_SEED_ = 1;      // 1= default value (program start)
unsigned long _OLD_SEED_  = 1;


//--------------------------------------------

unsigned long Krand()               // custom random function
{
  _RAND_SEED_ = _RAND_SEED_ * 1103515245 + 12345;
  return (_RAND_SEED_ % (RAND_MAX + 1));
}

//--------------------------------------------


void sKrand(unsigned long seed)      // seeds for a new random series
{

  if (seed==0)                      // 0: a "real" randomized random seed
  {
    seed = ( CurrentTick()*BatteryLevel() );
  }  // substitute to time(0)

  else
  if (seed==-1) {                   // -1: restore last random series
    seed = _OLD_SEED_;
  }

  _OLD_SEED_  = seed;
  _RAND_SEED_ = seed;               // patch for rand_ function
}


//---------------------------------------------------------
// Mersenne Twister by Makoto Matsumoto & Takuji Nishimura
//---------------------------------------------------------
const int N_MTW = 25;
const int M_MTW = 7;
const unsigned long A_MTW[] = {0, 0x8ebfd028};

unsigned long y_MTW[25];
int i_MTW;
char virgin_MTW=TRUE, new_MTW=FALSE;

void InitMTW() {
  i_MTW = N_MTW+1;
}

//--------------------------------------------

unsigned long sMrand(unsigned long x) {
  unsigned long r, s;
  int k;
  if (x==1) InitMTW();

  else {
    InitMTW();
    if (x==0) {x=Random(CurrentTick());}
    for (k=0 ; k<25 ; ++k) {
      InitMTW();
      r = 9+x;
      s = 3402+x;
      for (k=0 ; k<25 ; ++k) {
           r = 509845221 * r + 3;
           s *= s + 1;
           y_MTW[k] = s + (r >> 10);
      }
      new_MTW=TRUE;
    }
  }
}

//--------------------------------------------

unsigned long Mrand(void) {
  unsigned long r, s, e;
  int k;

  if (virgin_MTW) {InitMTW(); virgin_MTW=FALSE; }

  if (i_MTW >= N_MTW) {
      if ((i_MTW > N_MTW) && !new_MTW) {
         r = 9;
         s = 3402;
         for (k=0 ; k<N_MTW ; ++k) {
           r = 509845221 * r + 3;
           s *= s + 1;
           y_MTW[k] = s + (r >> 10);
         }
      }
      for (k=0 ; k<N_MTW-M_MTW ; ++k)
         y_MTW[k] = y_MTW[k+M_MTW] ^ (y_MTW[k] >> 1) ^ A_MTW[y_MTW[k] & 1];
      for (; k<N_MTW ; ++k)
         y_MTW[k] = y_MTW[k+(M_MTW-N_MTW)] ^ (y_MTW[k] >> 1) ^ A_MTW[y_MTW[k] & 1];
      i_MTW = 0;
  }

  new_MTW=false;
  e = y_MTW[i_MTW++];
  e ^= (e << 7)  & 0x2b5b2500;
  e ^= (e << 15) & 0xdb8b0000;
  e ^= (e >> 16);
  return e;
}

//---------------------------------------------------------
// Grafik RNG Test
//---------------------------------------------------------

struct target{
    byte X;
    byte Y;
    bool sleep;
    unsigned long wakeTick;
};

#define DrawTarget(x){ \
if(x.sleep) {}         \
else{                  \
CircleOut(x.X,x.Y,1);  \
CircleOut(x.X,x.Y,3);  \
CircleOut(x.X,x.Y,5);  \
    }                  \
}

#define SleepTarget(x){                         \
x.sleep = true;                                 \
x.wakeTick = CurrentTick() + 2000;   \
}

#define WakeTarget(x){ \
x.sleep = false;       \
x.X = Mrand()%100;     \
x.Y = Mrand()% 63;     \
}



task main(){
  long time0;
  target MyTarget;
  MyTarget.X = 30;
  MyTarget.Y = 30;
  MyTarget.sleep = 0;
  bool slp = false;
  time0=CurrentTick();
  while((CurrentTick()-time0)<1000*60*5) {  // do it for 5 min

    if(slp)
    {
      slp = false;
      SleepTarget(MyTarget);
    }
    else{slp=1;}
    DrawTarget(MyTarget);
    Wait(1000);
    WakeTarget(MyTarget);
  }
  PlayTone(1000,100); Wait(200);
  PlayTone(1000,100); Wait(100);
  while(1);
}

[timpattinson]

1 picture says more than 1000 words (1 word= 4 byte)^^

Thanks for taking time to improve my code
EDIT: 3989/4
= 997.25
So close

[afanofosc]

Using the new rand() function with modulus (i.e., only change to Doc's code was to get rid of all that MRand stuff and call rand() instead):

nxtimage_12.gif (4.13 KiB) Viewed 8538 times

The new rand() function generates a single line of NBC code. The modulus adds one more. You can get a million rand()%16 loop iterations in 222900 ms.

Doc, when are you going to run the infinite monkey test using the new rand() function?

John Hansen

[HaWe]

As I already wrote, any Lehmer or K&R LCG would be better than the old Lego rand. That the old Lego fw rand() was really poor was why I opened this thread, and I tried to prove that.
The LCGs all are very well documented, so that this RNG issue is now clear and obvious so far, I'm convinced of that and I don't need to convince anybody of that.
Now to have a new rand function based on a LCG is a big leap for Legokind.

But LCGs of what kind ever do not meet the requirements of stochastical filters, as you probably know since you already read the references, too, as I assume ^^

[timpattinson]

When will the new rand() be available I have BricxCC 3.3.8.9 and the latest test release, plus the EFW that came with the test release/