Mersenne Twister, Public domain

General FreeBASIC programming questions.
srvaldez
Posts: 3379
Joined: Sep 25, 2005 21:54

Mersenne Twister, Public domain

Post by srvaldez »

I found this at https://github.com/ESultanik/mtwister
I will leave the code cleaning as an exercise, but it works as is.

Code: Select all

/' An implementation of the MT19937 Algorithm for the Mersenne Twister
   by Evan Sultanik.  Based upon the pseudocode in: M. Matsumoto and
   T. Nishimura, "Mersenne Twister: A 623-dimensionally
   equidistributed uniform pseudorandom number generator," ACM
   Transactions on Modeling and Computer Simulation Vol. 8, No. 1,
   January pp.3-30 1998.
   
   https://github.com/ESultanik/mtwister
   
   License Public domain. Have fun!
 '/
#pragma once

#include once "crt/long.bi"

const UPPER_MASK = &h80000000
const LOWER_MASK = &h7fffffff
const TEMPERING_MASK_B = &h9d2c5680
const TEMPERING_MASK_C = &hefc60000
#define __MTWISTER_H
const STATE_VECTOR_LENGTH = 624
const STATE_VECTOR_M = 397 ' changes to STATE_VECTOR_LENGTH also require changes to this

type tagMTRand
	mt(0 to 623) as culong
	index as long
end type

type MTRand as tagMTRand
declare function seedRand(byval seed as culong) as MTRand
declare function genRandLong(byval rand as MTRand ptr) as culong
declare function genRand(byval rand as MTRand ptr) as double

private sub m_seedRand(byval rand as MTRand ptr, byval seed as culong)
  /' set initial seeds to mt[STATE_VECTOR_LENGTH] using the generator
     from Line 25 of Table 1 in: Donald Knuth, "The Art of Computer
     Programming," Vol. 2 (2nd Ed.) pp.102.
   '/
	dim as long i
	rand->mt(0) = seed and &hffffffff
	For i = 1 To STATE_VECTOR_LENGTH - 1	
		rand->mt(i) = (6069 * rand->mt((i - 1))) and &hffffffff
	next
	rand->index=i
end sub

/'
  Creates a new random number generator from a given seed.
'/
private function seedRand(byval seed as culong) as MTRand
	dim rand as MTRand
	m_seedRand(@rand, seed)
	return rand
end function

/'
   Generates a pseudo-randomly generated long.
 '/
private function genRandLong(byval rand as MTRand ptr) as culong
	dim y as culong
	static mag(0 to 1) as culong = {&h0, &h9908b0df} ' mag[x] = x * 0x9908b0df for x = 0,1
	if (rand->index >= 624) orelse (rand->index < 0) then
		' generate STATE_VECTOR_LENGTH words at a time
		dim kk as long
		if (rand->index >= (624 + 1)) orelse (rand->index < 0) then
			m_seedRand(rand, 4357)
		end if
		kk = 0
		Do While kk<STATE_VECTOR_LENGTH-STATE_VECTOR_M
			y = (rand->mt(kk) and &h80000000) or (rand->mt((kk + 1)) and &h7fffffff)
			rand->mt(kk) = (rand->mt((kk + 397)) xor (y shr 1)) xor mag((y and &h1))
			kk += 1
		loop
		Do While kk<STATE_VECTOR_LENGTH-1
			y = (rand->mt(kk) and &h80000000) or (rand->mt((kk + 1)) and &h7fffffff)
			rand->mt(kk) = (rand->mt((kk + (397 - 624))) xor (y shr 1)) xor mag((y and &h1))
			kk += 1
		loop
		y = (rand->mt((624 - 1)) and &h80000000) or (rand->mt(0) and &h7fffffff)
		rand->mt((624 - 1)) = (rand->mt((397 - 1)) xor (y shr 1)) xor mag((y and &h1))
		rand->index = 0
	end if
	y = rand->mt(rand->index)
	rand->index += 1
	y xor= y shr 11
	y xor= (y shl 7) and &h9d2c5680
	y xor= (y shl 15) and &hefc60000
	y xor= y shr 18
	return y
end function

