FreeBASIC's PRNG #2

[jj2007]

I thought you had PractRand.
Anyway, here it is: PractRand_0.94.zip

Thanks, I had a version that throws exceptions. Now I've run yours and got exceptions again. It seems that PractRand doesn't check for EOF. Oh well... but I got results:

Code: Select all

RNG_test using PractRand version 0.93
RNG = RNG_stdin, seed = 0x4a459127
test set = normal, folding = standard(unknown format)

rng=RNG_stdin, seed=0x4a459127
length= 4 megabytes (2^22 bytes), time= 2.2 seconds
  Test Name                         Raw       Processed     Evaluation
  [Low1/32]BCFN(2+0,13-9,T)         R= +10.5  p =  1.8e-3   unusual
  ...and 98 test result(s) without anomalies

rng=RNG_stdin, seed=0x4a459127
length= 8 megabytes (2^23 bytes), time= 11.0 seconds
  no anomalies in 107 test result(s)

rng=RNG_stdin, seed=0x4a459127
length= 16 megabytes (2^24 bytes), time= 23.4 seconds
  no anomalies in 119 test result(s)

rng=RNG_stdin, seed=0x4a459127
length= 32 megabytes (2^25 bytes), time= 39.6 seconds
  Test Name                         Raw       Processed     Evaluation
  [Low4/32]DC6-9x1Bytes-1           R=  +6.0  p =  3.8e-3   unusual
  ...and 129 test result(s) without anomalies

Note that I used my version 0.93. No crashes until it reaches EOF, and this time I chose a file with exactly 32MB, so it stops there. Your version 0.94 gives me error reading standard input using the same commandline as for version 0.93:

Code: Select all

RNG_test.exe stdin32 <MbRandPcg32.dat

Btw how can your machine be so fast???

length= 2 gigabytes (2^31 bytes), time= 22.7 seconds

[dafhi]

Think I remember deltarho mentioning he had a nice processor. I was thinking just on the smaller periods "that seems kinda fast"

I'm working on a low bits generator i can use to test different algorihms.

[deltarho[1859]]

It seems that PractRand doesn't check for EOF.

I don't have that problem because I don't read a file.

Btw how can your machine be so fast???

Because I am not ready a file.

I do not use 'RNG_test stdin32 <whatever>'.

I use 'My_RNG | RNG_test stdin32'.

The data is piped into RNG_test from MY_RNG.

In my last but one post I gave a link to a post in the thread 'asm .members'. There are two code blocks. The first code block is my piping program. What I do is create a 1MB string and fill it with 262144 * Ulongs. I then Print the string. The Do/Loop has no exit. I then fill the string with fresh data and pipe again. A 4TB test with piping would take over 13 hours on my Intel i7 3.9GHz (with turbo). That is just the 4TB test - we have to add the time for the 2TB, 1TB, 512GB tests ... and so on. I may be wrong but I think PractRand is written to accept no more than 16TB but the most that I have piped is 2TB with my CMWC4096.

[jj2007]

Thanks for the explanation. Yes, that's what I expected - interesting that piping is faster than reading from a file. Fortunately my own (ultra fast) algo passes ENT but fails early with PractRand - I don't want to torture my notebook ;-)

Is there evidence that the new multiplier is better than the old one? And why do we keep several less performant algos in FB if Pcg32 beats them all hands down? Seriously, what is the logic behind that? Wouldn't it be better to abandon the other ones, or do they have other advantages?

[deltarho[1859]]

Is there evidence that the new multiplier is better than the old one?

No. If they are from the same family of primes then there should be no difference. It could be that Melissa O'Neill doesn't even know about it. O'Neill did not know of PractRand until after she published the PCG family. O'Neill was using Small Crush and Big Crush from the TestU01 suite. She now thinks that PractRand is the best thing since sliced bread. <smile>

There are not many multipliers for 64-bit LCGs. Here are two more multipliers that I have found: 3202034522624059733 and 3935559000370003845 from a paper by Pierre L'Ecuyer (1999) which included the one that I found on a different website.

And why do we keep several less performant algos in FB if Pcg32 beats them all hands down? Seriously, what is the logic behind that? Wouldn't it be better to abandon the other ones, or do they have other advantages?

