1
I know it’s rediscovering the wheel, but I want to know how I can implement in my own way a nice linear congruent generator to draw integer numbers from 0 to 2^B-1. How to choose constants, how many bits to the sequence, how to make that sequence be well shuffled, etc.
2
First, know that Gcls (linear congruent/congruent generators) are generators of recursive sequences of integer elements in which applies to each element to discover the next of the sequence two linear operations, one multiplication and one addition, and then restrict the result to an interval by calculating split remainder. This is used to generate pseudo-random numbers by forming well-shuffled sequences. How?
Given a m ∈ ℕ*, one f ∈ {1,2,...,m-1} and a t ∈ {0,1,2,...,m-1}, we have a recursion s[n] = ( s[n-1] * f + t ) % m. Adopting an initial value (seed) s[0] ∈ {0,1,...,m-1}, we can from it calculate s[1] = ( s[0] * f + t ) % m, of it calculate s[2] = ( s[1] * f + t ) % m and so on. Each time a new element of the sequence is obtained, it can be modified so that it becomes a more suitable integer as a pseudo-random result. It is common to ignore the least significant bits because they usually have a clear pattern of behavior.
Visual studio’s C library has already used recursion s[n] = ( s[n-1]*214013+2531011 )%4294967296 returning the index bits 16 to 30. One way to implement this is to use a 32-bit integer global state variable int RandomState ; to store the current element of the sequence, implement a Seed definition void srand( int seed ){ RandomState = seed ; } and the draw that updates the state to the next element and takes the correct bits.
int rand(){
//RandomSeed = RandomSeed * 214013 + 2531011 ;
RandomSeed = RandomSeed * 0x343FD + 0x269EC3 ;
return ( RandomSeed & 0x7FFFFFFF )>>16 ;
}
It was not necessary to calculate the rest of the division because, in addition to the 32-bit integer type already automatically limit to 4294967296 numbers in the sequence (although the signal of the high binary codes render the interval of incorrect, not from 0 to 4294967295 but from -2147483648 to 2147483647), the operator & is to undo the index bit 31, which takes the signal and makes it possible with the offset to finish picking up the index bits 16 to 30. With this, the draw results in numbers from 0 to 32767.
Of course, not only in this but in any GCL, if always use the same seed the sequence will always be the same, it does not simulate a random behavior. Therefore it is recommended to use them in such a way that it is guaranteed different seeds with a hidden mechanism, as the use of a clock function, because at each time interval of the granularity of the function the result of the clock is different and each time the result of it changes (soon the seed changes) the sequence also changes.
It can also be used as seed a kind of "genetic code" that determines which will be the sequence so that a procedure with this sequence always forms the same results to form something fixed associated with this "genetic code" and varied for the different codes.
https://en.wikipedia.org/wiki/Linear_congruential_generator https://en.wikipedia.org/wiki/Linear_congruential_generator
First, know that Gcls have a state variable s which houses an integer number, usually with a bit number b greater than B. After determining an initial value (seed), each time this variable is updated a formula such as s=s*f+t where f is an entire factor and t is an integer term, both constant throughout the sequence and defined its choice, but the choice of this pair strongly affects the quality of the results.
It is recommended that these two constants meet the f%4 == 1 and t%2 == 1 so that any number accepted by the whole type of s is achievable in any sequence and it is possible for the draw to be reasonably realistic (look for the best and worst pairs of f and t accordingly).
Another alternative is f%8 == 3 || f%8 == 5 and t == 0 to take out the need of addition, but requires odd seed to function well, in this case the sequence only has odd numbers (half of the numbers supported by the integer type used), only one half or the other (depending on the seed) of these odd numbers are accessible (i.e., there is the partition of its distinct sequences, so the number of reachable values per sequence is 1/4 of those supported by the integer type used) and apparently the most realistic dozens of bit draws have t%2 == 1.
Second, note that the first bits usually follow a very predictable pattern, so the reason to recommend state variable with more bits than the ones you want to draw is exactly to move them. If the pseudo-random sequence is a variable of b bits, then to draw B bits must move b-B bits to the right. Remember to pay attention to the integer signal when scrolling. High values of b to have a larger, less repetitive sequence, but this little increases the local quality potential of the draw.
Third, to draw with this quality one does not have a convention of how to measure the realism of the draws, but it is a fact that Lics simulate discrete uniform models and this type of distribution has characteristics such as mean value, mean deviation, maximum deviation, standard deviation... therefore these properties should be taken into account to the maximum. From there, you can define whether you think the chosen constants are good enough or not, decide whether to use or check others.
One way to measure quality is to experimentally calculate consecutive sequence values, use them to estimate some derivatives (from the index function that corresponds to the sequence element) with finite differences and then compare with the expected values to estimate a kind of deviation from them.
For example, one can calculate s[0]=0, s[1]=s[0]*f+t, s[2]=s[1]*f+t and from there every time a new s[n]=s[n-1]*f+t you have a sequence s[n-3], s[n-2], s[n-1], s[n], which allows up to third order of estimated derivatives using finite differences. In case of b small bits (such as 16), one can calculate the entire sequence in reasonable time with executables that do not have major performance problems, such as compiling code in C maybe without even needing optimization.
The sequence being in s of single values of 0 to 2^b-1 pseudo-randomly distributed, one can obtain another sequence by applying a factor and a term while maintaining uniform properties. Apply v[n]=s[n]/2^b+1/2^(b+1) results in sequence in v of single discrete values evenly distributed between 0 and 1. The sequence v has average 0.5 and average absolute deviation of approximately 0.25. But derivatives have average 0 and miscellaneous absolute mean deviations as a function of b, f and t, the first order of which is expected 291539/320760 ≈ 0.909, second 1919/1200 ≈ 1.599 and third 142/135 ≈ 1.052. If, for example, you calculate derivatives 0.9, 1.6, 1.06 you realize that the relative error of each is 0.0098, 0.0005, 0.0077, therefore one can admit maximum error of 0.98% or minimum quality 100%-0.98% (99.02%).
To apply this quality criterion, the formula may be used 11*v[n]-18*v[n-1]+9*v[n-2]-2*v[n-3] to estimate 6x the derivative in the first order (expected 6*291539/320760 ≈ 5.453), 2*v[n]-5*v[n-1]+4*v[n-2]-v[n-3] to estimate 1x drift to second (≈1.599) and v[n]-3*v[n-1]+3*v[n-2]-v[n-3] to estimate 1x drift to third (≈1.052). With these formulas, using b=4, f=9, t=3 ou t=5 The GCL quality is 89.98% throughout the sequence, but b=4, f=1, t=1 ou t=15 has quality 23.45%.
This C code has the function qc implemented to calculate the quality of the draw, receiving f%4 == 1 and t%2 == 1 and using global variables m defined each time a b=3,4,5,6,7,8 in a loop, thus calculating the qualities and printing the best settings of f and t: https://ideone.com/gHg3T1. Realize that for each b travels f = 1, 5, 9, ..., 2^b-3, internally t = 1, 3, 5, ..., 2^b-1 and more internally (in the function) draws the entire possible sequence, which weighs heavily if using dozens of bits.
But this C code has been updated to VC++ code that accepts any f and t. Using
# include <time.h>
# include <stdio.h>
typedef __int8 Char8 ;
typedef double Float64 ;
typedef __int64 SInt64 ;
typedef __int32 SInt32 ;
typedef unsigned __int8 UInt8 ;
typedef unsigned __int32 UInt32 ;
typedef unsigned __int64 UInt64 ;
and implementing the function
Float64 RandomQuality_Percent( UInt32 Factor , UInt32 Term , SInt64 Bits , UInt32 *MinState , UInt32 *MaxState , Float64 *StatesCount ){
const UInt32 CBits=UInt8( Bits>0 ?( Bits<32 ? Bits : 32 ): 0 ) , AndArg=UInt32( (1LL<<CBits)-1 ) ;
const Float64 StvFactor=1.0/(1LL<<CBits) , StvTerm=( Term==0 ? 0.0 : 0.5*StvFactor ) , MaxStatesCount=AndArg+1.0 ;
UInt32 State=( *MinState = *MaxState = (Factor+Term)&AndArg ) , FinalState=(State*Factor+Term)&AndArg ;
Float64 Value0 , Value1=StvFactor+StvTerm , Value2=State*StvFactor+StvTerm , Temporary ;
Float64 Mrd0=(*StatesCount=0) , Mrd1=0 , Mrd2=0 , Value3=(State=FinalState)*StvFactor+StvTerm ;
do {
(*StatesCount)++ ;
Value0 = Value1 ;
Value1 = Value2 ;
Value2 = Value3 ;
Value3 = ( State=(State*Factor+Term)&AndArg )*StvFactor+StvTerm ;
if( State < *MinState ) *MinState = State ;
if( State > *MaxState ) *MaxState = State ;
Temporary = 1100/5.453404414515526 * Value3 ;
Temporary -= 1800/5.453404414515526 * Value2 ;
Temporary += 900/5.453404414515526 * Value1 ;
Temporary -= 200/5.453404414515526 * Value0 ;
( Temporary<0 ? Mrd0-=Temporary : Mrd0+=Temporary ) ;
Temporary = 200/1.5991666666666667 * Value3 ;
Temporary -= 500/1.5991666666666667 * Value2 ;
Temporary += 400/1.5991666666666667 * Value1 ;
Temporary -= 100/1.5991666666666667 * Value0 ;
( Temporary<0 ? Mrd1-=Temporary : Mrd1+=Temporary ) ;
Temporary = 100/1.0518518518518519 * Value3 ;
Temporary -= 300/1.0518518518518519 * Value2 ;
Temporary += 300/1.0518518518518519 * Value1 ;
Temporary -= 100/1.0518518518518519 * Value0 ;
( Temporary<0 ? Mrd2-=Temporary : Mrd2+=Temporary ) ;
} while(( State != FinalState )&&( MaxStatesCount != *StatesCount )) ;
if( (Mrd0/=*StatesCount) > 100 ) Mrd0 = 200 - Mrd0 ;
if( (Mrd1/=*StatesCount) > 100 ) Mrd1 = 200 - Mrd1 ;
Temporary = ( Mrd0<Mrd1 ? Mrd0 : Mrd1 ) ;
if( (Mrd2/=*StatesCount) > 100 ) Mrd2 = 200 - Mrd2 ;
return ( Temporary<Mrd2 ? Temporary : Mrd2 ) ;
}
to calculate the quality using seed 1 and also store minimum, maximum values of the sequence and counting of elements of it, one can thus obtain more data and for more configurations. This function main search for the best in b=3,4,5,...,12,13.
SInt32 main( UInt32 argsLength , Char8** args ){
Float64 TimeFactor=0 ;
for( UInt8 Bits=2 ; Bits<=13 ; Bits++ ){
UInt32 StartTime=time(0) ;
printf( " >>> BITS=%u Time=?????.??s" , UInt32(Bits) ) ;
UInt32 MaxValuePer4=1<<( Bits-2 ) , MaxValuePer2=MaxValuePer4<<1 ;
UInt64 MaxValuePlusOne = 2LL * MaxValuePer2 ;
Float64 Count=Float64(MaxValuePlusOne)*MaxValuePlusOne , BestQuality=0 , BestStatesCount[999] ;
UInt32 BestLength=0 , BestFactor[999] , BestTerm[999] , BestMinState[999] , BestMaxState[999] ;
for( UInt32 BlockFactor=0 ; BlockFactor<MaxValuePer4 ; BlockFactor++ ){
for( UInt32 BlockTerm=0 ; BlockTerm<MaxValuePer2 ; BlockTerm++ ){
for( UInt32 PieceFactor=0 ; PieceFactor<4 ; PieceFactor++ ){
UInt32 Factor = 4 * BlockFactor + PieceFactor ;
for( UInt32 PieceTerm=0 ; PieceTerm<2 ; PieceTerm++ ){
UInt32 Term = 2 * BlockTerm + PieceTerm ;
printf( "\b\b\b\b\b\b\b\b\b%8.2lfs" , TimeFactor * Count-- ) ;
Float64 Quality = 1e-9*UInt64( 0.5+1e9*RandomQuality_Percent( Factor , Term , Bits , BestMinState+BestLength , BestMaxState+BestLength , BestStatesCount+BestLength ) ) ;
if(( Quality >= BestQuality )&&( Quality < 100 )&&( BestStatesCount[BestLength] >= 4 )){
if( Quality > BestQuality ){
BestQuality = Quality ;
BestLength = 0 ;
}
BestFactor[ BestLength ] = Factor ;
BestTerm[ BestLength++ ] = Term ;
}
}
}
}
}
TimeFactor = Float64( time(0)-StartTime )/(( MaxValuePlusOne )*Float64( MaxValuePlusOne )) ;
printf( "\b\b\b\b\b\b\b\b\b\b\b\b\b\bQuality=%.9lf%%\n" , BestQuality ) ;
for( UInt32 index=0 ; index<BestLength ; index++ ){
printf( " --> FACTOR(%.4lf)=%u " , BestFactor[index]/Float64(MaxValuePlusOne) , BestFactor[index] ) ;
printf( "TERM(%.4lf)=%u " , BestTerm[index]/Float64(MaxValuePlusOne) , BestTerm[index] ) ;
printf( "MinState=%u MaxState=%u StatesCount=%.0lf\n" , BestMinState[index] , BestMaxState[index] , BestStatesCount[index] ) ;
}
}
}
Notice that there are more loops than necessary, because there is the composition of term and factor with two components, the block and the piece. It is to stabilize the time count, because there are neighboring values that take different times. In addition, the execution takes time. It would be a good thing to parallelize the inner cycles.
The result printed on the console was as follows.
>>> BITS=2 Quality=91.685846491%
--> FACTOR(0.2500)=1 TERM(0.2500)=1 MinState=0 MaxState=3 StatesCount=4
--> FACTOR(0.2500)=1 TERM(0.7500)=3 MinState=0 MaxState=3 StatesCount=4
>>> BITS=3 Quality=91.685846491%
--> FACTOR(0.1250)=1 TERM(0.2500)=2 MinState=0 MaxState=7 StatesCount=8
--> FACTOR(0.1250)=1 TERM(0.7500)=6 MinState=1 MaxState=7 StatesCount=4
--> FACTOR(0.6250)=5 TERM(0.2500)=2 MinState=1 MaxState=7 StatesCount=4
--> FACTOR(0.6250)=5 TERM(0.7500)=6 MinState=1 MaxState=7 StatesCount=4
>>> BITS=4 Quality=91.685846491%
--> FACTOR(0.0625)=1 TERM(0.2500)=4 MinState=0 MaxState=15 StatesCount=16
--> FACTOR(0.0625)=1 TERM(0.7500)=12 MinState=1 MaxState=13 StatesCount=4
--> FACTOR(0.3125)=5 TERM(0.0000)=0 MinState=1 MaxState=13 StatesCount=4
--> FACTOR(0.3125)=5 TERM(0.5000)=8 MinState=1 MaxState=13 StatesCount=4
--> FACTOR(0.5625)=9 TERM(0.2500)=4 MinState=1 MaxState=13 StatesCount=4
--> FACTOR(0.5625)=9 TERM(0.7500)=12 MinState=1 MaxState=13 StatesCount=4
--> FACTOR(0.8125)=13 TERM(0.0000)=0 MinState=1 MaxState=13 StatesCount=4
--> FACTOR(0.8125)=13 TERM(0.5000)=8 MinState=1 MaxState=13 StatesCount=4
>>> BITS=5 Quality=98.987676056%
--> FACTOR(0.9062)=29 TERM(0.1562)=5 MinState=0 MaxState=31 StatesCount=32
--> FACTOR(0.9062)=29 TERM(0.7188)=23 MinState=0 MaxState=31 StatesCount=32
>>> BITS=6 Quality=99.145568860%
--> FACTOR(0.5781)=37 TERM(0.0469)=3 MinState=0 MaxState=63 StatesCount=64
--> FACTOR(0.5781)=37 TERM(0.5156)=33 MinState=0 MaxState=63 StatesCount=64
>>> BITS=7 Quality=99.850182387%
--> FACTOR(0.1641)=21 TERM(0.3828)=49 MinState=0 MaxState=127 StatesCount=128
--> FACTOR(0.1641)=21 TERM(0.7734)=99 MinState=0 MaxState=127 StatesCount=128
--> FACTOR(0.7891)=101 TERM(0.3828)=49 MinState=0 MaxState=127 StatesCount=128
--> FACTOR(0.7891)=101 TERM(0.3984)=51 MinState=0 MaxState=127 StatesCount=128
>>> BITS=8 Quality=99.958983357%
--> FACTOR(0.5820)=149 TERM(0.2461)=63 MinState=0 MaxState=255 StatesCount=256
--> FACTOR(0.5820)=149 TERM(0.3320)=85 MinState=0 MaxState=255 StatesCount=256
>>> BITS=9 Quality=99.985729809%
--> FACTOR(0.9629)=493 TERM(0.4043)=207 MinState=0 MaxState=511 StatesCount=512
--> FACTOR(0.9629)=493 TERM(0.5566)=285 MinState=0 MaxState=511 StatesCount=512
>>> BITS=10 Quality=99.996862277%
--> FACTOR(0.3018)=309 TERM(0.3291)=337 MinState=0 MaxState=1023 StatesCount=1024
--> FACTOR(0.3018)=309 TERM(0.9717)=995 MinState=0 MaxState=1023 StatesCount=1024
>>> BITS=11 Quality=99.997289187%
--> FACTOR(0.4126)=845 TERM(0.4468)=915 MinState=0 MaxState=2047 StatesCount=2048
--> FACTOR(0.4126)=845 TERM(0.9653)=1977 MinState=0 MaxState=2047 StatesCount=2048
>>> BITS=12 Quality=99.999219599%
--> FACTOR(0.5071)=2077 TERM(0.6321)=2589 MinState=0 MaxState=4095 StatesCount=4096
--> FACTOR(0.5071)=2077 TERM(0.8748)=3583 MinState=0 MaxState=4095 StatesCount=4096
>>> BITS=13 Quality=99.999718263%
--> FACTOR(0.8717)=7141 TERM(0.2313)=1895 MinState=0 MaxState=8191 StatesCount=8192
--> FACTOR(0.8717)=7141 TERM(0.6403)=5245 MinState=0 MaxState=8191 StatesCount=8192
With a similar program, I was able for a few hours to track values of f form configurations b=32, t=±1, f%4==1 of sequences such that the quality in each complete one I could calculate in about a minute.
With 64 bits, the quality is not calculated in the entire sequence nor in centuries. No other reasonable criterion for measuring was defined. Soon there was no such attempt. Also, even limited to 32 bits no standard was found for values of t favoring higher quality sequence encounters. Because of this, increment was accepted (t=1) and decrease (t=-1 or t=4294967295).
Examining each f=5°,5¹,5²,5³,...,(5^759)%4294967296 with t ∈ {1,-1}, until then the best formula found was s[n+1] = (2335288753*s[n]-1)%4294967296 = 0x8B31ADB1*s[n]+0xFFFFFFFF with the 53rd value of f, t=-1 and quality 99.9999999751%.
I can only speak for myself. I found the question and the answer interesting, but if the person does not know anything about linear congruent generators and wants to know the algorithm will not be able to extract anything. Just my opinion, I may be mistaken, but I could introduce the subject before I get to the point.
I noticed an error in the code that changes some conclusions about t=0. I fixed everything and added more details and textual structure, I’ll just wait to see if you finish running the program I did to present the result.
Another PRNG algorithm doesn’t interest you, like Accorn?
Browser other questions tagged random random-numbers
You are not signed in. Login or sign up in order to post.
If someone posts a better, simpler and more effective answer, I accept as a solution.
– RHER WOLF