Générateur de Nombres Aléatoires "Mersenne Twister"

[grabiller]

Hello,

vu que j'en ai besoin pour mes expérimentations avec smallpt, j'ai porté un générateur de nombre aléatoire 'Mersenne Twister' basé sur la librairie Boruvka.

Ce genre de code est difficile à tester, d'être sûr qu'il marche correctement.

Je le poste ici au cas où vous y détecteriez des erreurs, pour que vous puissiez le tester vous-même, voir si il marche comme attendu.

Je ne suis pas encore trop sûr de comment créer un code de test efficace permettant de dire si il marche réellement selon l'algorithme. Si vous avez des suggestions n'hésitez pas.

Code : Tout sélectionner

; ============================================================================
;
;  Copyright (c) 2013, Guy Rabiller, RADFAC.
;  All rights reserved, worldwide.
;
;  Redistribution  and  use  in  source  and  binary  forms,  with or  without
;  modification, are permitted provided that the following conditions are met:
;
;  - Redistributions of  source code  must retain  the above copyright notice,
;    this list of conditions and the following disclaimer.
;  - Redistributions in binary form must reproduce the above copyright notice,
;    this list of conditions and the following disclaimer in the documentation
;    and/or other materials provided with the distribution.
;  - Neither the name of  RADFAC nor the names of its contributors may be used
;    to  endorse  or  promote  products  derived  from  this  software without
;    specific prior written permission.
;
;  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
;  AND ANY EXPRESS OR IMPLIED WARRANTIES,  INCLUDING,  BUT NOT LIMITED TO, THE
;  IMPLIED WARRANTIES OF  MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
;  ARE DISCLAIMED.  IN NO EVENT SHALL THE  COPYRIGHT HOLDER OR CONTRIBUTORS BE
;  LIABLE  FOR  ANY  DIRECT,  INDIRECT,  INCIDENTAL,  SPECIAL,  EXEMPLARY,  OR
;  CONSEQUENTIAL  DAMAGES  (INCLUDING,  BUT  NOT  LIMITED  TO,  PROCUREMENT OF
;  SUBSTITUTE GOODS OR SERVICES;  LOSS OF USE,  DATA,  OR PROFITS; OR BUSINESS
;  INTERRUPTION) HOWEVER CAUSED  AND ON ANY  THEORY OF  LIABILITY,  WHETHER IN
;  CONTRACT,  STRICT LIABILITY,  OR TORT  (INCLUDING  NEGLIGENCE OR OTHERWISE)
;  ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,  EVEN IF ADVISED OF THE
;  POSSIBILITY OF SUCH DAMAGE.
;
;  For permission, contact copy@radfac.com.
;
; ============================================================================
;- raafal.rand.mersenne.twister.pbi
; ............................................................................
;- raafal Mersenne Twister Random Generator Class.
; ============================================================================
;  2013/08/27 | Guy Rabiller
;  - creation
;  - based on https://github.com/danfis/boruvka
; ============================================================================
; RandMT - Mersenne Twister Random Number Generator
; ----------------------------------------------------------------------------
; 
; The Mersenne Twister is an algorithm for generating random numbers.  It
; was designed with consideration of the flaws in various other generators.
; The period, 2^19937-1, and the order of equidistribution, 623 dimensions,
; are far greater.  The generator is also fast; it avoids multiplication and
; division, and it benefits from caches and pipelines.  For more information
; see the inventors' web page at
; http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
; 
; Reference:
; M. Matsumoto and T. Nishimura, "Mersenne Twister: A 623-Dimensionally
; Equidistributed Uniform Pseudo-Random Number Generator", ACM Transactions on
; Modeling and Computer Simulation, Vol. 8, No. 1, January 1998, pp 3-30.
; ============================================================================
EnableExplicit