I agree, to a point. I would get rid off FB #1, FB #2 and FB #4. FB #5 should be kept because it is cryptographic - I use it in PCG32II and other projects. Yes, it is slow but if we are only requesting a few numbers it is fast enough. I would also keep FB #3 because it has a very large period and there are uses for such even though it fails PractRand at 256GB. Having said that it could be replaced by my CMWC4096 which has an even larger period, is faster and got to 2TB with PractRand.

PCG32II will work on Windows and Linux and srvaldez is looking at it with regard Mac OS. PCG32II is also thread safe - none of the FB are. If ever it got in FB there would be no need to mention that I had anything to do with it - I am not into that. O'Neill should be credited, of course.

Added: O'Neill does know about Pierre L'Ecuyer's paper mentioned above - it is Ref #25 in her September 5, 2014 paper.

[jj2007]

Have you seen this thread at Hacker News?

[deltarho[1859]]

Have you seen this thread at Hacker News?

Yes. Vigna and O'Neill have been at each other's throats for some time. On the other hand, O'Neill gets on quite well with Vigna's partner in crime, David Blackman. I have seen some of Vigna's comments about O'Neill's work and I found some of them a tad vitriolic. O'Neill's comments at her blog about Vigna's work have been entirely professional and take some of Vigna's algorithms to the cleaners.

Vigna reckoned that no one really needed a PRNG with a period greater than 128-bit. This is when the PCG family was sailing through Practrand whilst Vigna's efforts were failing at 64GB. He then got involved in 256-bit generators and PractRand tore them to shreds as well. The term 'sour grapes' comes to mind.

I have read O'Neill's paper and it is not 'classical' in the sense of papers by L'Ecuyer and Marsaglia and I suppose I would describe it as a modern read. That aside she is clever and is capable of admitting her mistakes; a rare condition amongst many academics who will hold their ground even when standing on thin ice. <smile>

[jj2007]

Yes, that thread is fun to read ;-)

Btw how can your machine be so fast???

length= 2 gigabytes (2^31 bytes), time= 22.7 seconds

I've caught up a little bit:

Code: Select all

rng=RNG_stdin32, seed=unknown
length= 2 gigabytes (2^31 bytes), time= 47.3 seconds
  no anomalies in 204 test result(s)
...
rng=RNG_stdin32, seed=unknown
length= 32 gigabytes (2^35 bytes), time= 681 seconds
  no anomalies in 251 test result(s)  

Still, with my Core i5 I am not fit to indulge in that hobby of yours <grin>
But you forced me to review my WritePipe command, which is a good thing. Now I wonder why my ultrafast (7*pcg32.FB) Rand() fails so early with PractRand, when ENT gave it a clap on the shoulder...

[deltarho[1859]]

One test method protocol:

1) Chi-square: With a failure then do not bother with 2) and 3) because they will result in failure. If a pass then move onto 2)
2) Ent: With a failure then do not bother with 3) because that will result in failure. If a pass then move on to 3)
3) Practrand: With a failure then since we are here because of a 2) pass then the tested generator may still be OK provided it is not used for large amounts of requested random numbers. If a pass then break out the champagne. <smile>

Examples:

A LCG may pass 1) and 2) but will fail 3). For samples less than 4KB a good LCG is very fast and the random numbers will be of good quality. 64-bit is better than 32-bit. If speed is more important than quality then a LCG can be used for larger samples but I don't know at what point the quality should not be ignored.

Mersenne Twister passes 1) and 2) but fails PractRand at 256GB. 128GB of random bytes, a safe limit, is 34,359,738,368 32-bit values. That is a lot of values before we need to be concerned about any bias creeping in. On the other hand, CMWC4096 is faster, has a larger period and passes PractRand up to at least 2TB.

Some will disagree with my LCG comments with regard 4KB but, even though I am a purist, I am not dismissing LCGs.

[paul doe]

Have you seen this thread at Hacker News?

Interesting. There's a link in that discussion that points to an implementation of a PRNG that I've never seen before. It's here.

I ported the code to FreeBasic to test it:

Code: Select all

