Announcement

Collapse

Forum Guidelines

This forum is for finished source code that is working properly. If you have questions about this or any other source code, please post it in one of the Discussion Forums, not here.
See more
See less

Mersenne Twister PRNG for PBWIN70 & PBCC30

Collapse
X
 
  • Filter
  • Time
  • Show
Clear All
new posts

  • Mersenne Twister PRNG for PBWIN70 & PBCC30

    Code:
    #IF 0
    -- MT.BAS --
     [revised]
       
    The Mersenne Twister pseudo-random number generator has a 
    period of (2^19937-1).  More detailed information is available at:
     [url="http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html"]http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html[/url] 
          
    The PB code provided below incorporates several widely discussed 
    improvements to the original Mersenne Twister implementation supplied 
    by Makoto Matsumoto and Takuji Nishimura. Contributors to this public 
    discussion include Shawn Cokus, Richard Wagner, and Eric Landry. 
       
    George Marsaglia's DIEHARD utility is usually considered a standard 
    tool for examining the output of PRNGs. DIEHARD can be found at:
     [url="http://stat.fsu.edu/pub/diehard/"]http://stat.fsu.edu/pub/diehard/[/url] 
    The code below includes a routine that will create the 10MB file 
    required for running the entire suite of Diehard tests. Specifying 
    the output file name (by defining the $DHFILE equate) causes the 
    file to be created.
       
    Note that this PRNG has an enormous period, but it possesses no 
    cryptographic security. 
       
    The following code compiles with either PBWIN70+ or PBCC30+.
       
    This implementation replaces the original, which was posted in 2002. 
    Advantages of this new version include enhanced efficiency, no 
    reliance on global data, and the ability to run multiple generators 
    simultaneously.
       
    -- Greg Turgeon  1/2005
    gturgeon at verizon dot net
    #ENDIF 
       
    #COMPILE EXE
    #DIM ALL
    #REGISTER NONE
    '============
       
    '-- Utility macros --
       
    #IF %def(%PB_WIN32)
    '--------------------
    MACRO eol=$CR
    '--------------------
    MACRO ShowText(T)
       msgbox T
    END MACRO
    #ELSE
    '--------------------
    MACRO eol=$CRLF
    '--------------------
    MACRO ShowText(T)
       stdout T
    END MACRO
    #ENDIF
       
    '--------------------
    MACRO zbs(x)=string$(x,0)
       
    '--------------------
    MACRO EnterCC
    #IF %DEF(%PB_CC32)
    LOCAL launched&
    if (cursory = 1) and (cursorx = 1) then launched& = -1
    #ENDIF
    END MACRO
       
    '--------------------
    MACRO ExitCC
    #IF %DEF(%PB_CC32)
    if launched& then
       input flush
       print "Press any key to end"
       waitkey$
    end if
    #ENDIF
    END MACRO
       
    '--------------------
    MACRO FUNCTION shiftl(x,shiftval)
    MACROTEMP retval
    LOCAL retval AS DWORD
    retval = x
    !  shl   retval, shiftval
    END MACRO = retval
       
    '--------------------
    MACRO FUNCTION shiftr(x,shiftval)
    MACROTEMP retval
    LOCAL retval AS DWORD
    retval = x
    !  shr   retval, shiftval
    END MACRO = retval
       
       
       
    '-- MT19937 equates & macros --
    %N        = 624
    %M        = 397
       
    %MATRIX_A = &h9908b0df???
    %UMASK    = &h80000000???
    %LMASK    = &h7fffffff???
       
    '--------------------
    MACRO MixBits(v,w)=(((v) AND %UMASK) OR ((w) AND %LMASK))
       
    '--------------------
    MACRO FUNCTION Twist(u,v,w)
    MACROTEMP x, y
    LOCAL x AS DWORD, y AS DWORD
    x = MixBits(v,w)
    !  shr   x, 1
    y = ((-(w AND 1)) AND %MATRIX_A)
    END MACRO = (u XOR x XOR y)
       
       
    TYPE MT_CONTEXT
       pState   AS DWORD PTR
       State    AS STRING * (4*(%N+1))
       Seed     AS DWORD
       Idx      AS LONG
    END TYPE
       
    DECLARE FUNCTION InitMT&(Ctx AS MT_CONTEXT)
    DECLARE FUNCTION RandomMT???(Ctx AS MT_CONTEXT)
       
    '-- Unrem and supply file name for Diehard data file
    '$DHFILE = "f:\diehard\mt.$$$"
       
    #IF %def($DHFILE)
    DECLARE FUNCTION MakeDiehardFile&()
    #ENDIF
       
    '====================
    FUNCTION PBMain&()
    REGISTER i&, j&
    LOCAL t$, ctx AS MT_CONTEXT
    EnterCC
       
    ctx.Seed = 5489
    InitMT ctx
       
    t = ""
    for i = 1 to 5
       for j = 1 to 5
          t = t + using$("############", RandomMT(ctx)) + "  "
       next j
       t = t + eol
    next i
    ShowText(t)
       
    #IF %def($DHFILE)
    MakeDiehardFile
    #ENDIF
       
    ExitCC
    END FUNCTION
       
       
    '====================
    FUNCTION InitMT&(Ctx AS MT_CONTEXT)
    LOCAL i???, x???, y???, pd AS DWORD PTR
       
    '-- Default seed if none supplied
    if Ctx.Seed = 0 then Ctx.Seed = 5489
    '-- Zero initial state & save pointer to it
    Ctx.State  = zbs(sizeof(Ctx.State))
    Ctx.pState = varptr(Ctx.State)
       
    pd  = Ctx.pState
    @pd = Ctx.Seed
    for i = 1 to %N-1
       x = @pd
       y = x XOR shiftr(x,30)
       'y = (y * 1812433253) + i
       'Be sure to handle PB DWORD 
       'multiply this way
       !  mov   edx, y
       !  mov   eax, 1812433253
       !  mul   edx
       !  mov   x, eax
       incr pd
       @pd = x + i
    next i
       
    '-- Force reload on first call to RandomMT()
    Ctx.Idx = %N+1
    END FUNCTION
       
       
    '====================
    FUNCTION RandomMT???(Ctx AS MT_CONTEXT)
    REGISTER i&
    LOCAL x???, pd AS DWORD PTR
       
    if Ctx.Idx >= %N then 'reload state array
       'if RandomMT() called w/o seeding
       if [email protected][0] = 0 then
          Ctx.Seed = 5489 : InitMT(Ctx)
       end if
       
       'local copy for speed
       pd = Ctx.pState
       
       for i = 0 to (%N-%M)-1
          @pd = Twist(@pd[%M],@pd[0],@pd[1])
          incr pd
       next i
       
       for i = (%N-%M) to (%N-2)
          @pd = Twist(@pd[%M-%N],@pd[0],@pd[1])
          incr pd
       next i
       
       @pd = Twist(@pd[%M-%N],@pd[0],[email protected][0])
       Ctx.Idx = 0
    end if
       
    x = [email protected][Ctx.Idx] : incr Ctx.Idx
       
    x = x XOR  shiftr(x,11)
    x = x XOR (shiftl(x, 7) AND &h9D2C5680)
    x = x XOR (shiftl(x,15) AND &hEFC60000)
    function = x XOR shiftr(x,18)
    END FUNCTION
       
    #IF %def($DHFILE)
    '-- NOTE: Routine performs no file access error checking.
    '============
    FUNCTION MakeDiehardFile&()
    REGISTER i&, j&
    LOCAL  t$, buffer$, outfile&
    LOCAL ctx AS MT_CONTEXT, pd AS DWORD PTR
       
    '-- Truncate any existing $DHFILE
    outfile = freefile
    open $DHFILE for binary as outfile base = 0
    seteof outfile
    close outfile
       
    ctx.Seed = 5489
    InitMT ctx
       
    for i = 1 to 10
       buffer = zbs(1000000)
       pd = strptr(buffer)
       j = len(buffer)\sizeof(@pd)
       do while j
          @pd = RandomMT(ctx) : incr pd
          decr j
       loop
       open $DHFILE for binary as outfile base = 0
       seek outfile, lof(outfile)
       put$ outfile, buffer
       if i = 10 then
          t = "Size of "+ $DHFILE + ":"+ str$(lof(outfile))
          ShowText(t)
       end if
       close outfile
    next i
    END FUNCTION
    #ENDIF
    '-- end MT.BAS --

    [This message has been edited by Greg Turgeon (edited February 08, 2005).]

  • #2
    I was running some random #s using the Mersenne Twister code here
    but there seems to be a problem. Tho the numbers vary, they don't
    pass any of the random # tests (like Diehard) and some numbers
    repeat often.

    Has anyone else tried this code with success? Maybe I'm just
    missing something simple. I increased the limits of the 2 loops
    and wrote the resulting DWORDS to a file as shown in the code
    below. The DWORDS produced are not random as generated. Note: I
    did convert them to binary (four-byte DWORDS) and tested them with
    ENT and Diehard, but it's clear just looking at them that they're
    not random.

    The only code I changed is shown below in the PBMAIN function:

    Code:
    FUNCTION PBMAIN&()
    REGISTER i&, j&
    LOCAL t$
    EnterCC
    
    REDIM STATE(%N)
    STATE_REMAINING = -1
    CALL SeedMT(4357)
    OPEN "c:\mersenne2.txt" FOR OUTPUT AS #1  'I create a file
    
    FOR i& = 1 TO 500            'increase loop size
       FOR j& = 1 TO 50          '     "     "    "
         PRINT #1,  RandomMT     'print the random #s to it          
       NEXT j&
    '   t$ = t$ + $CRLF
    NEXT i&
    CLOSE
    
    ShowText(t$)
    
    '============
    ExitPBMain:
    ExitCC
    END FUNCTION 'Main
    [This message has been edited by John Gleason (edited January 26, 2005).]

    Comment


    • #3
      John--

      As I looked at the original code (posted more than two years ago), I realized how dated it had become. I've now replaced it with a version I did for a recent project. Please dump the previous version and replace it with the new code. (I also suspect that the problem you were having resulted from a formatting change that occurred while posting the previous code--a change in placement of parentheses, which affected the order of operations.)


      [This message has been edited by Greg Turgeon (edited January 29, 2005).]

      Comment


      • #4
        Code:
        #IF 0
        -- MT.BAS --
         [revised]
        
        The Mersenne Twister pseudo-RANDOM number generator has a
        period OF (2^19937-1).  More detailed information is available AT:
         [url="http://www.math.sci.hiroshima-u.ac.jp/~m-MAT/MT/emt.html"]http://www.math.sci.hiroshima-u.ac.jp/~m-MAT/MT/emt.html[/url]
        
        The PB code provided below incorporates several widely discussed
        improvements TO the original Mersenne Twister implementation supplied
        by Makoto Matsumoto AND Takuji Nishimura. Contributors TO this public
        discussion include Shawn Cokus, Richard Wagner, AND Eric Landry.
        
        George Marsaglia's DIEHARD utility is usually considered a standard
        tool FOR examining the OUTPUT OF PRNGs. DIEHARD can be found AT:
         [url="http://stat.fsu.edu/pub/diehard/"]http://stat.fsu.edu/pub/diehard/[/url]
        The code below includes a routine that will CREATE the 10MB file
        required FOR running the entire suite OF Diehard tests. Specifying
        the OUTPUT file NAME (by defining the $DHFILE equate) causes the
        file TO be created.
        
        Note that this PRNG has an enormous period, but it possesses no
        cryptographic security.
        
        The following code compiles WITH either PBWIN70+ OR PBCC30+.
        
        This implementation replaces the original, which was posted IN 2002.
        Advantages OF this NEW version include enhanced efficiency, no
        reliance ON GLOBAL DATA, and the ability to run multiple generators
        simultaneously.
        
        -- Greg Turgeon  1/2005
        gturgeon AT verizon dot net
        ===========================
        This version OF Greg's code includes a macro to do the inner loop code with asm
        FOR a speed boost. I tested it only WITH PBWIN70+ but it probably should run
        IN the PBCC too, AS I tried NOT TO change ANY OF the CC code.
        -- John Gleason  2/2005
        
        #ENDIF
        
        #COMPILE EXE
        #DIM ALL
        #REGISTER NONE
        '============
        '-- MT19937 equates & macros --
        %N        = (624-1)
        %Nm1      = %N - 1
        %M        = 397
        %NmM      = %N - %M
        %NmMp1    = %NmM + 1
        %MmN      = %M - %N
        %MATRIX_A = &h9908b0df???
        %UMASK    = &h80000000???
        %LMASK    = &h7fffffff???
        %MTSEED   = 399495995      'this is the seed dword integer--1 to 4,294,967,295
        
        TYPE MT_CONTEXT
           pState   AS DWORD PTR
           STATE    AS STRING * (4*(%N+1))
           Seed     AS DWORD
           Idx      AS LONG
        END TYPE
        
        
        '-- Utility macros --
        MACRO mRandomMT
        
        IF Ctx.Idx > %N THEN 'reload state array
            
           pd = Ctx.pState
           FOR i = 0 TO %NmM
               !mov eax, pd
               !mov ecx, %M          ;index of twist1
               !mov edx, [eax+ecx*4] ;x4 = twist1 in edx
               !mov ecx, [eax]       ;put twist2 in ecx here to avoid dependency--for speed
               !mov twist1, edx
               !mov edx, [eax+4]     ;twist3 is second dword (index 1)
               !mov twist3, edx
                     'now we have all the parameters needed for Twist
                     '1st we MixBits of twist2 and twist3
               !and ecx, %UMASK      ;twist2 still is in ecx, so AND it
               !and edx, %LMASK      ;twist3 still is in edx, so AND it
               !or  ecx, edx
               !mov eax, twist3      ;look ahead for no dependency
               !shr ecx, 1           ;ecx is now x properly shifted
               !and eax, 1           ;twist3 AND 1
               !neg eax              ;-(twist3 AND 1)
               !and eax, %MATRIX_A   ;eax is now y
               !mov edx, twist1
               !xor eax, ecx         ;1st xor
               !xor eax, edx         ;2nd (and last) xor
                     'eax now has final result of the Twist macro
               !mov ecx, pd
               !mov [ecx], eax       ;so move it to the proper context string position
               !add ecx, 4           ;move up on context string 4--dword length
               !mov pd,  ecx         ;save it for next loop
                     'and we're done! Looks like we did it in 24 ticks--not too shabby
           NEXT i
        
           FOR i =  %NmMp1 TO %Nm1
               !mov eax, pd
               !mov ecx, %MmN        ;index of twist1
               !mov edx, [eax+ecx*4] ;x4 = twist1 in edx
               !mov ecx, [eax]       ;put twist2 in ecx here to avoid dependency--for speed
               !mov twist1, edx
               !mov edx, [eax+4]     ;twist3 is second dword (index 1)
               !mov twist3, edx
                     'now we have all the parameters needed for Twist
                     '1st MixBits of twist2 and twist3
               !and ecx, %UMASK      ;twist2 still is in ecx, so AND it
               !and edx, %LMASK      ;twist3 still is in edx, so AND it
               !or  ecx, edx
               !mov eax, twist3      ;look ahead for no dependency
               !shr ecx, 1           ;ecx is now x properly shifted
               !and eax, 1           ;twist3 AND 1
               !neg eax              ;-(twist3 AND 1)
               !and eax, %MATRIX_A   ;eax is now y
               !mov edx, twist1
               !xor eax, ecx         ;1st xor
               !xor eax, edx         ;2nd (and last) xor
                     'eax now has final result of the Twist macro
               !mov ecx, pd
               !mov [ecx], eax       ;so move it to the proper context string position
               !add ecx, 4           ;move up on context string 4--dword length
               !mov pd,  ecx         ;save it for next loop
           NEXT i
        
           twist3 = [email protected][0]
               !mov eax, pd
               !mov ecx, %MmN        ;index of twist1
               !mov edx, [eax+ecx*4] ;x4 = twist1 in edx
               !mov ecx, [eax]       ;put twist2 in ecx here to avoid dependency--for speed
               !mov twist1, edx
               !mov edx, twist3      ;twist3 is defined above as [email protected][0]
                     'now we have all the parameters needed for Twist
                     '1st MixBits of twist2 and twist3
               !and ecx, %UMASK      ;twist2 still is in ecx, so AND it
               !and edx, %LMASK      ;twist3 still is in edx, so AND it
               !or  ecx, edx
               !mov eax, twist3      ;look ahead for no dependency
               !shr ecx, 1           ;ecx is now x properly shifted
               !and eax, 1           ;twist3 AND 1
               !neg eax              ;-(twist3 AND 1)
               !and eax, %MATRIX_A   ;eax is now y
               !mov edx, twist1
               !xor eax, ecx         ;1st xor
               !xor eax, edx         ;2nd (and last) xor
                     'eax now has final result of the Twist macro
               !mov ecx, pd
               !mov [ecx], eax       ;so move it to the proper context string position
           Ctx.Idx = 0
        END IF
        x = [email protected][Ctx.Idx]
        
        INCR Ctx.Idx
        
               !mov eax, x
               !mov ecx, x
               !shr eax, 11       ;x = x XOR  shiftr(x,11)
               !xor eax, ecx
        '------------------------
               !mov ecx, eax
               !shl eax, 7
               !and eax, &h9D2C5680    ;x = x XOR (shiftl(x, 7) AND &h9D2C5680)
               !xor eax, ecx
        '------------------------
               !mov ecx, eax
               !shl eax, 15
               !and eax, &hEFC60000   ;x = x XOR (shiftl(x,15) AND &hEFC60000)
               !xor eax, ecx
        '------------------------
               !mov ecx, eax
               !shr eax, 18
               !xor eax, ecx           ;x = x XOR shiftr(x,18)
        '      !mov x, eax
        
        END MACRO
        
        #IF %DEF(%PB_WIN32)
        '--------------------
        MACRO eol=$CR
        '--------------------
        MACRO ShowText(T)
           MSGBOX T
        END MACRO
        #ELSE
        '--------------------
        MACRO eol=$CRLF
        '--------------------
        MACRO ShowText(T)
           stdout T
        END MACRO
        #ENDIF
        
        '--------------------
        MACRO zbs(x)=STRING$(x,0)
        
        '--------------------
        MACRO EnterCC
        #IF %DEF(%PB_CC32)
        LOCAL launched&
        IF (cursory = 1) AND (cursorx = 1) THEN launched& = -1
        #ENDIF
        END MACRO
        
        '--------------------
        MACRO ExitCC
        #IF %DEF(%PB_CC32)
        IF launched& THEN
           INPUT FLUSH
           PRINT "Press any key to end"
           waitkey$
        END IF
        #ENDIF
        END MACRO
        
        '--------------------
        MACRO FUNCTION shiftl(x,shiftval)
        MACROTEMP retval
        LOCAL retval AS DWORD
        retval = x
        !  shl   retval, shiftval
        END MACRO = retval
        
        '--------------------
        MACRO FUNCTION shiftr(x,shiftval)
        MACROTEMP retval
        LOCAL retval AS DWORD
        retval = x
        !  shr   retval, shiftval
        END MACRO = retval
        
        '--------------------
        MACRO MixBits(v,w)=(((v) AND %UMASK) OR ((w) AND %LMASK))
        
        '--------------------
        MACRO FUNCTION Twist(u,v,w)
        MACROTEMP x, y
        LOCAL x AS DWORD, y AS DWORD
        x = MixBits(v,w)
        !  shr   x, 1
        y = ((-(w AND 1)) AND %MATRIX_A)
        END MACRO = (u XOR x XOR y)
        
        
        DECLARE FUNCTION InitMT&(Ctx AS MT_CONTEXT)
        DECLARE FUNCTION RandomMT???(Ctx AS MT_CONTEXT)
        
        '-- Unrem and supply file name for Diehard data file
        $DHFILE = "c:\mersenneTwist1x.dat"
        
        #IF %DEF($DHFILE)
        DECLARE FUNCTION MakeDiehardFile&()
        #ENDIF
        
        '====================
        FUNCTION PBMAIN&()
        REGISTER i&, j&
        LOCAL t$, ctx AS MT_CONTEXT
        'EnterCC
        
        ctx.Seed = %MTSEED
        InitMT ctx
        
        t = ""
        FOR i = 1 TO 5
           FOR j = 1 TO 5
              t = t + USING$("############", RandomMT(ctx)) + "  "
           NEXT j
           t = t + eol
        NEXT i
        ShowText(t)
        
        #IF %DEF($DHFILE)
        MakeDiehardFile
        #ENDIF
        
        ExitCC
        END FUNCTION
        
        
        '====================
        FUNCTION InitMT&(Ctx AS MT_CONTEXT)
        LOCAL i???, x???, y???, pd AS DWORD PTR
        
        '-- Default seed if none supplied
        IF Ctx.Seed = 0 THEN Ctx.Seed = %MTSEED
        '-- Zero initial state & save pointer to it
        Ctx.State  = zbs(SIZEOF(Ctx.State))
        Ctx.pState = VARPTR(Ctx.State)
        
        pd  = Ctx.pState
        @pd = Ctx.Seed
        FOR i = 1 TO %N
           x = @pd
           y = x XOR shiftr(x,30)
           'y = (y * 1812433253) + i
           'Be sure to handle PB DWORD
           'multiply this way
           !  mov   edx, y
           !  mov   eax, 1812433253
           !  mul   edx
           !  mov   x, eax
           INCR pd
           @pd = x + i
        NEXT i
        
        '-- Force reload on first call to RandomMT()
        Ctx.Idx = %N+1
        END FUNCTION
        
        
        '====================
        FUNCTION RandomMT???(Ctx AS MT_CONTEXT)
        REGISTER i AS LONG
        LOCAL x???, pd AS DWORD PTR
        
        IF Ctx.Idx > %N THEN 'reload state array
        
           'local copy for speed
           pd = Ctx.pState
        
           FOR i = 0 TO %NmM
              @pd = Twist(@pd[%M],@pd[0],@pd[1])
              INCR pd
           NEXT i
        
           FOR i =  %NmMp1 TO %Nm1
              @pd = Twist(@pd[%MmN],@pd[0],@pd[1])
              INCR pd
           NEXT i
        
           @pd = Twist(@pd[%MmN],@pd[0],[email protected][0])
           Ctx.Idx = 0
        END IF
        
        x = [email protected][Ctx.Idx]
        INCR Ctx.Idx
        
        x = x XOR  shiftr(x,11)
        x = x XOR (shiftl(x, 7) AND &h9D2C5680)
        x = x XOR (shiftl(x,15) AND &hEFC60000)
        FUNCTION = x XOR shiftr(x,18)
        
        END FUNCTION
        
        #IF %DEF($DHFILE)
        '-- NOTE: Routine performs no file access error checking.
        '============
        FUNCTION MakeDiehardFile&()
        'REGISTER i&, j&
        LOCAL ii AS LONG, j AS LONG
        LOCAL i AS LONG, i2 AS LONG
        LOCAL  t$, buffer$, outfile&
        LOCAL ctx AS MT_CONTEXT, pd AS DWORD PTR
        LOCAL x AS DWORD, pd2 AS DWORD PTR, basePd2 AS LONG
        LOCAL twist1 AS DWORD, twist2 AS DWORD, twist3 AS DWORD
        LOCAL aTime1 AS DOUBLE, aTime2 AS DOUBLE, macroTimeTot AS DOUBLE, functionTimeTot AS DOUBLE
        LOCAL pupEbx AS LONG 'push/pop register substitute to save 3 ticks / loop
        '-- Truncate any existing $DHFILE
        outfile = FREEFILE
        OPEN $DHFILE FOR BINARY AS outfile BASE = 0
        SETEOF outfile
        CLOSE outfile
        
        ctx.Seed = %MTSEED
        InitMT ctx
        
           buffer = zbs(1000000)
           pd2 = STRPTR(buffer)
           basePd2 = pd2
           j = LEN(buffer)\SIZEOF(@pd2)
           OPEN $DHFILE FOR BINARY AS outfile BASE = 0
        
           aTime1 = TIMER
        FOR i2 = 1 TO 40                       'this loop uses the new asm macro
            pd2 = basePd2
            FOR ii = 1 TO j
               mRandomMT
               !mov ecx,   pd2  ;buffer pointer
               !mov [ecx], eax  ;write the randomMT dword to buffer
               !add ecx, 4      ;incr the buffer pointer 4
               !mov pd2, ecx    ;and save it for the next loop
            NEXT
            '   PUT outfile, , buffer          'uncomment this line to write to disk
        NEXT
           aTime2 = TIMER
           macroTimeTot = aTime2 - aTime1
        
           aTime1 = TIMER
        FOR i2 = 1 TO 40                       'this loop uses the original function
           pd = basePd2
           j = LEN(buffer)\SIZEOF(@pd)
           DO WHILE j
              @pd = RandomMT(ctx) : INCR pd
              DECR j
           LOOP
             '  PUT outfile, , buffer          'uncomment this line to write to disk
        NEXT
           aTime2 = TIMER
           functionTimeTot = aTime2 - aTime1
        
              t = "Size of "+ $DHFILE + ":"+ STR$(LOF(outfile)) & $CRLF & _
                        "Macro sec: " & FORMAT$(macroTimeTot, "0.00") & ";  Function sec: " & FORMAT$(functionTimeTot, "0.00")
              ShowText(t)
           CLOSE outfile
        END FUNCTION
        #ENDIF
        '-- end MT.BAS --

        Comment


        • #5
          Code:
          '-- Well, if you want speed, do the tempering in advance and buffer the result.
          '-- And add MMX, if available (as it almost surely will be).
          #COMPILE EXE
          #DIM ALL
          #REGISTER NONE
          '============
             
          '-- Utility macros --
             
          #IF %def(%PB_WIN32)
          '--------------------
          MACRO eol=$CR
          '--------------------
          MACRO ShowText(T)
             msgbox T
          END MACRO
          #ELSE
          '--------------------
          MACRO eol=$CRLF
          '--------------------
          MACRO ShowText(T)
             stdout T
          END MACRO
          #ENDIF
             
          '--------------------
          MACRO zbs(x)=string$(x,0)
             
          '--------------------
          MACRO EnterCC
          #IF %DEF(%PB_CC32)
          LOCAL launched&
          if (cursory = 1) and (cursorx = 1) then launched& = -1
          #ENDIF
          END MACRO
             
          '--------------------
          MACRO ExitCC
          #IF %DEF(%PB_CC32)
          if launched& then
             input flush
             print "Press any key to end"
             waitkey$
          end if
          #ENDIF
          END MACRO
             
          '--------------------
          MACRO FUNCTION shiftl(x,shiftval)
          MACROTEMP retval
          LOCAL retval AS DWORD
          retval = x
          !  shl   retval, shiftval
          END MACRO = retval
             
          '--------------------
          MACRO FUNCTION shiftr(x,shiftval)
          MACROTEMP retval
          LOCAL retval AS DWORD
          retval = x
          !  shr   retval, shiftval
          END MACRO = retval
             
             
          '-- MT19937 equates & macros --
          %N        = 624
          %M        = 397
             
          %MATRIX_A = &h9908b0df???
          %UMASK    = &h80000000???
          %LMASK    = &h7fffffff???
             
          '--------------------
          MACRO MixBits(v,w)=(((v) AND %UMASK) OR ((w) AND %LMASK))
             
          '--------------------
          MACRO FUNCTION Twist(u,v,w)
          MACROTEMP x, y
          LOCAL x AS DWORD, y AS DWORD
          x = MixBits(v,w)
          !  shr   x, 1
          y = ((-(w AND 1)) AND %MATRIX_A)
          END MACRO = (u XOR x XOR y)
             
             
          %BUFFERSIZE  = (4*(%N))
          TYPE MT_CONTEXT
             pState   AS DWORD PTR
             Idx      AS LONG
             Seed     AS DWORD
             State    AS STRING * %BUFFERSIZE
             Buffer   AS STRING * %BUFFERSIZE
             pBuffer  AS DWORD PTR
             HasMMX   AS LONG
          END TYPE
          %BUFFEROFFSET  = %BUFFERSIZE
             
          DECLARE FUNCTION InitMT&(Ctx AS MT_CONTEXT)
          DECLARE FUNCTION RandomMT???(Ctx AS MT_CONTEXT)
          DECLARE FUNCTION HasMMX&()
             
          '-- Unrem and supply file name for Diehard data file
          $DHFILE = "f:\diehard\mt.$$$"
             
          #IF %def($DHFILE)
          DECLARE FUNCTION MakeDiehardFile&()
          #ENDIF
             
          '====================
          FUNCTION PBMain&()
          REGISTER i&, j&
          LOCAL t$, ctx1 AS MT_CONTEXT
          EnterCC
             
          ctx1.Seed = 5489
          InitMT ctx1
             
          t = ""
          for i = 1 to 5
             for j = 1 to 5
                t = t + using$("############", RandomMT(ctx1))
             next j
             t = t + eol
          next i
          ShowText(t)
             
          #IF %def($DHFILE)
          MakeDiehardFile
          #ENDIF
             
          ExitPBMain:
          ExitCC
          END FUNCTION
             
             
          '====================
          FUNCTION InitMT&(Ctx AS MT_CONTEXT)
          LOCAL i???, x???, y???, pd AS DWORD PTR
             
          '-- Default seed if none supplied
          if Ctx.Seed = 0 then Ctx.Seed = 5489
          '-- Zero initial state & buffer save pointers to them
          Ctx.State  = zbs(sizeof(Ctx.State))
          Ctx.pState = varptr(Ctx.State)
          Ctx.Buffer  = zbs(sizeof(Ctx.Buffer))
          Ctx.pBuffer = varptr(Ctx.Buffer)
          Ctx.HasMMX  = HasMMX&
             
          pd  = Ctx.pState
          @pd = Ctx.Seed
          for i = 1 to %N-1
             x = @pd
             y = x XOR shiftr(x,30)
             'y = (y * 1812433253) + i
             'Be sure to handle PB DWORD
             'multiply this way
             !  mov   edx, y
             !  mov   eax, 1812433253
             !  mul   edx
             !  mov   x, eax
             incr pd
             @pd = x + i
          next i
             
          '-- Force reload on first call to RandomMT()
          Ctx.Idx = %N+1
          END FUNCTION
             
          '--------------------
          MACRO TemperEAX
          'x = x XOR  shiftr(x,11)
          !  mov         edx, eax
          !  shr         edx, 11
          !  xor         eax, edx
          'x = x XOR (shiftl(x,7) AND &h9D2C5680)
          !  mov         edx, eax
          !  shl         edx, 7
          !  and         edx, &h9D2C5680
          !  xor         eax, edx
          'x = x XOR (shiftl(x,15) AND &hEFC60000)
          !  mov         edx, eax
          !  shl         edx, 15
          !  and         edx, &hEFC60000
          !  xor         eax, edx
          'x = x XOR shiftr(x,18)
          !  mov         edx, eax
          !  shr         edx, 18
          !  xor         eax, edx
          END MACRO
             
          '--------------------
          MACRO TemperMMXregs
          '-- Uses mm3 & mm7 as scratch regs
          ''x = x XOR  shiftr(x,11)
          !  movq     mm3, mm0
          !  movq     mm7, mm4
          !  psrld    mm3, 11
          !  psrld    mm7, 11
          !  pxor     mm0, mm3
          !  pxor     mm4, mm7
          'x = x XOR (shiftl(x,7) AND &h9D2C5680)
          !  movq     mm3, mm0
          !  movq     mm7, mm4
          !  pslld    mm3, 7
          !  pslld    mm7, 7
          !  pand     mm3, QWORD PTR Tempering1
          !  pand     mm7, QWORD PTR Tempering1
          !  pxor     mm0, mm3
          !  pxor     mm4, mm7
          'x = x XOR (shiftl(x,15) AND &hEFC60000)
          !  movq     mm3, mm0
          !  movq     mm7, mm4
          !  pslld    mm3, 15
          !  pslld    mm7, 15
          !  pand     mm3, QWORD PTR Tempering2
          !  pand     mm7, QWORD PTR Tempering2
          !  pxor     mm0, mm3
          !  pxor     mm4, mm7
          'x = x XOR shiftr(x,18)
          !  movq     mm3, mm0
          !  movq     mm7, mm4
          !  psrld    mm3, 18
          !  psrld    mm7, 18
          !  pxor     mm0, mm3
          !  pxor     mm4, mm7
          END MACRO
             
             
          '====================
          FUNCTION RandomMT???(Ctx AS MT_CONTEXT)
             
          if Ctx.Idx >= %N then 'reload state array
             'if RandomMT2() called w/o seeding
             if [email protected][0] = 0 then
                Ctx.Seed = 5489 : InitMT(Ctx)
             end if
             
             if Ctx.HasMMX then
                !  push        ebx
                !  push        esi
                !  push        edi
             
                '  for i = 0 to (%N-%M)
                '     @pd = Twist(@pd[%M],@pd[0],@pd[1]) : incr pd
                '  next i
                !  mov         esi, ctx
                !  push        esi         ;save ctx[0] for last operation
                !  mov         esi, [esi]
                !  push        esi         ;save [ctx[0]]for last operation
                !  mov         edi, (4*%M)
                !  mov         ecx, 56                 ;0-223
                LOOP1:
                !  movq        mm1, [esi+4]
                !  movq        mm0, [esi]
                !  movq        mm5, [esi+12]
                !  movq        mm4, [esi+8]
             
                !  movq        mm2, mm1    ;copy
                !  movq        mm6, mm5    ;copy
             
                !  pand        mm0, DoubleUMask
                !  pand        mm4, DoubleUMask
                !  pand        mm1, DoubleLMask
                !  pand        mm5, DoubleLMask
                !  por         mm0, mm1
                !  por         mm4, mm5
                !  psrld       mm0, 1
                !  psrld       mm4, 1
             
                !  pand        mm2, Double1
                !  pand        mm6, Double1
                !  pcmpeqd     mm2, Double1
                !  pcmpeqd     mm6, Double1
                !  pand        mm2, DoubleMatrixA
                !  pand        mm6, DoubleMatrixA
             
                !  pxor        mm0, mm2
                !  pxor        mm4, mm6
                !  pxor        mm0, [esi+edi]
                !  pxor        mm4, [esi+edi+8]
                !  movq        [esi], mm0
                !  movq        [esi+8], mm4
             
                   TemperMMXregs
             
                !  movq        [esi+%BUFFEROFFSET], mm0
                !  add         esi, 8
                !  movq        [esi+%BUFFEROFFSET], mm4
                !  add         esi, 8
                !  dec         ecx
                !  jnz         LOOP1
             
                '-- three remaining
                !  mov         ecx, 3         ;224-226
                Loop2:
                !   mov        eax, [esi]
                !   mov        ebx, [esi+4]
                !   mov        edx, ebx
             
                !  and         eax, %UMASK
                !  and         ebx, %LMASK
                !  or          eax, ebx
                !  shr         eax, 1
             
                !  and         edx, 1
                !  neg         edx
                !  sbb         edx, edx
                !  and         edx, %MATRIX_A
                !  xor         eax, edx
                !  xor         eax, [esi+edi]
                !  mov         [esi], eax
                   TemperEAX
                !  mov         [esi+%BUFFEROFFSET], eax
                !  add         esi, 4
                !  dec         ecx
                !  jnz         Loop2
             
                '  for i = (%N-%M)+1 to (%N-1)
                '     @pd = Twist(@pd[%M-%N],@pd[0],@pd[1]) : incr pd
                '  next i
                !  mov         edi, 4*(%M-%N)
                !  mov         ecx, 99           ;227-622
                LOOP3:
                !  movq        mm1, [esi+4]
                !  movq        mm0, [esi]
                !  movq        mm5, [esi+12]
                !  movq        mm4, [esi+8]
             
                !  movq        mm2, mm1    ;copy
                !  movq        mm6, mm5    ;copy
             
                !  pand        mm0, DoubleUMask
                !  pand        mm4, DoubleUMask
                !  pand        mm1, DoubleLMask
                !  pand        mm5, DoubleLMask
                !  por         mm0, mm1
                !  por         mm4, mm5
                !  psrld       mm0, 1
                !  psrld       mm4, 1
             
                !  pand        mm2, Double1
                !  pand        mm6, Double1
                !  pcmpeqd     mm2, Double1
                !  pcmpeqd     mm6, Double1
                !  pand        mm2, DoubleMatrixA
                !  pand        mm6, DoubleMatrixA
             
                !  pxor        mm0, mm2
                !  pxor        mm4, mm6
                !  pxor        mm0, [esi+edi]
                !  pxor        mm4, [esi+edi+8]
                !  movq        [esi], mm0
                !  movq        [esi+8], mm4
             
                   TemperMMXregs
             
                !  movq        [esi+%BUFFEROFFSET], mm0
                !  add         esi, 8
                !  movq        [esi+%BUFFEROFFSET], mm4
                !  add         esi, 8
             
                !  dec         ecx
                !  jnz         LOOP3
             
             
                '  @pd = Twist(@pd[%M-%N],@pd[0],[email protected][0])
                !  mov         eax, [esi]  ;last value of 0-623
                !  pop         ebx         ;edx = [ctx[0]]
                !  mov         ebx, [ebx]
                !  mov         edx, ebx    ; save
             
                !  and         eax, %UMASK
                !  and         ebx, %LMASK
                !  or          eax, ebx
                !  shr         eax, 1
             
                !  and         edx, 1
                !  neg         edx
                !  sbb         edx, edx
                !  and         edx, %MATRIX_A
             
                !  xor         eax, edx
                !  xor         eax, [esi+edi]
                !  mov         [esi], eax
                   TemperEAX
                !  mov         [esi+%BUFFEROFFSET], eax
             
                '  Ctx.Idx = 0
                !  pop         esi      ;esi = ctx[0]
                !  xor         ebx, ebx
                !  mov         [esi][4], ebx
                !  pop         edi
                !  pop         esi
                !  pop         ebx
                !  emms
             
             else ' no MMX
             
                LOCAL i&, x???, pds AS DWORD PTR, pdb AS DWORD PTR
                pds = Ctx.pState : pdb = Ctx.pBuffer
             
                for i = 0 to (%N-%M)-1
                   @pds = Twist(@pds[%M],@pds[0],@pds[1])
                   x = @pds XOR shiftr(@pds,11)
                   x = x XOR (shiftl(x, 7) AND &h9D2C5680)
                   x = x XOR (shiftl(x,15) AND &hEFC60000)
                   @pdb = x XOR shiftr(x,18)
                   incr pds : incr pdb
                next i
                for i = (%N-%M) to (%N-2)
                   @pds = Twist(@pds[%M-%N],@pds[0],@pds[1])
                   x = @pds XOR shiftr(@pds,11)
                   x = x XOR (shiftl(x, 7) AND &h9D2C5680)
                   x = x XOR (shiftl(x,15) AND &hEFC60000)
                   @pdb = x XOR shiftr(x,18)
                   incr pds : incr pdb
                next i
                @pds = Twist(@pds[%M-%N],@pds[0],[email protected][0])
                x = @pds XOR shiftr(@pds,11)
                x = x XOR (shiftl(x, 7) AND &h9D2C5680)
                x = x XOR (shiftl(x,15) AND &hEFC60000)
                @pdb = x XOR shiftr(x,18)
             
                Ctx.Idx = 0
             
             end if
          end if
             
          function = [email protected][Ctx.Idx]
          incr Ctx.Idx
          exit function
             
          DoubleMatrixA:
          ! DD &h9908b0df, &h9908b0df
          DoubleUMask:
          ! DD &h80000000, &h80000000
          DoubleLMask:
          ! DD &h7fffffff, &h7fffffff
          Double1:
          ! DD &h00000001, &h00000001
          Tempering1:
          ! DD &h9D2C5680, &h9D2C5680
          Tempering2:
          ! DD &hEFC60000, &hEFC60000
          END FUNCTION
             
             
          #IF %def($DHFILE)
          '-- NOTE: Routine performs no file access error checking.
          '============
          FUNCTION MakeDiehardFile&()
          REGISTER i&, j&
          LOCAL  t$, buffer$, outfile&
          LOCAL ctx AS MT_CONTEXT, pd AS DWORD PTR
             
          '-- Truncate any existing $DHFILE
          outfile = freefile
          open $DHFILE for binary as outfile base = 0
          seteof outfile
          close outfile
             
          ctx.Seed = 5489
          InitMT ctx
             
          for i = 1 to 10
             buffer = zbs(1000000)
             pd = strptr(buffer)
             j = len(buffer)\sizeof(@pd)
             do while j
                @pd = RandomMT(ctx) : incr pd
                decr j
             loop
             open $DHFILE for binary as outfile base = 0
             seek outfile, lof(outfile)
             put$ outfile, buffer
             if i = 10 then
                t = "Size of "+ $DHFILE + ":"+ str$(lof(outfile))
                ShowText(t)
             end if
             close outfile
          next i
          END FUNCTION
          #ENDIF
             
          '===================
          FUNCTION HasMMX&()
          !  mov      eax,  1
          !  cpuid
          !  xor      eax,  eax
          !  test     edx,  &h800000       ;bit 23
          '-- Rem next line to force HasMMX = %FALSE
          !  setnz    al
          !  mov      function, eax
          END FUNCTION
          [This message has been edited by Greg Turgeon (edited February 07, 2005).]

          Comment


          • #6
            Greg, should we expect the MMX version and your original
            version to produce the same sequence given the same seed?
            I get divergence around byte 900 when I compare the two
            starting with the same seed value.


            ------------------

            Comment


            • #7
              should we expect the MMX version and your original
              version to produce the same sequence given the same seed?
              Yes. And you will, now.

              ------------------
              -- gturgeon at verizon dot net --

              Comment

              Working...
              X