/'
   Generates a pseudo-randomly generated double in the range [0..1].
 '/
function genRand(byval rand as MTRand ptr) as double
	return cdbl(cdbl(genRandLong((rand))) / cast(culong, &hffffffff))
end function

'' end of mtwister

'' simple test
dim as MTRand r = seedRand(1337)
for i as long=0 to 20
	print genRand(@r)
next
deltarho[1859]
Posts: 4313
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: Mersenne Twister, Public domain

Post by deltarho[1859] »

Mersenne Twister (MT)

1) Slow compared to modern generators.
2) Fails PractRand at 256GB. Modern generators get to 16TB.
3) State vector is huge compared with modern generators.
4) The above implementation seeds the state vector with a Linear Congruential Generator which in turn is seeded with 32-bit. This limits the number of entry points to MT's period to a tiny fraction of the possible entry points.

There is a version of MT on this forum which took Greg Lyon's FreeBASIC implementation and replaced the BASIC engine with assembler. For fixed seeding a string is used and hashed enough times using SHA512 to fully populate the state vector. For random seeding the state vector is fully populated using cryptographic random numbers.

The resulting code is called RndMT and published by 'muggins' here in July 2017.

When benchmarking generators I no longer include RndMT.

Many implementations of MT in programming languages, applications and you name it have been moving away from MT for some time.

MT: RIP

Sorry, srvaldez.
srvaldez
Posts: 3379
Joined: Sep 25, 2005 21:54

Re: Mersenne Twister, Public domain

Post by srvaldez »

@deltarho[1859]
have checked out Shishua PRNG https://www.mdeditor.tw/pl/pij1/zh-hk ?
also ChaCha https://cr.yp.to/chacha.html
deltarho[1859]
Posts: 4313
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: Mersenne Twister, Public domain

Post by deltarho[1859] »

On the forum is code for the following that I have written.

PCG32II
MsWs
CryptoRndII
CMWC4096
xoroshiro128**
xoshiro256**
drSquares4

They are all PractRand resistant and some are very fast.

Did you not find one to your liking assuming that you have checked some of them out?

With regard Shishua I am a bit wary of people who talk about Chi squared and Die Hard although PractRand was eventually discovered. Something of PractRand was learnt: "It cannot rate high-quality generators, making comparisons between them impossible." I have said as much myself.

I look forward to your implementation of Shishua. Coding the engine is one thing - giving the user some useful functions is another.
srvaldez
Posts: 3379
Joined: Sep 25, 2005 21:54

Re: Mersenne Twister, Public domain

Post by srvaldez »

deltarho[1859] wrote: Did you not find one to your liking assuming that you have checked some of them out?
...
I look forward to your implementation of Shishua. Coding the engine is one thing - giving the user some useful functions is another.
I gather that random numbers are used in games, statistics, cryptography and scientific simulations but my level of programming does not fit in any of them, actually I do almost no programming these days
translating Shishua to FB might be an interesting learning experience but it will take time because I am not familiar with the _mm256 instructions used in the C source, I will try and see if it will build with gcc
deltarho[1859]
Posts: 4313
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: Mersenne Twister, Public domain

Post by deltarho[1859] »

The _mm256 instructions are from the AVX (Advanced Vector Extensions) architecture. They are effectively the 'rich man's' Streaming SIMD Extesions where we move from 128-bit to 256-bit and, seven years ago, 512-bit. I, too, am not familiar with AVX.

We are not talking mind-boggling floating-point granularity but the ability to do many integer calculations, for example, in parallel.
I do almost no programming these days
That is a shame for us especially since you joined the forum when FreeBASIC was still in short trousers. Image
srvaldez
Posts: 3379
Joined: Sep 25, 2005 21:54

Re: Mersenne Twister, Public domain

Post by srvaldez »

@deltarho[1859]
does the CPU in your PC support AVX2 ?
Wikipedia gives a list https://en.wikipedia.org/wiki/Advanced_ ... _with_AVX2
but besides that, how do you do memory alignment in FB ?
I think it needs to be 64-bit aligned
deltarho[1859]
Posts: 4313
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: Mersenne Twister, Public domain

