Un meilleur RNG…

En parcourant le web, j’ai lu de nombreuses critiques concernant les générateurs de nombres aléatoires livrés avec les systèmes d’exploitation, les langages de programmation…

Et notre couple Randomize/Rnd que vaut-il ?

Tester un RNG (Random Number Generator) nécessite de lui faire passer une batterie de tests statistiques qui permettent de mesurer leur qualité.
le plus célèbre jeu de tests est DieHards tests de G. Marsaglia mais je me suis amusé à tester notre générateur avec le programme de fourmilab (explications sur les résultats).

Résultats

Entropy = 6.960137 bits per byte.

Optimum compression would reduce the size
of this 11469916 byte file by 12 percent.

Chi square distribution for 11469916 samples is 91767249.86, and randomly
would exceed this value less than 0.01 percent of the times.

Arithmetic mean value of data bytes is 111.6449 (127.5 = random).
Monte Carlo value for Pi is 3.304316894 (error 5.18 percent).
Serial correlation coefficient is -0.036392 (totally uncorrelated = 0.0).

Les résultats démontrent que notre RNG n’est pas à la hauteur de nos espérances…

Si vous avez besoin d’un bon générateur de nombres aléatoires, vous avez deux possibilités :

Récupérer un code VBA existant sur le web
Depuis la fin des années 90, il existe d’excellents algorithmes qui utilisent nativement des variables ‘unsigned long’ voir des ‘unsigned long long’ ce qui ne facilite pas le portage en VBA avec notre type ‘long’…

J’ai quand même péché sur le web deux RNG de qualité.

  • Le fameux ‘Mersenne twister’ porté ici en VBA.
  • et ISAAC de qualité cryptographique

Ces deux RNG codés en VBA passent sans problème les tests qualité de Fourmilab.

Exemple avec ISAAC :

Entropy = 7.999985 bits per byte.

Optimum compression would reduce the size
of this 11469916 byte file by 0 percent.

Chi square distribution for 11469916 samples is 236.35, and randomly
would exceed this value 79.31 percent of the times.

Arithmetic mean value of data bytes is 127.4950 (127.5 = random).
Monte Carlo value for Pi is 3.141571792 (error 0.00 percent).
Serial correlation coefficient is -0.000191 (totally uncorrelated = 0.0).

En adapter un…
Je dis bien ‘adapter’ et non ‘programmer’ car se sont des algorithmes très pointus.
Je me suis d’abord penché sur des sources qui utilisaient des variables ‘uint32′ dans le but de simplifier l’adaptation.
Leur qualité dépassait rarement celle de notre Rnd et j’ai du me résoudre à explorer les sources ‘modernes’ qui exploitent les ‘uint64′.
Je me suis intéressé à la méthode MWC ‘Multiply with carry’ de G. Marsaglia tout particulièrement à l’algorithme MWC256.

Voici donc mon adaptation de son code :

'Random number generator - MWC256 - Multiply with Carry
'Author : George Marsaglia
' Reference : http://school.anhb.uwa.edu.au/personalpages/kwessen/shared/Marsaglia03.html
' Period about 2^8222. Uses a static array Q[256] and an initial carry 'c',
' the Q array filled with 256 random 32-bit integers in the calling program
' and an initial carry c<809430660 for the multiply-with-carry operation
'Adapted from C to VBA : Philben - v1.0 - 03/2012
Public Function MWC256_UINT32(Optional ByVal Start As Boolean = False) As Double
   Const ULL As Double = 2 ^ 64
   Const UL As Double = 2 ^ 32
   Static Q(1 To 256) As Variant, c As Variant, a As Variant, i As Integer
   Dim t As Variant, v As Variant

   If Start Or i = 0 Then
      Randomize

      'c < 809430660 - 0 <= Rnd < 1 - Marsaglia : c=362436
     c = CDec(CLng((809430660 - 2 ^ 12) * Rnd() + 2 ^ 11))
      a = CDec(809430660)

      For i = 1 To 256
         Q(i) = CLng(Rnd() * (2 ^ 31 - 1))
      Next i

      'warm up
     Do
         MWC256_UINT32
      Loop Until i = 256
   End If

   i = i + 1
   If i > 256 Then i = 1

   v = a * Q(i) + c
   t = CDec(v - Int(v / ULL) * ULL)
   c = CDec(Int(t / UL))
   Q(i) = v - Int(v / UL) * UL

   MWC256_UINT32 = Q(i)
End Function

'Signed integer : (-2^31 to 2^31-1) ?
'author : philben
Public Function MWC256_INT32(Optional ByVal Start As Boolean = False) As Long
   MWC256_INT32 = CLng(MWC256_UINT32(Start) - &H7FFFFFFF)
End Function

'Double : (0.0, 1.0) ?
'author : philben - v1.01
Public Function MWC256_Dbl(Optional ByVal Start As Boolean = False) As Double
   MWC256_Dbl = MWC256_UINT32(Start) * 2.32830643653869E-10
End Function

et le résultat des tests :

Entropy = 7.999985 bits per byte.

Optimum compression would reduce the size
of this 11469916 byte file by 0 percent.

Chi square distribution for 11469916 samples is 245.85, and randomly
would exceed this value 64.83 percent of the times.

Arithmetic mean value of data bytes is 127.4586 (127.5 = random).
Monte Carlo value for Pi is 3.144400759 (error 0.09 percent).
Serial correlation coefficient is -0.000303 (totally uncorrelated = 0.0).

Il semble que cette adaptation fonctionne bien et à l’avantage d’être plus rapide que les deux précédentes mais beaucoup plus lente que son équivalent en C et que Rnd…
Pour avoir un ordre d’idée, compter moins de 2 secondes par million de nombres générés.

Un autre inconvénient de mon portage est qu’il n’est validé que par moi (!) et je ne suis pas particulièrement fier de l’initialisation de l’algorithme…

@+

Philippe

Laisser un commentaire