random number issues

[mark]

Out of curiosity I grepped your codebase for usage of rand() and have found some problems there, you are probably unaware of.

There is one minor issue and a rather large one.

1. This is a minor concern, which is also probably known: rand() is not a mathematically state of the art random function, there are faster and "more random" generators like SFMT. Respecting the fact that gemrb runs on many platforms, it might still be valid to use rand() but as SFMT is written in plain C, I don't see a reason why it should not run on gemrb's platforms (it's not an external library,so there is no dependency either).
   
2. This is the major concern: In the codebase there are several occurences of formulas like 'pos = rand() % count' or 'time = 500 + 500 * (rand() % 20)' but these do not guarantee an equal outcome (sadly they are all over the net and quickly written as oneliners and seem so elegant). Especially in a dice game like the DnD series you should provide a random function with unbiased outcome.
   

Here are some details: You basically want a random number from a certain interval, e.g. [1,...,6] for a d6 but rand() delivers double values from 0 to RAND_MAX. So the outcome of rand() has to be mapped e.g. to your range of {1, ..., 6}. This can be done by modulo operation or by division with RAND_MAX and multiplication of some value - either way adds bias to some numbers from your range though.

 

This is why:

You divide the range from 0 to RAND_MAX into piles of a certain size (in this case 6) and whatever number comes out of rand(), you mapped it to a number in your range.

The problem is that this only works if the amount of rand() numbers is a multiple of your range, Otherwise there is a pile of rand() numbers that contains less numbers than the other piles, so the probability of chosing this pile is lower, which results into {1,...,6} numbers that are preferred in the total outcome.

 

RAND_MAX is mostly the prime number 2^32 -1, so if you're counting rand() % count, then count will never be a divisor or RAND_MAX and the outcome is skewed.

Let's say, rand() would create 2^32 possible integers (instead of doubles, but it comes to the same result), then rand % count would only be unbiased if count were a power of 2.

Let's compare a d6 and a d8 roll:

For RAND_MAX / 6 = 715827882.667 you are screwed, because the number of piles is not even; your d6 is flawed.

For RAND_MAX / 8 = 536870912 you are lucky because 2^3 = 8 fits well into RAND_MAX.

 

So, whatever RNG you use, you must not use the formulas you are using right now, because you cannot guaranty unbiased outcome.

 

One solution is to use rejection sampling, which means: divide your random numbers into piles of equal size like you want to do by intuition with your modulo operation - but: if one pile is smaller than the others and rand() delivers a number from this pile, reject it and ask for another one until you get such a number.

Then you have a uniform distribution of numbers and the game becomes fair.

 

If you want additional information on this, look at http://eternallyconf...w_art_rand.aspx and

http://www.azillionm...qed/random.html and other pages.

 

So alltogether I can really recommend SFMT as rand() replacement and strongly suggest to replace the skewed random formulas all over the codebase. Have a look at

https://github.com/V...mmon/rng_sfmt.h

https://github.com/V...on/rng_sfmt.cpp

There you will find a well commented implementation used in another opensource game, that uses rejection sampling and SFMT but you can also use the code without SFMT if you choose to do so.

So in the end you'll have a little RNG class where the constructor seeds your real rand function and one method that delivers an unbiased random number.

 

If you use C++11, you may also use the classes from <random>, where you instanciate one Mersenne Twister generator and a uniform_distribution object to address issue 1 and 2 at the same time. It is 2-6 times slower though than SFMT but still better than rand() and the skewed distribution functions.

[SyntaxError]

I'm sure we would happily accept a patch from you (I don't think its likely anybody else cares enough to do anything) with the following caveats:

 

1. any new dependancies should be optional

 

2. the code must still compile with older C++ compilers (I think C++ 98, not sure). Trust me I've been annoyed plenty by this requirement myself.

[lynx]

Seems a bit pointless to me, as the largest modulo rand call we have uses 1000 and the rest are much smaller. Comparing to the scale of RAND_MAX, the differences are probably negligable. And that 1000 is used just as a trivial rate-limiter, so the small skew is irrelevant.

 

Also, what did the original use? Probably plain srand and rand too (with some odd code for lucky rolls), but again, I doubt it matters.

 

Nice writeup though!

[Avenger]

i saw either this:

nNumber = rand() * nRange / RAND_MAX;

or this:

nNumber = rand()%nRange

I don't think we need better.

[mark]

i saw either this:

nNumber = rand() * nRange / RAND_MAX;

or this:

nNumber = rand()%nRange

I don't think we need better.

but those are the stochastically biased formulas, I would not use them.

I wrote a little demo code with those formulas and made some control samples: I generated several hundred thousand random numbers, counted the outcome and the bias was observable, some numbers came consistently some tenthousand times less often than the other numbers.

So probably every tenth roll or more is wrong.

I also found some information in the rand manpage which says some creepy things regarding rand():

 

The versions of rand() and srand() in the Linux C Library use the same random

number generator as random(3) and srandom(3), so the lower-order bits should be

as random as the higher-order bits. However, on older rand() implementations,

and on current implementations on different systems, the lower-order bits are

much less random than the higher-order bits. Do not use this function in appli‐

cations intended to be portable when good randomness is needed. (Use random(3)

instead.)

So depending on the system it's even worse because rand() itself is written poorly.

 

I'd like to write a patch and include some SFMT files, but I have no overview of the project; some parts seem to be scripted, others not and all sources must use the new RNG if we want improvement. I did not find e.g. the dice rolling code itself (just picking something that interested me in the first place).

If you can give me a short info about where you'd like to have the new code and where the most important things are and what I have to know about the config of your build system (all of which might need changes for the new stuff), then I can start.

 

Edit: oh, and I need to know if some threading happens somewhere