/'
  Pseudo-Random Number Generator
  
  This is a PRNG found on this thread:
    https://stackoverflow.com/a/27160892
  
  According to its author, it 'passes all general purpose
  statistical tests, but is not cryptographically secure'.
  
  Original code:
  
  static unsigned long long rng_a, rng_b, rng_c, rng_counter;
  unsigned long long rng64() {
      unsigned long long tmp = rng_a + rng_b + rng_counter++;
      rng_a = rng_b ^ (rng_b >> 12);
      rng_b = rng_c + (rng_c << 3);
      rng_c = ((rng_c << 25) | (rng_c >> (64-25))) + tmp;
      return tmp;
  }
  void seed(unsigned long long s) {
      rng_a = rng_b = rng_c = s; rng_counter = 1;
      for (int i = 0; i < 12; i++) rng64();
  }
  
  Wrapped in a class for convenience. Dangling data is not my cup-of-tea,
  and even less if it's STATIC dangling data =D
'/
type PRNG
  public:
    declare constructor()
    declare constructor( byval as ulongint )
    
    declare sub seed( byval as ulongint )
    declare function rnd64() as ulongint
    
  private:
    m_rng_a as ulongint
    m_rng_b as ulongint
    m_rng_c as ulongint
    m_rng_counter as ulongint
end type

constructor PRNG()
  seed( culngint( timer() * ( 2 ^ 64 ) ) ) '' Whatever
end constructor

constructor PRNG( byval s as ulongint )
  seed( s )
end constructor

sub PRNG.seed( byval s as ulongint )
  m_rng_a = s
  m_rng_b = s
  m_rng_c = s
  m_rng_counter = 1
  
  for i as integer = 1 to 12
    rnd64()
  next
end sub

function PRNG.rnd64() as ulongint
  dim as ulongint tmp = m_rng_a + m_rng_b + m_rng_counter
  m_rng_counter += 1
  m_rng_a = m_rng_b xor ( m_rng_b shr 12 )
  m_rng_b = m_rng_c + ( m_rng_c shl 3 )
  m_rng_c = ( ( m_rng_c shl 25 ) or ( m_rng_c shr ( 64 - 25 ) ) ) + tmp
  
  return( tmp )
end function

/'
  Unit test
'/
/'
  This seed was used to test integrity and equivalence with the C code
  You can use whatever you wish (instantiating it with no parameters
  will use the timer() function).
'/
var rn = PRNG( 23 )

for i as integer = 1 to 100
  ? rn.rnd64()
next

sleep()

@deltarho[1859]: What do you think of this one? Testing it as rigurously as you do will take me ages here in crapland =D
According to the person that posted it, it 'passes all general purpose statistical tests, but is not cryptographically secure'. Your opinion on this one?

[deltarho[1859]]

According to the person that posted it, it 'passes all general purpose statistical tests, but is not cryptographically secure'. Your opinion on this one?

In view of your comments, I went straight into PractRand but, just in case, I started at 1KB. I got a clean sweep up to 512KB but then got a spectacular failure at 1MB.

I will now go back to 1) above and work forward.

[paul doe]

In view of your comments, I went straight into PractRand but, just in case, I started at 1KB. I got a clean sweep up to 512KB but then got a spectacular failure at 1MB.

Yeah, I thought so >=D
Don't bother with the others, please. I was only curious by the claims made by the SO post. Since I'm not a crypto/rng specialist, I thought I'd ask someone who is =D

[deltarho[1859]]

Don't bother with the others, please.

I am geared up to do tests quickly.

I'm not too good at following my own advice and used ENT, on 16MB of data (2MB of 64-bits) since it has a chi-square test.

Code: Select all

Chi square distribution for 16777216 samples is 201.98, and randomly
would exceed this value 99.39 percent of the times.

Arithmetic mean value of data bytes is 127.5385 (127.5 = random).
Monte Carlo value for Pi is 3.138319764 (error 0.10 percent).
Serial correlation coefficient is -0.000154 (totally uncorrelated = 0.0).

The last three metrics are OK but the Chi-Square tells us that there is very strong evidence that the data is not random.

I would not write it off. it is better than a LCG and up to 512KB should be safe. It is fast - 520MHz on my kit.

[paul doe]

I would not write it off. it is better than a LCG and up to 512KB should be safe. It is fast - 520MHz on my kit.

Excellent, thanks for your troubles. Scouring the web, I found this one:

Middle Square Weyl Sequence RNG

The paper describes a brand-new PRNG (it was published this very year), so perhaps not enough tests have been performed yet. According to its author, it passes the TestU01 BigCrunch without issues. This site has a C implementation available (32-bit granularity only for now, sadly):