; ============================================================================
;- Class (CRandMT)
; ============================================================================
;{
; ----------------------------------------------------------------------------
;- Public
; ----------------------------------------------------------------------------
;{
DeclareModule CRandMT
  ; ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  ;  Interface
  ; ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  Interface Instance
    ; ...[ Object ]...........................................................
    Delete()
    ; ...[ Seeds ]............................................................
     ; Reseeds generator.
    Reseed( seed.l )
    ; Similar to *Reseed() but more seeds are used.
    Reseed2( *seed.Long, seedlen.l )
    ; Similar to *Reseed() function but seed is chosen automatically from
    ; /dev/urandom if available or from time() and clock().
    ReseedAuto()
    ; ...[ Generate ].........................................................
    ; Returns randomly generated number in range [from, upto).
    Rand         .d ( from.d, upto.d )
    ; Returns number between [0-1) real interval.
    Rand01       .d (                )
    ; Returns number between [0-1] real interval.
    Rand01Closed .d (                )
    ; Returns number between (0-1) real interval.
    Rand01Open   .d (                )
    ; Returns number between [0-1) real interval with 53-bit resolution.
    Rand01_53    .d (                )
    ; Returns random integer number in interval [0, 2^32-1].
    RandInt      .l (                )
    ; Returns number from a normal (Gaussian) distribution.
    RandNormal   .d (                )
  EndInterface
  ; ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  ;  Constructors
  ; ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  Declare.i New( seed.l )
  Declare.i New2( *seed.Long, seedlen.l )
  Declare.i NewAuto()
EndDeclareModule
;}
; ----------------------------------------------------------------------------
;- Private
; ----------------------------------------------------------------------------
;{
Module CRandMT
  ; ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  ;  Constants
  ; ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  ;{
  #BOR_RAND_MT_N = 624 ; Length of state vector
  #BOR_RAND_MT_M = 397 ; Period
  #MmN           = #BOR_RAND_MT_M - #BOR_RAND_MT_N
  #S32_MAX       = 2147483647
  ;}
  ; ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  ;  Structures
  ; ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  ;{
  Structure Longs
    l.l[0]
  EndStructure
  Structure big_seed_t
    l.l[#BOR_RAND_MT_N]
  EndStructure
  ;}
  ; ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  ;  Instance
  ; ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  ;{
  Structure Instance_t
    *vtable
    state .l[#BOR_RAND_MT_N] ; Internal state
    *next .Long              ; Next value from state[]
    left  .l                 ; Number of values left before reload needed
  EndStructure
  ;}
  ; ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  ;  Macros (Private)
  ; ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  ;{
  Macro hiBit( u_ )
    ( u & $80000000 )
  EndMacro
  Macro loBit( u_ )
    ( u & $00000001 )
  EndMacro
  Macro loBits( u_ )
    ( u & $7fffffff )
  EndMacro
  Macro mixBits( u_, v_ )
    ( hiBit(u_) | loBits(v_) )
  EndMacro
  Macro twist( m_, s0_, s1_ )
    ( m_ ! (mixBits(s0, s1) >> 1) ! magic(s1) )
  EndMacro
  ;}
  ; ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  ;  Implementation (Private)
  ; ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  ;{
  Procedure.q BOR_MAX( x.q, y.q )
    If x > y
      ProcedureReturn x
    Else
      ProcedureReturn y
    EndIf
  EndProcedure
  Procedure.q magic( u.q )
    ; loBit(u) ? 0x9908b0dfUL : 0x0UL
    If loBit(u)
      ProcedureReturn $9908b0df
    Else
      ProcedureReturn $0
    EndIf
  EndProcedure
  Procedure Init( *Me.Instance_t, seed.q )
    ; Initialize generator state with seed
    ; See Knuth TAOCP Vol 2, 3rd Ed, p.106 for multiplier.
    ; In previous versions, most significant bits (MSBs) of the seed affect
    ; only MSBs of the state array.  Modified 9 Jan 2002 by Makoto Matsumoto.
    
    ; ---[ Local Variables ]--------------------------------------------------
    Protected i.l=1,j.l=0,q.q
    
    ; ---[ Init ]-------------------------------------------------------------
    *Me\state[0] = seed & $ffffffff
    While i < #BOR_RAND_MT_N
      *Me\state[i] = ( 1812433253*( *Me\state[j] ! (*Me\state[j] >> 30) ) + i ) & $ffffffff
       i + 1 : j + 1
    Wend
    
  EndProcedure
  Procedure Reload( *Me.Instance_t )
    ; Generate N new values in state
    ; Made clearer and faster by Matthew Bellew (matthew.bellew@home.com)
    
    ; ---[ Local Variables ]--------------------------------------------------
    Protected *p .Long  = @*Me\state[0]
    Protected *p_.Longs = *p
    Protected  i.l
    
    ; ---[ Reload ]-----------------------------------------------------------
    i = #BOR_RAND_MT_N - #BOR_RAND_MT_M
    While i >= 0
      *p\l = twist( *p_\l[#BOR_RAND_MT_M], *p_\l[0], *p_\l[1] )
      i - 1 : *p  + SizeOf(Long) : *p_ = *p
    Wend
    i = #BOR_RAND_MT_M
    While i > 0
      *p\l = twist( *p_\l[#MmN], *p_\l[0], *p_\l[1] )
      i - 1 : *p  + SizeOf(Long) : *p_ = *p
    Wend
    *p\l = twist( *p_\l[#MmN], *p_\l[0], *Me\state[0] )
    
    *Me\left = #BOR_RAND_MT_N
    *Me\next = @*Me\state[0]

  EndProcedure
  ;}
  ; ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  ;  Implementation
  ; ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  ;{
  ; ~~~[ Object ]~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  Procedure Delete( *Me.Instance_t )
    
    ; ---[ Release Instance ]-------------------------------------------------
    FreeMemory(*Me)
    
  EndProcedure
  ; ~~~[ Seeds ]~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  Procedure Reseed( *Me.Instance_t, seed.l )
    
    Init  ( *Me, seed )
    Reload( *Me )
    
  EndProcedure
  Procedure Reseed2( *Me.Instance_t, *seed.Long, seedlen.l )
    
    ; Seed the generator with an array of uint32's.
    ; There are 2^19937-1 possible initial states. This function allows
    ; all of those to be accessed by providing at least 19937 bits (with a
    ; default seed length of N = 624 uint32's). Any bits above the lower 32
    ; in each element are discarded.
    
    ; ---[ Local Variables ]--------------------------------------------------
    Protected i.l=1,j.l=0,k.l=BOR_MAX(#BOR_RAND_MT_N,seedlen)
    Protected *seed_.Longs = *seed
    
    Init( *Me, 19650218 )
    
    While k
      *Me\state[i] ! ( ( *Me\state[i-1] ! ( *Me\state[i-1] >> 30 ) )*1664525 )
      *Me\state[i] + ( ( *seed_\l[j] & $ffffffff ) + j )
      *Me\state[i] & $ffffffff
      i + 1
      j + 1
      If i >= #BOR_RAND_MT_N
        *Me\state[0] = *Me\state[#BOR_RAND_MT_N - 1]
        i = 1
      EndIf
      If j >= seedlen
        j = 0
      EndIf
      k - 1
    Wend
    
    k = #BOR_RAND_MT_N - 1
    While k
      *Me\state[i] ! ( ( *Me\state[i-1] ! ( *Me\state[i-1] >> 30 ) )*1566083941 )
      *Me\state[i] - i
      *Me\state[i] & $ffffffff
      i + 1
      If i >= #BOR_RAND_MT_N
        *Me\state[0] = *Me\state[#BOR_RAND_MT_N - 1]
        i = 1
      EndIf
      k - 1
    Wend
    
    *Me\state[0] = $80000000 ; MSB is 1, assuring non-zero initial array.

    Reload( *Me )
    
  EndProcedure
  Procedure ReseedAuto( *Me.Instance_t )
    
    ; ---[ Local Variables ]--------------------------------------------------
    Protected bigSeed.big_seed_t,i.l=#BOR_RAND_MT_N
    
    ; ---[ Generate Big Seed ]------------------------------------------------
    While i > 0
      i - 1
      bigSeed\l[i] = Random(#S32_MAX,0)
    Wend
    
    ; ---[ Reseed ]-----------------------------------------------------------
    Reseed2( *Me, @bigSeed\l[0], #BOR_RAND_MT_N )
    
  EndProcedure
  ; ~~~[ Generate ]~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  Procedure.d Rand( *Me.Instance_t, from.d, upto.d )
    
    ; ---[ Retrieve Interface ]-----------------------------------------------
    Protected Me.Instance = *Me
    
    ; ---[ Generate ]---------------------------------------------------------
    Protected val.d = Me\Rand01()
    val * ( upto - from )
    val + from
    ProcedureReturn val
    
  EndProcedure
  Procedure.d Rand01( *Me.Instance_t )
    
    ; ---[ Retrieve Interface ]-----------------------------------------------
    Protected Me.Instance = *Me
    
    ; ---[ Generate ]---------------------------------------------------------
    CompilerIf #PB_Compiler_Processor = #PB_Processor_x86
      ProcedureReturn Me\RandInt()*(1.0/4294967296.0)+0.5
    CompilerElse
      ProcedureReturn Me\RandInt()*(1.0/4294967296.0)
    CompilerEndIf
    
  EndProcedure
  Procedure.d Rand01Closed( *Me.Instance_t )
    
    ; ---[ Retrieve Interface ]-----------------------------------------------
    Protected Me.Instance = *Me
    
    ; ---[ Generate ]---------------------------------------------------------
    CompilerIf #PB_Compiler_Processor = #PB_Processor_x86
      ProcedureReturn Me\RandInt()*(1.0/4294967295.0)+0.5
    CompilerElse
      ProcedureReturn Me\RandInt()*(1.0/4294967295.0)
    CompilerEndIf
    
  EndProcedure
  Procedure.d Rand01Open( *Me.Instance_t )
    
    ; ---[ Retrieve Interface ]-----------------------------------------------
    Protected Me.Instance = *Me
    
    ; ---[ Generate ]---------------------------------------------------------
    CompilerIf #PB_Compiler_Processor = #PB_Processor_x86
      ProcedureReturn ( Me\RandInt() + 0.5 )*(1.0/4294967296.0)+0.5
    CompilerElse
      ProcedureReturn ( Me\RandInt() + 0.5 )*(1.0/4294967296.0)
    CompilerEndIf
    
  EndProcedure
  Procedure.d Rand01_53( *Me.Instance_t )
    
    ; ---[ Retrieve Interface ]-----------------------------------------------
    Protected Me.Instance = *Me
    
    ; ---[ Generate ]---------------------------------------------------------
    Protected a.q = Me\RandInt() >> 5
    Protected b.q = Me\RandInt() >> 6
    CompilerIf #PB_Compiler_Processor = #PB_Processor_x86
      ProcedureReturn ( a*67108864.0 + b )*( 1.0/9007199254740992.0 )+0.5 ; by Isaku Wada
    CompilerElse
      ProcedureReturn ( a*67108864.0 + b )*( 1.0/9007199254740992.0 ) ; by Isaku Wada
    CompilerEndIf
    
  EndProcedure
  Procedure.q RandInt( *Me.Instance_t )
    ; Pull a 32-bit integer from the generator state.
    ; Every other access function simply transforms the numbers extracted here.
    
    ; ---[ Local Variable ]---------------------------------------------------
    Protected s1.q
    
    ; ---[ Check Reload ]-----------------------------------------------------
    If *Me\left = 0
      Reload( *Me )
    EndIf
    *Me\left - 1
    
    ; ---[ Generate ]---------------------------------------------------------
    s1 = *Me\next\l : *Me\next + SizeOf(Long)
    s1 ! (  s1 >> 11 )
    s1 ! ( (s1 <<  7) & $9d2c5680 )
    s1 ! ( (s1 << 15) & $efc60000 )
    s1 ! (  s1 >> 18 )
    ProcedureReturn s1
    
  EndProcedure
  Procedure.d RandNormal( *Me.Instance_t, mean.d, stddev.d )
    ; Return a real number from a normal (Gaussian) distribution with given
    ; mean and standard deviation by polar form of Box-Muller transformation.
    
    ; ---[ Retrieve Interface ]-----------------------------------------------
    Protected Me.Instance = *Me
    
    ; ---[ Local Variables ]--------------------------------------------------
    Protected x.d,y.d,r.d,s.d
    
    ; ---[ Generate ]---------------------------------------------------------
    Repeat
      x = 2.0*Me\Rand01() - 1.0
      y = 2.0*Me\Rand01() - 1.0
      r = x*x + y*y
    Until r < 1.0 And r <> 0.0
    s = Sqr( -2.0 * Log(r) / r )
    ProcedureReturn mean + x*s*stddev
    
  EndProcedure
  ;}
  ; ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  ;  VTable
  ; ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  ;{
  DataSection
    VRandMT:
    ; ---[ Object ]-----------------------------------------------------------
    Data.i @Delete()
    ; ---[ Seeds ]------------------------------------------------------------
    Data.i @Reseed()
    Data.i @Reseed2()
    Data.i @ReseedAuto()
    ; ---[ Generate ]---------------------------------------------------------
    Data.i @Rand()
    Data.i @Rand01()
    Data.i @Rand01Closed()
    Data.i @Rand01Open()
    Data.i @Rand01_53()
    Data.i @RandInt()
    Data.i @RandNormal()
  EndDataSection
  ;}
  ; ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  ;  Constructors
  ; ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  ;{
  Procedure.i Nes( *Me.Instance_t, seed.l )
    
    ; ---[ Sanity Check ]-----------------------------------------------------
    If #Null = *Me : ProcedureReturn #Null : EndIf
    
    ; ---[ Init VTable ]------------------------------------------------------
    *Me\vtable = ?VRandMT
    
    ; ---[ Set Members ]------------------------------------------------------
    Reseed( *Me, seed )
    
    ; ---[ Return Initialized Object ]----------------------------------------
    ProcedureReturn *Me
    
  EndProcedure
  Procedure.i New( seed.l )
    
    ; ---[ Allocate Instance Memory ]-----------------------------------------
    Protected *p.Instance_t = AllocateMemory(SizeOf(Instance_t))
    
    ; ---[ Init Instance ]----------------------------------------------------
    ProcedureReturn Nes( *p, seed )
    
  EndProcedure
  Procedure.i Nes2( *Me.Instance_t, *seed.Long, seedlen.l )
    
    ; ---[ Sanity Check ]-----------------------------------------------------
    If #Null = *Me : ProcedureReturn #Null : EndIf
    
    ; ---[ Init VTable ]------------------------------------------------------
    *Me\vtable = ?VRandMT
    
    ; ---[ Set Members ]------------------------------------------------------
    Reseed2( *Me, *seed, seedlen )
    
    ; ---[ Return Initialized Object ]----------------------------------------
    ProcedureReturn *Me
    
  EndProcedure
  Procedure.i New2( *seed.Long, seedlen.l )
    
    ; ---[ Allocate Instance Memory ]-----------------------------------------
    Protected *p.Instance_t = AllocateMemory(SizeOf(Instance_t))
    
    ; ---[ Init Instance ]----------------------------------------------------
    ProcedureReturn Nes2( *p, *seed, seedlen )
    
  EndProcedure
  Procedure.i NesAuto( *Me.Instance_t )
    
    ; ---[ Sanity Check ]-----------------------------------------------------
    If #Null = *Me : ProcedureReturn #Null : EndIf
    
    ; ---[ Init VTable ]------------------------------------------------------
    *Me\vtable = ?VRandMT
    
    ; ---[ Set Members ]------------------------------------------------------
    ReseedAuto( *Me )
    
    ; ---[ Return Initialized Object ]----------------------------------------
    ProcedureReturn *Me
    
  EndProcedure
  Procedure.i NewAuto()
    
    ; ---[ Allocate Instance Memory ]-----------------------------------------
    Protected *p.Instance_t = AllocateMemory(SizeOf(Instance_t))
    
    ; ---[ Init Instance ]----------------------------------------------------
    ProcedureReturn NesAuto( *p )
    
  EndProcedure
  ;}
EndModule
;}
;}

; ============================================================================
;- TEST CODE (Mainfile)
; ============================================================================
;{
CompilerIf #PB_Compiler_IsMainFile
  
  Procedure.q getULong( *source.Long )
    ; from: http://code.google.com/p/purebasic-extension-for-adobe-air/source/browse/trunk/pureair/native/src/Unsigned.pb
    ;- Reads 4 bytes from the specified memory address,
    ; and returns the value as *unsigned* integer
    ; (minimum = 0, maximum = 4294967295).
    If *source\l < 0
      ProcedureReturn *source\l + $100000000
    Else
      ProcedureReturn *source\l
    EndIf
  EndProcedure
  Procedure.q fromULong( source.l )
    ; ; from: http://code.google.com/p/purebasic-extension-for-adobe-air/source/browse/trunk/pureair/native/src/Unsigned.pb
    ProcedureReturn getULong(@source)
  EndProcedure

  Procedure test( r.CRandMT::Instance )
    Protected i
    
    ;Debug "    Generate 10 numbers in [0,1):"
    ;For i=0 To 9
      ;Debug " "+StrD(r\Rand01())
    ;Next
    ;Debug ""
    
    Debug "    Generate 10 numbers in [0,1]:"
    For i=0 To 9
      Debug " "+StrD(r\Rand01Closed())
    Next
    Debug ""
    
    ;Debug "    Generate 10 numbers in (0,1):"
    ;For i=0 To 9
      ;Debug " "+StrD(r\Rand01Open())
    ;Next
    ;Debug ""
    
    Debug "    Generate 5 integers in:"
    For i=0 To 4
      Debug " "+StrU(fromULong(r\RandInt()))
    Next
    Debug ""
    
    ;Debug "    Generate 10 numbers in [0,1.5):"
    ;For i=0 To 9
      ;Debug " "+StrD(r\Rand(0.0,1.5))
    ;Next
    ;Debug ""
    
    ;Debug "    Generate 10 numbers in [1.1,1.2):"
    ;For i=0 To 9
      ;Debug " "+StrD(r\Rand(1.1,1.2))
    ;Next
    ;Debug ""
    
    ;Debug "    Generate 10 numbers in [-0.5,0.2):"
    ;For i=0 To 9
      ;Debug " "+StrD(r\Rand(-0.5,0.2))
    ;Next
    ;Debug ""
    
  EndProcedure
  
  Global rmt.CRandMT::Instance
  
  rmt = CRandMT::NewAuto()
  Global rMax.d = 0.0
  Global rMin.d = 1.0
  Global i.l,rr.d
  For i = 0 To 6294
    rr = rmt\Rand01()
    If rr > rMax : rMax = rr : EndIf
    If rr < rMin : rMin = rr : EndIf
  Next
  Debug "rMax: "+ StrD(rMax)
  Debug "rMin: "+ StrD(rMin)
  rmt\Delete()
  
  Debug "Auto:"
  rmt = CRandMT::NewAuto()
  test( rmt )
  rmt\Delete()
  
  Debug "7654:"
  rmt = CRandMT::New(7654)
  test( rmt )
  rmt\Delete()
  
  Debug "1111:"
  rmt = CRandMT::New(1111)
  test( rmt )
  
  Debug "1111 -> 2222:"
  rmt\Reseed(2222)
  test( rmt )
  
  Debug "2222 -> auto:"
  rmt\ReseedAuto()
  test( rmt )
  rmt\Delete()

CompilerEndIf
;}

; ============================================================================
;  EOF
; ============================================================================

[SPH]

Pourrais tu m'en dire plus sur la routine mersenne twister ? (a quoi ca sert notamment)

[grabiller]

Pourrais tu m'en dire plus sur la routine mersenne twister ? (a quoi ca sert notamment)

Mersenne Twister Wikipedia

En résumé, je cite:
../.. Développé par Makoto Matsumoto et Takuji Nishimura en 1997, le Mersenne Twister est un générateur de nombres pseudo-aléatoires particulièrement réputé pour sa qualité ../..
../.. L'algorithme Mersenne Twister a été optimisé pour être utilisé dans le cadre de simulations de Monte-Carlo dans un grand nombre de domaines ../..

Quant à son algorithme lui-même, sur un plan mathématique, je serai bien incapable de te l'expliquer

[SPH]

Pourrais tu m'en dire plus sur la routine mersenne twister ? (a quoi ca sert notamment)

Mersenne Twister Wikipedia

En résumé, je cite:
../.. Développé par Makoto Matsumoto et Takuji Nishimura en 1997, le Mersenne Twister est un générateur de nombres pseudo-aléatoires particulièrement réputé pour sa qualité ../..
../.. L'algorithme Mersenne Twister a été optimisé pour être utilisé dans le cadre de simulations de Monte-Carlo dans un grand nombre de domaines ../..

Quant à son algorithme lui-même, sur un plan mathématique, je serai bien incapable de te l'expliquer

J'y suis deja allé sur wikipedia. Mais comme je n'ai rien compris, je te demande d'expliquer avec des mots de tous les jours...

[grabiller]

Ok, mais que ne comprends-tu pas dans:

../.. Développé par Makoto Matsumoto et Takuji Nishimura en 1997, le Mersenne Twister est un générateur de nombres pseudo-aléatoires particulièrement réputé pour sa qualité ../..
../.. L'algorithme Mersenne Twister a été optimisé pour être utilisé dans le cadre de simulations de Monte-Carlo dans un grand nombre de domaines ../..

Que veux-tu savoir de plus, exactement ?

[SPH]

c'est surtout ce que tu vas en faire qui m'interroge

[G-Rom]

c'est surtout ce que tu vas en faire qui m'interroge

Pour le raytracing, c'est un générateur de nombre aléatoire.
Le problème de l'algorithme qu'il utilise , c'est que l'image n'est jamais à 100% fini , en gros son raytracer lance des rayons n'importe ou ( random )
la fonction de random est donc importante pour son problème (réduction du grain) , d'ou l'importation de cet algo je suppose.
Ceci dit, PB est codé en C/C++ , l'algo si dessus est standard d’après ce que j'ai pu lire sur wiki , je suppose que PB aussi du coup.

[comtois]

J'ai récupéré ton code sur le forum anglais pour ne plus avoir l'erreur avec setULong().
Maintenant j'ai une erreur à la 299
Random() : Max ne peut pas être négatif.
En effet la doc indique que la valeur maxi est 2147483647, et tu as #U32_MAX = 4294967295

C'est bizarre que tu n'aies pas cette erreur. Testé sous 5.20 beta 14

[grabiller]

J'ai récupéré ton code sur le forum anglais pour ne plus avoir l'erreur avec setULong().
Maintenant j'ai une erreur à la 299
Random() : Max ne peut pas être négatif.
En effet la doc indique que la valeur maxi est 2147483647, et tu as #U32_MAX = 4294967295
C'est bizarre que tu n'aies pas cette erreur. Testé sous 5.20 beta 14

Je viens de corriger le code ici aussi.

Oui effectivement, c'est bizarre, je n'ai pas cette erreur.

Peux-tu m'en dire plus sur ta config ? Tu es sous Windows ? x86/x64 ?

Je n'ai pas encore tester le code avec la beta 14, je vais essayer aujourd'hui. Avec la beta 13 en tous cas, je n'ai pas cette erreur. Par contre je viens de découvrir un bug en 32 bits et avec les structures du type "Long". Je vais checker sur la b14 aussi.

Merci pour ton retour.

[grabiller]

c'est surtout ce que tu vas en faire qui m'interroge

Et bien effectivement, comme le dit G-Rom, la qualité - ainsi que la rapidité - d'un générateur de nombres aléatoires est important dans le cas d'un path tracer "Monte Carlo".

Car l'idée de cet algorithme est de lancer des rayons aléatoirement, la moyenne de leurs résultats donnant la solution idéale (mais il faut un grand nombre de rayons pour que cela marche).

Or le seul défaut dans l'image avec un path tracer 'brute force' (sans optimisation ou techniques d'accélération) c'est ce fameux 'bruit' dont l’homogénéité dépend uniquement du générateur de nombres aléatoires.

Meilleur est ce générateur (dans le sens plus il tend à être réellement 'aléatoire'), meilleurs (plus homogène) sera ce bruit et donc plus agréable à l'oeil et plus facile à traiter avec des algorithmes de réduction de bruit.

Mersenne Twister est une sorte de standard par rapport au Monte Carlo car il est de grande qualité tout en étant assez rapide pour être utilisé avec un moteur de rendu.

[grabiller]

J'ai récupéré ton code sur le forum anglais pour ne plus avoir l'erreur avec setULong().
Maintenant j'ai une erreur à la 299
Random() : Max ne peut pas être négatif.
En effet la doc indique que la valeur maxi est 2147483647, et tu as #U32_MAX = 4294967295

C'est bizarre que tu n'aies pas cette erreur. Testé sous 5.20 beta 14

Je viens de tester avec la beta 14.

Alors effectivement, j'obtiens cette erreur mais en x86 seulement. En x64 le code marche correctement.

Je vais utiliser 2147483647 comme valeur max, c'est pas très grave c'est juste pour générer une valeur de seed aléatoire, mais cela ne remet pas en cause la génération des nombres eux-même.

Par contre en x86 les résultats ne sont pas bon dutout, je vais réintroduire les toULong/fromULong que j'avais retiré pensant que cela marchait sans (mais je n'avais pas testé sous x86..).

[grabiller]

Je viens de mettre à jour le code, cela devrait maintenant marcher aussi bien en x64 que en x86..

[djes]

Juste question d'algorithmie, dans la plupart des cas je ne vois pas l'intérêt d'utiliser ces générateurs de pseudos nombres aléatoires, sauf dans le contexte de simulations où il n'y a pas d'autres solutions (problèmes de maths, génétique...). La plupart du temps, il faut privilégier l'aspect rapidité, et l'aléatoire (simulé ici) ne sert à rien d'autre qu'à bouffer des cycles machines et à ralentir le programme !

A mon avis, et si éventuellement on tient à cet aspect, il est beaucoup plus intelligent et efficace de simuler une répartition aléatoire en ayant un tableau fixe de valeurs tirées au hasard, qu'on peut décliner à de multiples échelles au besoin.

Exemple avec quatre valeurs : [0, 1, 2, 3] -> [2, 1, 0, 3] (calculé une fois au début du programme)

On a besoin de quatre valeurs ? On ressort le tableau de base : [2, 1, 0, 3]

On a besoin d'une décimale ? µ = 1/4 + 1/4/2 = 0.25 + 0.125 -> [0.625, 0.175, 0.125, 0.875] ;(ça serait plus évident avec 10 valeurs, mais j'ai la flemme de taper ! )
[2, 2.625, 1.625, 0.625, 3.625,
1, 2.175, 1.175, 0.175, 3.175,
0, 2.125, 1.125, 0.125, 3.125,
3, 2.875, 1.875, 0.875, 3.875]

Vous voyez le principe. Un autre avantage que la rapidité est qu'on a toujours le même résultat, ce qui permet d'exempter son programme de bugs "aléatoires".

Edit : p'tite faute dans le tableau, sacrévindiou.

[PAPIPP]

@ grabiller

Pourquoi utiliser un générateur de nombre aléatoire alors que le générateur de PB est excellent

Voici un test du générateur de PB
Remarque :
La loi statistique est une loi uniforme c'est-à-dire que pour tout nombre compris entre un mini et un max la probabilité de chaque apparition est identique.
Voici un programme qui teste le générateur de PB. Il permet de vérifier que le générateur est fiable et qu’il génère bien une loi uniforme

Code : Tout sélectionner

;Loi uniforme d'un tirage par random ()
; la probabilité d'apparition d'un nombre entre A et B est de f(x)=1/(B-A)
; la moyenne est la Somme(x)/nb occurences=(A+A+1+A+2 ..... B-1+B)/nb(A+B)*n/2n=A+B/2
nb_tir=1000000  ; Nombre de tirage max
nmax=30; pour tirage statistiqie loi uniforme entre 0 et nmax (au mini 30 zones pour  vérifier la fiabilité du random de PB)
sx.q=0
sxx.q=0
Dim tab.l(nmax)
For i=0 To nb_tir
  X=Random(nmax)
  tab(x)+1
  sx+X
  sxx+(x*x)
Next
Moy.d=sx/nb_tir
Debug "moyenne="+StrD(Moy)
Debug "Somme_des carré="+Str(sxx)
smoyc.q=moy*moy*nb_tir
sigma2.D=(sxx-smoyc)/nb_tir
; Debug "sigma2="+StrD(sigma2)
sigma.d=Pow(sigma2,0.5)
Debug "Sigma="+StrD(sigma)
Debug "************ Table des résultats **************"
For i=0 To nmax
; Debug _n(i)+_n(tab(i))
Debug "i="+Str(i)+" tab(i)="+Str(tab(i))
Next
Debug " valeur théorique à trouver =" +Str((nb_tir+1)/(nmax+1))

et en utilisant le test de Wilbert

Code : Tout sélectionner

If OpenWindow(0, 0, 0, 512, 512, "Random generator", #PB_Window_SystemMenu | #PB_Window_ScreenCentered)
  If CreateImage(0, 512, 512) And StartDrawing(ImageOutput(0))
    
    For i = 0  To 131072
      R.l = Random(2147483647)
      ; 			R.l = Random(261121)
      X.l = R & 511
      Y.l = (R >> 16) & 511
      Plot(X, Y)
;       If (i%512)=0
;         Debug _n(R)+_n(X)+_n(Y)
;       EndIf
    Next
    
    StopDrawing() 
    ImageGadget(0, 0, 0, 200, 200, ImageID(0))
  EndIf
  
  Repeat
    Event = WaitWindowEvent()
  Until Event = #PB_Event_CloseWindow
EndIf

A+

[G-Rom]

Juste question d'algorithmie, dans la plupart des cas je ne vois pas l'intérêt d'utiliser ces générateurs de pseudos nombres aléatoires, sauf dans le contexte de simulations où il n'y a pas d'autres solutions (problèmes de maths, génétique...). La plupart du temps, il faut privilégier l'aspect rapidité, et l'aléatoire (simulé ici) ne sert à rien d'autre qu'à bouffer des cycles machines et à ralentir le programme !

A mon avis, et si éventuellement on tient à cet aspect, il est beaucoup plus intelligent et efficace de simuler une répartition aléatoire en ayant un tableau fixe de valeurs tirées au hasard, qu'on peut décliner à de multiples échelles au besoin.

Exemple avec quatre valeurs : [0, 1, 2, 3] -> [2, 1, 0, 3] (calculé une fois au début du programme)

On a besoin de quatre valeurs ? On ressort le tableau de base : [2, 1, 0, 3]

On a besoin d'une décimale ? µ = 1/4 + 1/4/2 = 0.25 + 0.125 -> [0.625, 0.175, 0.125, 0.875] ;(ça serait plus évident avec 10 valeurs, mais j'ai la flemme de taper ! )
[2, 2.625, 1.625, 0.625, 3.625,
1, 2.175, 1.175, 0.175, 3.175,
0, 2.125, 1.125, 0.125, 3.125,
3, 3.875, 1.875, 0.875, 3.875]

Vous voyez le principe. Un autre avantage que la rapidité est qu'on a toujours le même résultat, ce qui permet d'exempter son programme de bugs "aléatoires".

Pas bête du tout , ca me rappelle les optimisations des cos/sin à l'époque des démomakers, on vois bien les "vieux" ici
Mais je pense la solution tout à fait valable