[AnonymousHero]

So depending on the system it's even worse because rand() itself is written poorly.

The systems mentioned are all but obsolete, and so aren't very relevant. However, I'm sure any attempt at moving towards unbiased "random-number-between-0-and-N" procedures (etc.) would be welcomed. I don't see how that would need a new implementation of rand(), given the manual page, though.

 

(Just an interested observer.)

[lynx]

Avenger was replying to me, saying what he saw while reverse engineering the original, not what we have here.

You're still hypotethising about the errors and see what I wrote above — we're not using arbitrary moduli(?) often and especially in the game rule case, it is eventually mostly the dice rolls you'd expect.

 

Can you give us a table or graph of deviations from the expected probabilities for eg. 1d10 and 1d100 (for a 32-bit RAND_MAX and large sample size)? It should prove whether this whole thing is "just" an academic endeavour or not. It would be a shame for you to waste time if it turns out the differences are negligible.

 

That rand caveat was likely written in the 90's. Most of the dice-like randomness goes through Interface::Roll. Threading is only present for the audio plugins so far.

[mark]

it does not say "on old systems" it says "on old systems and different systems". Whatever that means, it does not imply that all current systems are safe.

[lynx]

Oh come on. Looking at the repo, the 2004 import already has this mentioned, so it's been there unchanged for at least 10 years:

http://git.kernel.org/cgit/docs/man-pages/man-pages.git/log/man3/rand.3

 

Judging from the comments, it was added in 1995, so my guess seemed to be spot on.

[mark]

regarding the wasting time concern - I already imlemented the SFMT wrapperclass with uniform distribution function some days ago and used half a day for testing the RNG quality, so that's not really an issue

The only thing left would be placing the code somewhere in the codebase and replace the rand() calls if you want it.

 

Currently I finished evaluating a chi^2 test.

I wanted to do this for all numbers from 0 to RAND_MAX at first but unfortunately this will require more than 4GB of RAM only for the random numbers, which I don't have available, and I wasn't able to link against STXXL (http://stxxl.sourceforge.net/). So I chi-squared the numbers produced by the formulas rand%size and sfmt(0,size-1) for different sizes.

 

Here are some results:

Chi^2 value pairs for rand%size, sfmt(0,size-1) for d{2,3,6,8,10,12,20,100}; p=10%, 90% in brackets.

The results should be in between, everything lower or higher is suspicious. Matching p=50% would be pretty perfect.

To compare for yourself, search for "Selected percentage points of the chi-square distribution" tables, take the chi^2 value and find it in the d{k-1} line. Or use http://statpages.org/pdfs.html

 

As a rule of thumb (according to Knuth) you can say

<1%, >99% not random

1-5, 95-99 suspect

5-10, 90-95 almost suspect

Everything beyond those values is marked by * and commented below (p has been computed for these values to tell how bad the result is).

 

The first 4 value pairs were for 1,000,000,000 samples each

The 5th for 100,000,000

The 6-7th for 500,000,000

 

Some hours later:

 

d2 (0.0158, 2.7055)

0.0363609, 0.297287

1.52194, 1.25628

0.661518, 0.637664

0.118287, 0.751089

1.41515, 0.321716

4.17808*, 1.4526

0.000245*, 0.632399

Two cases where rand is beyond 95% or below 2%.

 

d3 (0.2107, 4.6052)

2.1746, 0.7441

2.48713, 4.01729

2.23582, 1.50578

2.82781, 1.92041

5.64079*, 0.74178

0.249783, 1.19496

0.405172, 0.644316

One case where rand is beyond 95%

 

d6 (1.6103, 9.2364)

5.48917, 2.89016

14.5876*, 9.12873

4.98174, 6.13206

2.58365, 0.961907*

10.7035*, 11.8858*

3.24196, 3.01078

0.894795*, 3.32388

rand is above 99%, 95%, below 3%

sfmt is above 97%, below 4%

The ideal chi^2 should be near 4.3515:

6.06892, 5.33.. (average)

 

d8 (2.8331, 12.0170)

18.2556*, 6.91943

9.59598, 5.24599

4.86716, 16.8013*

8.61307, 3.94753

4.18284, 4.99872

5.95096, 3.44947

14.4946*, 8.16555

rand is above 99%, 96%

sfmt is above 98%

 

d10 (4.1682, 14.6837)

4.90288, 7.87709

6.97481, 5.71335

7.55181, 6.64293

5.94711, 7.61737

9.98804, 12.4981

7.41316, 6.91671

4.24675, 14.6781

Nothing suspicious here. The ideal chi^2 should be near 8.3428:

6.71779, 8.84909 (averages)

 

d12 (5.5778, 17.2750)

9.71548, 9.25721

12.5145, 13.315

6.89115, 10.0049

20.0306*, 16.0564

9.99703, 9.19753

5.50514*, 14.4397

13.6396, 16.516

rand is above 95%, below 10%

 

d20 (11.6509, 27.2036)

21.188, 12.6777

16.4749, 25.2602

17.7793, 29.2486*

29.2415*, 23.3705

19.988, 19.7884

22.455, 21.3371

24.4157, 21.8703

Both beyond 94%

 

d100 (81.4493, 117.4069)

93.205, 80.1217

80.824*, 115.596

88.8697, 93.8696

94.7682, 79.6047*

109.524, 91.0377

81.8716, 84.232

114.466, 130.344*

rand is below 10%

sfmt is below 10%, beyond 98%

Optimal chi^2 should be around 98.3341:

94.78979, 96.40081 (averages)

 

There are some results where both RNGs produced suspect numbers, but rand produced twice as much of those and in two cases rand's numbers were even (statistically) nonrandom. In the three cases where I calculated the average distance to p=50%, rand was also worse.