https://mswsrng.wixsite.com/rand

And this is my port of the code, wrapped in an immutable class for convenience:

Code: Select all

/'
  The original code can be found at:
    https://mswsrng.wixsite.com/rand
  
  uint64_t x = 0, w = 0, s = 0xb5ad4eceda1ce2a9;
  
  inline static uint32_t msws() {
  
     x *= x; 
     x += (w += s); 
     return x = (x>>32) | (x<<32);
}
'/
#define hiqword( x ) ( ( cast( ulongint, x ) and &hffffffff00000000 ) shr 32 )
#define loqword( x ) ( cast( ulongint, x ) and &h00000000ffffffff )
#define lsb &b1 '' Least-significant bit

type MsWs
  public:
    declare constructor()
    declare constructor( byval as ulongint )
    
    declare function one() as ulong
    declare property seed() as ulongint
    
  private:
    m_x as ulongint
    m_w as ulongint
    m_s as ulongint
end type

constructor MsWs()
  /'
    According to Bernard Widynski (the creator of this PRNG), the seed
    must be a number that has the upper 32-bits non-zero and the
    least significant bit set to 1. There are 2 ^ 63 numbers that satisfy
    this property, and 25000 such numbers are given in the C implementation.
    Here, I just use a random seed that satisfies that property.
  '/
  randomize( , 5 )
  
  do
    '' Need to use two random numbers as FB MT has 32-bit granularity
    m_s = ( culngint( rnd() * culngint( -1 ) ) shr 32 ) or _
      culngint( rnd() * culngint( -1 ) )
  loop until( ( hiqword( m_s ) <> 0 ) andAlso ( loqword( m_s ) and lsb ) )
  
  this.constructor( m_s )
end constructor

constructor MsWs( byval s as ulongint )
  '' You can of course provide your own seed if you know what you're doing =D
  m_s = s
end constructor

property MsWs.seed() as ulongint
  return( m_s )
end property

function MsWs.one() as ulong
  m_x *= m_x : m_w += m_s : m_x += m_w
  m_x = ( m_x shr 32 ) or ( m_x shl 32 )
  
  return( m_x )
end function

/'
  Unit test
'/
'var randomNumber = MsWs( &hb5ad4eceda1ce2a9 ) <-- To test equivalence with C code
var randomNumber = MsWs()

for i as integer = 1 to 100
  ? randomNumber.one()
next

sleep()

Sorry to waste your time if this is old news for you =D

[deltarho[1859]]

@paul doe

Sorry to waste your time if this is old news for you =D

I hadn't seen this.

I read the paper, looked at your 'function one() as ulong' and figured Widynski was having a laugh especially since von Neumann died in 1957.

The last three metrics of ENT were OK but I got a chi-square of 2.47%. This is equivalent to 97.53% at the other end of the scale giving strong evidence of non-randomness but not as bad as your earlier PRNG. However, I decided to give it the benefit of the doubt and let PractRand loose on it.

I wasn't expecting it to last long so started PractRand at 1KB. I got a lowest ranking anomaly at 128KB and 512KB. We are now bettering a LCG. It got to 64GB without further issue. We are now bettering Vigna and Blackman's xor* generators. A clean sweep followed to 256GB with only a lowest ranking anomaly. We are now betting Mersenne Twister.

By now I was sat in disbelief. The algorithm is tiny and does not look capable of what it was doing. I did a speed test and got 431MHz so it is no slouch. That, by the way, on my machine, is very nearly five times faster than Mersenne Twister.

Back at PractRand. The 512 GB test will take about 90 minutes and I could not see it getting past that. If it fails there then we could have a little gem safe to 256GB. Obviously, further testing is required but it is looking good so far. Had I started PractRand later it would not be showing three anomalies to date.

It has a similarity with PCG32II giving access to 2^63 sequences.

It has just sailed past 512GB! I have another three hours to wait for the 1TB result. My nerves cannot take it!

Two hours, 59.7 minutes later it walked past 1TB.

Well, well, well - what do you know? 1TB and only three lowest ranking anomalies.

Paul, you should get mentioned in dispatches for finding this.

I have booked PractRand for further tests using different seeds and sequences.

Well done!