Post by deltarho[1859] »

srvaldez wrote:does the CPU in your PC support AVX2 ?
No.

Here we have a guy with a vectorized version of PCG32 which outputs twice as fast as the standard C version. Unfortunately he has AVX-512 so that leaves me out. Image

Therein lies a problem. Writing vectorized generators may have a small target user base. This is probably why folks like Melissa O'Neill and Sebastiano Vigna do not have AVX versions of their generators. It may be a different story in the future when we all have AVX-512.
srvaldez
Posts: 3379
Joined: Sep 25, 2005 21:54

Re: Mersenne Twister, Public domain

Post by srvaldez »

my CPU supports AVX2 but not AVX-512
I am sure the Shishua PRNG could be implemented without the use of AVX2 but with slower performance
deltarho[1859]
Posts: 4313
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: Mersenne Twister, Public domain

Post by deltarho[1859] »

srvaldez wrote:I am sure the Shishua PRNG could be implemented without the use of AVX2 but with slower performance
My guess is that it will drop like a stone. It is a bit like a multi-threaded application specifically designed for a multi-core machine. Yes, it will work on a single core machine but will be very much slower.
srvaldez
Posts: 3379
Joined: Sep 25, 2005 21:54

Re: Mersenne Twister, Public domain

Post by srvaldez »

did a search on AVX2 emulation and found that Intel has an emulator https://software.intel.com/content/www/ ... lator.html
which is OK if you want to run a program that uses instruction that your CPU doesn't support, but what I really was after was a library that would emulate those instructions, found this https://blog.adafruit.com/2020/08/07/si ... n-library/
looks like the instructions for AVX2 and above are a work in progress but it looks promising
dodicat
Posts: 7983
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: Mersenne Twister, Public domain

Post by dodicat »

Hi srvaldez.
You must be a schoolteacher:
"I will leave the code cleaning as an exercise"
Your twister, because is regurgitates a long, is very fast for a random range, 10 times faster than fb twister with gcc optimisation.

Code: Select all

/' An implementation of the MT19937 Algorithm for the Mersenne Twister
   by Evan Sultanik.  Based upon the pseudocode in: M. Matsumoto and
   T. Nishimura, "Mersenne Twister: A 623-dimensionally
   equidistributed uniform pseudorandom number generator," ACM
   Transactions on Modeling and Computer Simulation Vol. 8, No. 1,
   January pp.3-30 1998.
   
   https://github.com/ESultanik/mtwister
   
   License Public domain. Have fun!
 '/
#pragma once

#include once "crt/long.bi"

const UPPER_MASK = &h80000000
const LOWER_MASK = &h7fffffff
const TEMPERING_MASK_B = &h9d2c5680
const TEMPERING_MASK_C = &hefc60000
#define __MTWISTER_H
const STATE_VECTOR_LENGTH = 624
const STATE_VECTOR_M = 397 ' changes to STATE_VECTOR_LENGTH also require changes to this

type tagMTRand
   mt(0 to 623) as culong
   index as long
end type

type MTRand as tagMTRand
declare function seedRand(byval seed as culong) as MTRand
declare function genRandLong(byval rand as MTRand ptr) as culong
declare function genRand(byval rand as MTRand ptr) as double

private sub m_seedRand(byval rand as MTRand ptr, byval seed as culong)
  /' set initial seeds to mt[STATE_VECTOR_LENGTH] using the generator
     from Line 25 of Table 1 in: Donald Knuth, "The Art of Computer
     Programming," Vol. 2 (2nd Ed.) pp.102.
   '/
   dim as long i
   rand->mt(0) = seed and &hffffffff
   For i = 1 To STATE_VECTOR_LENGTH - 1   
      rand->mt(i) = (6069 * rand->mt((i - 1))) and &hffffffff
   next
   rand->index=i
end sub

/'
  Creates a new random number generator from a given seed.
'/
private function seedRand(byval seed as culong) as MTRand
   dim rand as MTRand
   m_seedRand(@rand, seed)
   return rand
end function

/'
   Generates a pseudo-randomly generated long.
 '/
