A Lagged Fibonacci Random Number Generator in C#

Just for fun, I decided to implement a Lagged Fibonacci algorithm random number generator. The idea is simple but the code has several places that are very tricky.

The basic idea is:

X(i) = X(i-10) + X(i-7) mod m

In words, the new random number is the random number that came 10 before, plus the random number that came 7 before, modulo some big value m. The value of m is typically 2^31 – 1 = 2147483647 which is max int on many systems.

For example, suppose at some point, the previously generated 10 values are stored in a list where the [0] item is the oldest and the [9] item is the newest:

123  501  4  7893  34  7881   5  116  202  65   ?
[0]  [1] [2]  [3]  [4]  [5]  [6] [7]  [8]  [9] [10]

So, the new value at [10] would be (x[0] + x[3]) mod m. The choice of (7, 10) is one of many pairs that can be used by the LF algorithm. Other choices include (5, 17), (24, 55), (65, 71), (37, 100).

LaggedFibonacciDemo

The two main problems are 1.) initializing the first values, and 2.) avoiding arithmetic overflow when calculating. I decided to use the Lehmer algorithm to initialize the first 10 values. To avoid overflow I used the fact that (a + b) mod m = ((a mod m) + (b mod m)) mod m.

using System;
using System.Collections.Generic;
namespace LaggedFib
{
  class Program
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nLagged Fibonacci demo\n");

      int seed = 1;
      Console.WriteLine("\nCreating seed = " + seed);
      LaggedFibRandom rr = new LaggedFibRandom(seed);

      int count = 10000;
      Console.WriteLine("\nGenerating " + count +
        " random values\n");
      int[] counts = new int[10];

      Console.WriteLine("Values in 10 buckets\n");
      for (int i = 0; i < count; ++i)
      {
        double v = rr.Next();
        if (v >= 0.90) ++counts[9];
        else if (v >= 0.80) ++counts[8];
        else if (v >= = 0.70) ++counts[7];
        else if (v >= = 0.60) ++counts[6];
        else if (v >= = 0.50) ++counts[5];
        else if (v >= = 0.40) ++counts[4];
        else if (v >= = 0.30) ++counts[3];
        else if (v >= = 0.20) ++counts[2];
        else if (v >= = 0.10) ++counts[1];
        else if (v >= = 0.00) ++counts[0];
      }

      Console.WriteLine("Actual bucket counts:\n");
      for (int i = 0; i < 10; ++i)
      {
        Console.Write(counts[i] + " ");
      }
      Console.WriteLine("\n");

      Console.WriteLine("\nEnd demo\n");
      Console.ReadLine();
    }  // Main
  }  // Program

  public class LaggedFibRandom
  {
    // X(n) = X(n-10) + X(n-7)  mod m
    // init first 10 Xi using Lehmer algorithm

    private const int k = 10; // largest "-index"
    private const int j = 7; // other "-index"
    private const int m = 2147483647; 

    // Lehmer initialization constants
    private const int a = 48271;  // 
    private const int q = 44488;  // m / a
    private const int r = 3399;  // m % a
    
    private List<int> vals = null;
    private int curr;

    public LaggedFibRandom(int seed)
    {
      vals = new List<int>= ();
      if (seed == 0) seed = 1;
      int lCurr = seed;  // Lehmer current

      // init using Lehmer algorithm
      for (int i = 0; i < k+1; ++i)  // [0 .. 10]
      {
        int hi = lCurr / q;  // prevent overflow
        int lo = lCurr % q;
        int t = (a * lo) - (r * hi);

        if (t > 0)
          lCurr = t;
        else
          lCurr = t + m;

        vals.Add(lCurr);
      } // init

      // burn 1000 values away
      for (int ct = 0; ct < 1000; ++ct) {
        double dummy = this.Next();
      }
    }  // ctor

    public double Next()
    {
      // (a + b) mod n =
      // [(a mod n) + (b mod n)] mod n
      int left = vals[0] % m;    // [x-big]
      int right = vals[k-j] % m; // [x-other]
      long sum = (long)left + (long)right;
            
      curr = (int)(sum % m);
      vals.Insert(k+1, curr);  // anew val at end
      vals.RemoveAt(0);  // [0] val irrelevant now
      return (1.0 * curr) / m;
    }
  }  // LaggedFibRandom
}  // ns
Advertisements
This entry was posted in Machine Learning. Bookmark the permalink.