Linear congruential generator in C++

Question

I wrote a simple program (tried to implement the Linear congruential generator actually), but I'm not quite sure it works like it should.

I wanted to generate 250 number from [0,1] using my generator. However, it seems that instead of random numbers, I get equal values ..

How to improve it / what I did wrong?

Here's the code:

#include <iostream>
#include <cmath>

static const double A = 0.001342;
static const double C = 0.00025194;
static const double RAND_MAX = 1.0;

double rand()
{
    static double prev = 0;
    prev = A * prev + fmod(C, RAND_MAX);
    return prev;
}

int main(int argc, char **argv)
{
    for(int i=0; i<6; i++)
    std::cout << rand() << "\n";
    return 0;
}

And the output:

0.00025194
0.000252278
0.000252279
0.000252279
0.000252279
0.000252279

Switching to int instead of double, however gives some nice results:

#include <iostream>
#include <cmath>

static const int A = 5;
static const int C = 3;
static const int RAND_MAX = 8;

double rand()
{
    static int prev = 1;
    prev = A * prev + (C % RAND_MAX);
    return prev;
}

int main(int argc, char **argv)
{
    for(int i=0; i<100; i++)
    std::cout << rand() << "\n";
    return 0;
}

Output:

8
43
218
1093
5468
27343
136718
683593
3.41797e+06
1.70898e+07
8.54492e+07
4.27246e+08
2.13623e+09
2.09122e+09
1.86615e+09
7.40836e+08
-5.90786e+08
1.34104e+09
...

But I need it for generating random double numbers, greater than or equal to 0 and less than or equal to 1 :(

Answer 1

It's not the program, it's the choice of numbers.

prev is in the beginning equal to zero, so the first number becomes C.

Then, prev is equal to C, which makes prev A*C + C. However, A*C is so small, that when adding it as a floating point to the previous one, significant digits are shifted out and you're left with what you had before.

You can read more on What Every Computer Scientist Should Know About Floating-Point Arithmetic.

Answer 2

First, don't use floating point arithmetic for LCG's. Floating point is inherently imprecise, and can lead to undesirable behaviors such as fixed point convergence or interleaved short subcycles. With integer arithmetic, number theory can tell you what parameters are guaranteed to achieve full cycle, i.e., every value between 0 and M-1 (where M is your modulus) will be attained. See the Wikipedia article on LCGs for the full-cycle parameter requirements and a table of commonly used parameters.

Second, you have misinterpreted the LCG formula. It should be:

prev = (A * prev + C) % M;

Third, with your current choice of parameters you won't run into this but in general your intermediate calculations should be done with long to avoid overflow. The % operation will bring the answer back to an int, but if you stick with int arithmetic the multiplication may not yield the mathematically correct value which assures full cycle length and proper distributional behavior.