private function genRandLong(byval rand as MTRand ptr) as culong
   dim y as culong
   static mag(0 to 1) as culong = {&h0, &h9908b0df} ' mag[x] = x * 0x9908b0df for x = 0,1
   if (rand->index >= 624) orelse (rand->index < 0) then
      ' generate STATE_VECTOR_LENGTH words at a time
      dim kk as long
      if (rand->index >= (624 + 1)) orelse (rand->index < 0) then
         m_seedRand(rand, 4357)
      end if
      kk = 0
      Do While kk<STATE_VECTOR_LENGTH-STATE_VECTOR_M
         y = (rand->mt(kk) and &h80000000) or (rand->mt((kk + 1)) and &h7fffffff)
         rand->mt(kk) = (rand->mt((kk + 397)) xor (y shr 1)) xor mag((y and &h1))
         kk += 1
      loop
      Do While kk<STATE_VECTOR_LENGTH-1
         y = (rand->mt(kk) and &h80000000) or (rand->mt((kk + 1)) and &h7fffffff)
         rand->mt(kk) = (rand->mt((kk + (397 - 624))) xor (y shr 1)) xor mag((y and &h1))
         kk += 1
      loop
      y = (rand->mt((624 - 1)) and &h80000000) or (rand->mt(0) and &h7fffffff)
      rand->mt((624 - 1)) = (rand->mt((397 - 1)) xor (y shr 1)) xor mag((y and &h1))
      rand->index = 0
   end if
   y = rand->mt(rand->index)
   rand->index += 1
   y xor= y shr 11
   y xor= (y shl 7) and &h9d2c5680
   y xor= (y shl 15) and &hefc60000
   y xor= y shr 18
   return y
end function

dim shared as MTRand r

/'
   Generates a pseudo-randomly generated double in the range [0..1].
 '/
function genRand(byval rand as MTRand ptr) as double
   return cdbl(cdbl(genRandLong((rand))) / cast(culong, &hffffffff))
end function



#define rangeI(f,l) (genRandLong(@r) mod (((l)-(f))+1)) + (f)
#define range(f,l) Int(Rnd*(((l)+1)-(f))+(f))
'' end of mtwister

'' simple test
 
 
 dim as double t
 dim as long L
 for k as long=1 to 5
 r = seedRand(1337)
 t=timer
 for n as long=1 to 10000000
     L=rangeI(-10,10)
next

print timer-t,L,"srvaldez"
randomize 2,3
t=timer
 for n as long=1 to 10000000
     L=range(-10,10)
next

print timer-t,L,"fb"
print
next k
sleep

 
Personally I use random ranges far more often than decimal random numbers.
Er . . . sorry sir, I haven't cleaned the code yet.
srvaldez
Posts: 3379
Joined: Sep 25, 2005 21:54

Re: Mersenne Twister, Public domain

Post by srvaldez »

:smile:
deltarho[1859]
Posts: 4313
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: Mersenne Twister, Public domain

Post by deltarho[1859] »

@srvaldez

I reckon this AVX2 emulation lark is right up your street but it ain't mine. Image
srvaldez
Posts: 3379
Joined: Sep 25, 2005 21:54

Re: Mersenne Twister, Public domain

Post by srvaldez »

@dodicat
if you compile for loops with gcc, gcc will eliminate the loop entirely if possible with -O3
so you must add code that will prevent gcc from eliminating these timing loop, so here's your test with my mods

Code: Select all

'' simple test
 
 
 dim as double t, t2
 dim as long L
 dim as double sum
 print " time   factor rnd  sum"
 for k as long=1 to 5
 r = seedRand(1337)
 t=timer
 sum=0
 for n as integer=1 to 10000000
     L=rangeI(-10,10)
     sum+=L
next
t=timer-t
'print t,0.0,L,sum,"srvaldez"
print using "##.#### ##.## ### ###### srvaldez";t;1;L;sum
randomize 2,3
t2=timer
sum=0
 for n as integer=1 to 10000000
     L=range(-10,10)
     sum+=L
next
t2=timer-t2
print using "##.#### ##.## ### ###### fb";t2;t2/t;L;sum
print
next k
sleep
Post Reply