Shop OBEX P1 Docs P2 Docs Learn Events
Random/LFSR on P2 - Page 84 — Parallax Forums

Random/LFSR on P2

1818284868792

Comments

  • xoroshironotxoroshironot Posts: 309
    edited 2020-01-12 04:32
    evanh wrote: »
    They're only warnings so compile and link still succeeds. I haven't tried using it yet.
    Yes, I eventually got the compile to succeed, but then found the results were missing some entries that were present in 0.94... some of that was intentional, but some was not and needed the source to be hand modified.
    All the notes are in that discussion.
    I hope it is finalized soon.
  • xoroshironotxoroshironot Posts: 309
    edited 2020-01-12 17:45
    evanh wrote: »
    Strong scores but has a weak spot with rotated 8-bit window. Will need search for a better candidate.
    I am unable to replicate your 8-bit or 16-bit window results using this scrambler seeded { 1, 0, 0 }: Edit: I found you are using rotr instead of rotl, which allows me to replicate your results.
    state[2] += rotl(s0 + s1, 9) + 1;
    
    In the case of the 8-bit window, I am attempting an easy fail at 32MB in PractRand 0.94 by:
    1. Rotating the 16-bit output by 11 bits.
    2. Concatenating the low 8 bits with those from the next output shifted left 8 bits.
    3. Buffer the result from #2 and repeat.
    4. Using this code with a 16-bit buffer array:
    buffer[j] = (rotl(prng(), 11) & 0xff) | (rotl(prng(), 11) << 8);
    
    I also tried a single element char buffer:
    buffer[0] = rotl(prng(), 11) & 0xff;
    
    Both with same 4GB fail result.

    Do you see something obviously wrong in my approach (or in your code)?

    Edit: I found my issue. See above.
  • evanhevanh Posts: 16,007
    Heh, good. I would have been scratching my head to help. It feels like another life when I wrote that.

  • evanhevanh Posts: 16,007
    edited 2020-01-12 20:45
    I only just remembered why so many of my old scripts are wildly out of date. Tony's approach of using frequency distribution sped up the culling process a lot.

    1: Identify full period engine candidates.
    2: Run distribution tests on all those.
    3: Select a decent number of top distribution scoring to run gridding with Practrand.
    4: The best few of those can then be Crush tested.
  • xoroshironotxoroshironot Posts: 309
    edited 2020-01-12 23:21
    I ran lots of distribution tests on various candidates (attached to several of my previous posts), but none on the one candidate you focused on with the missing + s0.
    That simplest/fastest candidate using three 32-bit state variables is not expected to show any weakness until after 256TB... but I have not finished testing it to be sure it gets even that far.
    I was already at least vaguely aware that it was unsuitable for 16-bit state variables, which is another reason why Tony and I put little, if any, effort into discussing it.
    Using 16-bit state variables, I will stand by the [13 5 10 9] scrambler:
    state[2] += rotl(s0 + s1, 9) + s0 + 1; // toggle + 1 on/off each iteration to increase equidistribution
    
    ... which should fail at 8GB under the same conditions that caused the 'no + s0' version to fail at 32MB.
    At other windows, most fails would be at 16GB (or perhaps greater), without other changes.

    All of that being said, I have not searched for other ABCD candidates that might redeem the 'no + s0' version for 16-bit variable use.

    Evan, are you confident that the existing XORO32 32-bit output accessed as two 16-bit words cannot be used to feed the scrambler above, as Tony and I had hoped for creating a XOROACC extension?

    As far as Crush testing goes, all variants Tony and I discussed passed BigCrush.
  • evanhevanh Posts: 16,007
    edited 2020-01-13 01:05
    state[2] += rotl(s0 + s1, 9) + s0 + 1; // toggle + 1 on/off each iteration to increase equidistribution
    
    Ignoring the accumulation component, that fits in the the execution speed criteria of XORO32. Although, toggling the +1 would likely be dropped. That is effectively an extra state bit.
    Evan, are you confident that the existing XORO32 32-bit output accessed as two 16-bit words cannot be used to feed the scrambler above, as Tony and I had hoped for creating a XOROACC extension?
    I have no idea. The prop's first SIMD instruction. :) EDIT: Add a bunch of 2x16-bit SIMD's. Catch is they'd need to be single operand. The double operand instruction spaces are all but gone.

  • xoroshironotxoroshironot Posts: 309
    edited 2020-01-13 05:26
    evanh wrote: »
    state[2] += rotl(s0 + s1, 9) + s0 + 1; // toggle + 1 on/off each iteration to increase equidistribution
    
    Ignoring the accumulation component, that fits in the the execution speed criteria of XORO32. Although, toggling the +1 would likely be dropped. That is effectively an extra state bit.
    Evan, are you confident that the existing XORO32 32-bit output accessed as two 16-bit words cannot be used to feed the scrambler above, as Tony and I had hoped for creating a XOROACC extension?
    I have no idea. The prop's first SIMD instruction. :) EDIT: Add a bunch of 2x16-bit SIMD's. Catch is they'd need to be single operand. The double operand instruction spaces are all but gone.
    The extra state bit increases the period to 2^49-2^17 in order to achieve the 32-bit near-perfect equidistribution... a nice touch, but obviously only if possible.

    The SIMD on x86/x64 is part of the reason that this approach (without toggling the +1) might be extremely desirable, as a double sum is a good way to automatically invoke compilers to generate an LEA instruction (where the +1 is free) and the accumulator is a good way to defer the scrambler in the pipeline, since its result is not immediately required. To wit:
    Seba's benchmark with 'HARNESS_ADD' run on Intel Xeon W3690@3.8GHz
    xoroshiro128+  : 3.891 s, 6.169 GB/s, 0.771 words/ns, 1.297 ns/word
    xoroshiro128++ : 4.438 s, 5.408 GB/s, 0.676 words/ns, 1.479 ns/word
    xoroshiro128+a : 3.922 s, 6.120 GB/s, 0.765 words/ns, 1.307 ns/word
    xoroshiro128++a: 4.203 s, 5.710 GB/s, 0.714 words/ns, 1.401 ns/word
    
    Which, if properly vetted, might make the accumulator versions the fastest all-purpose (non-cryptographic) PRNGs available that also:
    1. Pass all known randomness tests.
    2. Have fixed long period that is an easy, mathematically provable, extension to the base xoroshiro engine properties.
    3. Achieve near-perfect 2 dimensional equidistribution.
    4. Produces small and easily comprehensible source code.
    5. Are non-trivially invertible.
    6. Fairly easily implementable in hardware.
    7. Have inherited Jump support using existing xoroshiro subroutine.
    8. Allows for extension to near-perfect 1 dimensional equidistribution at twice the output word size (by toggling the +1 on and off, at the expense of less speed / higher implementation complexity).
    9. Reduces sparse state (a.k.a. escape-from-zero) recovery time and occurrence substantially.

    The 64-bit variable 'a' versions benchmarked here (I'm using Tony's naming convention), that we have been discussing in 16-bit variable form, are using a slightly more CPU pipeline optimized version of the original code that I posted for the 32-bit variable version (which I just updated and linked back to).
    I will release the final code in a few months, if all goes well.

    Edit: Expanded documentation
  • evanhevanh Posts: 16,007
    edited 2020-01-13 10:15
    That's looking nice. Good even when up against algorithms that are using magic multiplies. Huh, I wonder if avoiding the use of multiplies allows tighter loop times. There's bound to be more execution stages needed for big multiplies. In a tight loop that'll have an impact even in fat x86 CPUs.

    Right, I've just been checking progress of the full set distribution run of Xoroshiro64plusacc, I've named it Xoroacc32+1 for now. It's about 2/3rd way through. After sorting, candidate [13 5 10 9] is sitting at position 380. That's well outside the normal group cut-off that I'd then run gridding on. It has sum-ranking of 1812.

    [6 2 3 7] is in the lead at the moment with sum-ranking of 118.

    SumRanking formula is: sumRank = pRank * 2 + zRank + nzRank
  • TonyB_TonyB_ Posts: 2,190
    edited 2020-01-14 10:55
    evanh wrote: »
    Right, I've just been checking progress of the full set distribution run of Xoroshiro64plusacc, I've named it Xoroacc32+1 for now. It's about 2/3rd way through. After sorting, candidate [13 5 10 9] is sitting at position 380. That's well outside the normal group cut-off that I'd then run gridding on. It has sum-ranking of 1812.

    [6 2 3 7] is in the lead at the moment with sum-ranking of 118.

    SumRanking formula is: sumRank = pRank * 2 + zRank + nzRank

    Using the name xoroshiro32++a for xoroshiro32++ plus accumulator for the time being at least, two new constants are needed: e for low word increment and f for high word left rotation. e must be odd for 2D (32-bit) equidistribution, but f > 0 is not vital. Chris, does f = 1 improve PractRand scores when lsb is bit 0?

    Evan, what exactly are you testing and for how many iterations? As defined in the previous paragraph, is it xoroshiro32++a [a,b,c,d,1,1] or xoroshiro32++a [a,b,c,d,1,0]? To achieve full equidistribution requires 2^32-1 * 65536 double-iterations, however at just * 1 [13,5,10,9] should be bettered by other [a,b,c,d] depending on initial state and acc.

    f = 1 is easy for hardware, but a little tricky in software. For f = 0, I think the code below is the shortest possible for a software emulation of the new XOROACC functionality (six extra instructions). Could you please test this to see if results are the same as in C?

    P.S.
    I prefer equally weighting, i.e. sumRank = pRank + zRank + nzRank.
    Note that xoroshiro32++a [a,b,c,d,e,f) has two a's and I'm open to other suggestions, e.g. xoroshiro32++acc [a,b,c,d,e,f].
    'XOROACC instruction pseudo-code:
    'acc_out[15:0]  := acc_in[31:16] + prn[15:0] + 1
    'acc_out[31:16] := acc_out[15:0] + prn[31:16]
    
    'P2 code equivalent
    'acc[31:16] = current acc word, acc[15:0] = previous acc word
    'N.B. rotl(prn[31:16],1) omitted for speed
    
    	xoro32	state			'iterate state twice, PRN in next S
    	mov	prn,0			'prn = XORO32 PRN
    	shr	acc,#16			'acc = {16'b0,acc[31:16]}
    	add	acc,prn			'acc = {16'bX,acc[31:16]+prn[15:0]}
    	add	acc,#1			'acc = {16'bX,acc[31:16]+prn[15:0]+1}
    	movbyts	acc,#%%1010		'acc = {acc[31:16]+prn[15:0]+1,acc[31:16]+prn[15:0]+1}
    	bitl	prn,#111100000		'prn = {prn[31:16],16'b0}
    	add	acc,prn			'acc = {acc[31:16]+prn[15:0]+1+prn[31:16],acc[31:16]+prn[15:0]+1}
    
    state	long	1
    prn	long	0
    acc	long	0
    
    'In future hardware, the above code could be replaced by
    
    	xoro32	state
    	xoroacc	acc,0			'preserve acc[31:16] for next xoroacc
    	
    state	long	1
    acc	long	0
    
  • xoroshironotxoroshironot Posts: 309
    edited 2020-01-13 19:17
    Evan, regarding your comments, what is even more surprising, is that doing the following 'magic' will still result in slightly faster execution than that of xoroshiro128++ or Lehmer:
    uint64_t s2 = state[2] * 1181783497276652981; // large prime multiplier
    uint64_t tmp = state[2] + rotl(s0 + s1, 31) + 1;
    
    I need something other than my current Xeons to confirm these results on a more modern architecture... still considering a modern Threadripper.
    The speed delta I am seeing between + and ++ matches the 15% that Seba reported, which provides some comfort.

    Thanks for looking at other constants... that you have found a 118 without using '+ s0' is very surprising to me!
    An 'ideal' sum the way you are doing it is equal to about 92, or 80 when using Tony's equal weighting, but could vary substantially from about 40 all the way up to about 160 on these OR a perfect PRNG, depending entirely on seed.

    Extra detail:
    The SumRanking formula I had been using in my code, as I recall, was something like: average(pranksum/N1 + zranksum/N2 + nzranksum/N3)
    N1, N2 and N3 are are either degrees of freedom (n-1) or the number of values added (n), resulting in a weighted average between 0.5 and 2.0.
    Since I was combining data from multiple randomly seeded runs, I was using 14, 24 and 51, respectively, since my column sums were decidely non-zero out to that range.
    The larger values are due a larger number of x in freq(x) since larger state and/or more random PRNGs will commonly bump into those values across multiple runs.
    Obviously none of this 'extra detail' is required for fast vetting of candidates.


  • TonyB_ wrote: »
    Using the name xoroshiro32++a for xoroshiro32++ plus accumulator for the time being at least, two new constants are needed: e for low word increment and f for high word left rotation. e must be non-zero for 2D (32-bit) equidistribution, but f > 0 is not vital. Chris, does f = 1 improve PractRand scores when lsb is bit 0?
    Tony, we found e must be odd, normally 1, with larger (usually prime, like AD65h) values only considered for improving multiple (more than 1) double-iterated stream freqs (which was a questionable pursuit, but not entirely without merit).
    I noticed f = 1 for rotation of the second iteration solved multiple issues simultaneously, primarily improving PractRand/freqs and enabling de-correlation of the second stream for double-iterated use.
    I later found that the rotation would suffice for the second iteration, in lieu of using a reverse-bits approach that would also (like rotl 1) superimpose the msb from the first iteration on the lsb for the second.
    Other non-zero f values near 1 might provide some minor benefit, but generally as f increases, we lose the benefit seen in PractRand, but I never looked a freqs in that case.
  • evanhevanh Posts: 16,007
    edited 2020-01-13 20:36
    Here's the automated distribution results. I had to append the .txt suffix for the forum attachment.

    EDIT: And a zip of the individual case outputs

    PS: Tony, these are all single iterated only. It doesn't look like I ever got the double iterating ported to new automation. Ah, found a long forgotten Bash variable to set to enable double iterating. I should make that something for the parameters I think.

    EDIT2: Also attached the source file for my algorithm implementation.
  • evanhevanh Posts: 16,007
    The number of iterations is still the non-accumulator full period. No automation changes were made to compensate for the accumulator.
  • evanhevanh Posts: 16,007
    edited 2020-01-13 20:51
    I've just modified the algorithm source. It doesn't affect the results, just reduces number of explicit masks applied when testing a non-native C word size.
  • xoroshironotxoroshironot Posts: 309
    edited 2020-01-14 01:04
    I looked at the TXT file and all of those look good to me if I am interpreting the data properly.

    Take a random example from half way down the list:
    Xoroacc32(16)+1 [2 11 3 8] 12 17 40 779 441 558 2557

    The pfChi is perfect... it should be about 12.
    The zfChi is only a little low... it should be about 20.
    The nzfChi is only a little low... it should be about 44.

    Edit: Each perfect value is approximately the total number of > 0 freqs summed, with each individual value that created the freqsum being about 1 (non-linear, but usually in the range 0.05 to 5).
  • evanhevanh Posts: 16,007
    edited 2020-01-13 22:58
    Gridding takes 40 mins I think per candidate. Err, it'll take longer with larger Practrand scores, maybe a few hours. I'll let you select what you want gridded.
  • xoroshironotxoroshironot Posts: 309
    edited 2020-01-14 20:27
    Since you are doing 'no + s0', in order to test a theory I would be interested in full grids on:
    [9 8 14 15]
    [13 9 8 7]
    [4 8 5 1] // should benefit from the same type of logic as rotl 1 being used on second iteration???

    I picked these from near the top of a chi square sum list I created using data from '(ideal freq sum - actual freq sum)^2/ideal freq sum' where ideal is the degrees of freedom I believe you are using.

    Edit: DF for zfChi is likely 20, not 22. Updated list.
    Sums less than 0.5 (which is kind of overly restrictive):
    Candidate	pfChi	zfChi	nzfChi	p	z	nz	s
    [11 2 6 11]	12	20	43	0.00	0.00	0.02	0.02
    [11 7 10 15]	12	19	44	0.00	0.05	0.00	0.05
    [12 10 13 2]	12	19	43	0.00	0.05	0.02	0.07
    [12 14 5 7]	12	19	43	0.00	0.05	0.02	0.07
    [14 8 9 11]	12	20	42	0.00	0.00	0.09	0.09
    [13 5 8 7]	11	19	44	0.08	0.05	0.00	0.13
    [9 8 14 15]	11	21	44	0.08	0.05	0.00	0.13
    [4 8 5 1]	12	21	42	0.00	0.05	0.09	0.14
    [2 11 3 0]	11	19	43	0.08	0.05	0.02	0.16
    [11 11 14 9]	11	19	45	0.08	0.05	0.02	0.16
    [7 10 10 4]	11	21	45	0.08	0.05	0.02	0.16
    [5 8 4 13]	13	19	45	0.08	0.05	0.02	0.16
    [14 3 13 11]	13	21	45	0.08	0.05	0.02	0.16
    [11 10 12 8]	13	20	42	0.08	0.00	0.09	0.17
    [13 5 10 2]	12	22	44	0.00	0.20	0.00	0.20
    [5 14 12 11]	12	20	41	0.00	0.00	0.20	0.20
    [11 1 14 3]	11	19	42	0.08	0.05	0.09	0.22
    [12 10 11 15]	11	18	44	0.08	0.20	0.00	0.28
    [5 14 12 7]	11	22	44	0.08	0.20	0.00	0.28
    [6 14 11 6]	11	20	41	0.08	0.00	0.20	0.29
    [4 8 5 13]	13	20	41	0.08	0.00	0.20	0.29
    [15 1 2 0]	13	20	47	0.08	0.00	0.20	0.29
    [2 1 15 12]	12	22	46	0.00	0.20	0.09	0.29
    [3 11 14 15]	11	18	45	0.08	0.20	0.02	0.31
    [3 11 14 10]	11	22	43	0.08	0.20	0.02	0.31
    [15 3 6 3]	13	18	45	0.08	0.20	0.02	0.31
    [12 10 11 7]	13	19	47	0.08	0.05	0.20	0.34
    [2 11 3 6]	10	20	43	0.33	0.00	0.02	0.36
    [4 1 9 14]	14	20	45	0.33	0.00	0.02	0.36
    [6 14 11 9]	11	18	42	0.08	0.20	0.09	0.37
    [6 2 11 5]	11	18	42	0.08	0.20	0.09	0.37
    [14 1 11 2]	11	18	46	0.08	0.20	0.09	0.37
    [11 11 14 4]	11	22	42	0.08	0.20	0.09	0.37
    [3 2 6 9]	13	18	46	0.08	0.20	0.09	0.37
    [13 3 14 6]	13	22	46	0.08	0.20	0.09	0.37
    [10 3 3 0]	10	19	45	0.33	0.05	0.02	0.41
    [10 3 3 3]	13	20	40	0.08	0.00	0.36	0.45
    [10 3 3 7]	14	21	42	0.33	0.05	0.09	0.47
    [1 2 8 5]	13	21	40	0.08	0.05	0.36	0.50
    
    This list might have a side effect of predicting that [13 5 10 2] might be about the best choice if '+ s0' were used.
    I am curious how [13 5 10 2] might have fared as the original XORO32.
    List was slightly incorrect.
  • evanhevanh Posts: 16,007
    edited 2020-01-15 00:20
    Here's gridding results of that whole list. Well, as it was 14 hours ago at least. And a quick summary of the grids. Wasn't as high scoring as I was expecting. EDIT: Added the individual Practrand reports as well.
    Scoring summary for single iterated sampling of Xoroacc32(16)+1
    
    Candidate	ExpMin	ExpAvg	ExpMax
    [11 2 6 11]	29	32.583	34
    [12 14 5 7]	28	32.416	34
    [13 3 14 6]	27	32.104	34
    [12 4 15 11]	27	32.062	34
    [5 14 12 7]	27	31.645	34
    [5 14 12 11]	27	31.562	34
    [11 11 14 4]	27	31.250	34
    [4 1 9 11]	27	30.854	33
    [6 14 11 6]	27	30.375	33
    [13 5 8 7]	26	31.520	34
    [4 8 5 13]	26	30.770	34
    [4 8 5 12]	26	30.583	34
    [13 9 8 7]	26	30.375	34
    [13 5 10 2]	25	31.708	34
    [4 8 5 1]	25	30.687	34
    [14 8 9 11]	25	30.166	34
    [14 8 9 8]	25	29.979	34
    [7 2 10 15]	24	30.833	33
    [9 8 14 15]	24	30.333	33
    [7 10 10 4]	24	29.354	33
    [10 3 3 7]	24	28.791	32
    [11 10 12 8]	23	31.541	34
    [12 10 13 2]	23	30.812	34
    [11 7 10 15]	23	28.604	33
    [8 2 1 9]	23	26.604	30
    [8 15 7 1]	20	22.166	25
    [6 2 5 2]	19	25.437	29
    [4 7 15 4]	18	27.583	33
    [1 2 8 5]	18	27.458	33
    [8 7 15 4]	18	21.604	27
    [14 3 13 11]	17	23.187	30
    [1 2 8 0]	17	22.854	33
    [3 11 14 10]	17	22.791	32
    [15 1 2 0]	17	20.020	27
    [2 1 15 6]	17	19.250	23
    [2 1 15 12]	17	19.250	23
    
  • evanhevanh Posts: 16,007
    edited 2020-01-15 00:41
    There's a problem with the C algorithm implementation. The scrambler is not being applied to the immediate result output. That means the first result is always bare state data.

    I was just setting up for a new run with +s0+1 when it sunk in. It really needs to return the accumulator state[2] directly, not s2. s2 becomes redundant.

    Here's updated algorithm:
    #include "xorotypes.h"
    
    static rngword_t   s[3] = { 1, 0, 0 }; // states 0 and 1 must not be seeded everywhere zero, state 2 should not to be used in isolation as a stream selector
    
    
    rngword_t nextrnd( void )
    {
    	rngword_t  s0 = s[0];
    	rngword_t  s1 = s[1];
    
    #if (ACCUM_SIZE == 16) || (ACCUM_SIZE == 32) || (ACCUM_SIZE == 64)
    	s[2] += rotl(s0 + s1, CONSTANT_D) + s0 + 1;
    	s1 ^= s0;
    	s[0] = rotl( s0, CONSTANT_A ) ^ s1 ^ (s1 << CONSTANT_B);
    #else
    	s[2] = (s[2] + rotl(s0 + s1, CONSTANT_D) + 1) & ACCUM_MASK;
    	s1 ^= s0;
    	s[0] = rotl( s0, CONSTANT_A ) ^ s1 ^ (s1 << CONSTANT_B) & ACCUM_MASK;
    #endif
    	s[1] = rotl(s1, CONSTANT_C);
    
    	return s[2];
    }
    
  • evanhevanh Posts: 16,007
    PS: I've renamed the just completed algorithm to Xoroshiro128+1a
    and this one is Xoroshiro128++1a

  • xoroshironotxoroshironot Posts: 309
    edited 2020-01-17 04:00
    Evan, thanks for that... it will take time to digest.

    Just in case I wasn't clear on the totals I'm using, look at this familiar table:
    #|Actual-Ideal|^2/Ideal																		
    #a	b	c	d	pfreq0	pfreq1	pfreq2	pfreq3	pfreq4	pfreq5	pfreq6	pfreq7	pfreq8	pfreq9	pfreq10	pfreq11	pfreq12	pfreq13	total
    #ideal				0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
    #scro				0	0	0	0	0	0	0	0	1	0	0	0	0	0	3
    
    The 'ideal' freq values of '0' are what cause a potential issue. If all of the freqs were '0' then the PRNG that produced them would be a perfectly uniform generator.
    However, 'uniform' <> 'random'. Truly random numbers are not completely uniform, otherwise you could predict their frequency distribution at some point in the period.
    This was great for xoroshiro32++, where it was known where to stop and look for uniformity... single and/or double iterated period.
    Therefore, that helped constants search, since it was a struggle to achieve full period randomness on freqs.
    However, now there will be an impact, since the entire period cannot be inspected directly.
    So now a new entry is required in the table:
    #a	b	c	d	pfreq0	pfreq1	pfreq2	pfreq3	pfreq4	pfreq5	pfreq6	pfreq7	pfreq8	pfreq9	pfreq10	pfreq11	pfreq12	pfreq13	total
    #normal				~1	~1	~1	~1	~1	~1	~1	~1	~1	~1	~1	~1	~1	~1	~14 (Actually closer to 12, see below).
    
    The above is an approximation, as are the following: as each '~1' is most likely to be in the range of 0.15 to 6, and the sum '~14' is likely to be in the range of 3 to 30.
    The actual value is normally one less to express Degrees of Freedom, but in this case it is approximately even one less than that due to the final point being so close to zero it causes skew... so effectively ~12 if Yates correction is applied.
    With only a single trial, high freq values near the end may carry excess weight, so stopping before actual '0's appear can minimize bias, which would eventually dissipate with the addition of more 'randomly seeded trials'* of the same PRNG. (see below)
    In that case you can sum the freqs and apply expected * #trials. For example (sum_of_pfreq0's - expected * count_of_pfreq0's)^2/(expected * count_of_pfreq0's).
    Normally, for long period PRNGs, you might just run your trial for a longer period and adjust your expected values accordingly.
    However, with these 'a' or 'acc' PRNGs, it is a known issue that each stream is correlated to the next, and subsequent streams are likely more strongly correlated with each other... and perhaps in freq results.
    *Thus the recommendation to not use consecutive streams for statistical evaluation.

    Edit: Clarify how degrees of freedom is calculated, and suggest why it is not the number of data points - 1 in all cases.
  • evanhevanh Posts: 16,007
    edited 2020-01-15 02:48
    Sorry, most of that is just going over my head. I do roughly get some of it but then the full talk just sinks me.

    I mostly followed Tony by asking for exact code that I could incorporate.

  • evanh wrote: »
    Wasn't as high scoring as I was expecting.
    Given what I just posted above, the odds of missing the best candidate, without testing all possible candidates, is very high.
    That being said, I am completely surprised with what [11 2 6 11] accomplished.

    This +1a design (as you call it) was intended to exceed xoroshiro+ in most all relevant considerations, which I believe it will do with most constants of sufficient linear complexity.
    The ++1a design (using '+ s0') was intended to do the same when compared with xoroshiro++, while offering the extended benefit of allowing the randomness period to be more predictable (e.g. 8GB for ++1a), thus compatible with a broader range of constants.
  • evanh wrote: »
    I mostly followed Tony by asking for exact code that I could incorporate.
    I did the same thing, if you look back through the posts over the last several months, and ended up porting most of your freq code after multiple discussions with Tony.
    Then I started looking at the freq data from high bits of xoroshiro128+ and **, which prompted me to agree with Tony that PractRand, Crush, etc. are insufficient for these evaluations.
    However, the status quo of how to interpret the freq data no longer applied, and after several false starts, it sunk in reasonably well.

    The biggest disappointment was finding that although the freqs for certain more complex variants would look good across multiple independent and consecutive streams, there was no getting around the fact that the streams are correlated, no matter how I tried to mask it... with the exception of stumbling onto rotl(xoroshiro32pp, 1) for use on the second iteration, which granted a solid doubling of randomness period (which I am still somewhat suspicious of, but cannot demonstrate any substantial correlation fault when [13, 5, 10, 9] was used).

  • evanhevanh Posts: 16,007
    What is a stream?
  • I was curious, so ran 16-bit aperture rotations on [11 2 6 11] +1a (no '+ s0') with the additional modification of applying rotl(, 1) to the second iteration:
          PractRand v0.94 options:  stdin -multithreaded -te 1 -tf 2 -tlmin 128KB
        |===00====01====02====03====04====05====06====07====08====09====10====11====12====13====14====15==
     16 |  16GB  32GB  32GB  32GB  32GB  32GB  32GB  32GB  32GB   8GB  16GB   8GB  16GB  16GB  16GB  32GB
    
    Much as expected.
  • TonyB_TonyB_ Posts: 2,190
    edited 2020-01-15 13:27
    evanh wrote: »
    What is a stream?
    The accumulator algorithm is existing XORO32 followed by new XOROACC. XORO32 on its own and XORO32 + XOROACC both generate a pseudo-32-bit PRN by concatenating two successive 16-bit ones. Let prn = XORO32 PRN and acc = new XOROACC PRN, then
    	prn = {prn1,prn0}
    	acc = {acc1,acc0}
    	
    	acc0 := acc1 + prn0 + e
    	acc1 := acc0 + prn1
    
    where e = new odd constant (and left rotation of prn1 by new constant f omitted for now for simplicity). Note that xoroshiro32++ state variable is not needed to calculate acc once prn is known.

    A stream is a unique sequence of 2^32-1 XOROACC iterations (i.e. 2^32-1 double iterations or 2^33-2 single iterations of xoroshiro32++) with a unique initial acc value. To be more precise, the high word of acc (acc1) is unique as it represents the current 16-bit output, whereas acc0 is the previous output. Thus there are 2^16 streams and the total XOROACC period is 2^48-2^16 double or 2^49-2^17 single iterations. This is the number of iterations needed for state and acc1 to repeat their initial values at the same time.

    Each stream is correlated with all the others, i.e. the difference between the outputs of any two streams for a particular iteration i is the same for all i. One stream leads automatically to the next one and there wil be perfect 1D (16-bit) equidistribution, including all zeroes, over the total period. However, for 2D (32-bit) equidistribution that is near perfect, the minimum requirement is that an odd value (e) must be added to either acc0 or acc1 just once for each stream output and e = 1 is the smallest and easiest value to use.

    To achieve 2D equidistribution, each stream must be output in full, with low and high words repeated but swapped over after 2^32-1 single iterations. However, the "-1" means acc will be one less than it should be at this midway point of the double iteration period and so adding 1 will correct this error. This + e adjustment could be done at any time during every other 2^32-1 single iteration period, but that would mean testing the state all the time. It is quicker to do the + e once every double iteration, which has the effect of interleaving all the streams by continually jumping from one to another. Each stream will be output fully over the total period.

    Evan, I'm not quite sure what you've been testing but I think you should be doing the following:

    Amend the pfreq/zfreq/nzfreq C program used to create the binary data files that showed [13,5,10,9] was the best overall candidate when combining PractRand scores and frequency distributions. You will need to add code for acc as shown above and use the 32-bit acc value as the index to the 4GB array instead of prn. There must be 2^32-1 double iterations.

    I still prefer equally weighting of pfreq/zfreq/nzfreq.

    Bear in mind that different initial acc values will affect how closely the distributions match the ideal ones, therefore more than one initial acc should be tested. Other [a,b,c,d] will be better than [13,5,10,9] for init_acc = X, but could [13,5,10,9] do just as well if not better when init_acc = Y?

    Do pfreq/zfreq/nzfreq ever repeat for different init_acc or is there always one best value? If the latter, do Chi-square scores fall away smoothly either side of this peak init_acc? Probably the only practicable way to test this is to look at xoroshiro16++ with all possible 256 acc1 values.
  • xoroshironotxoroshironot Posts: 309
    edited 2020-01-15 20:12
    Some clarification on equidistribution as it relates to the full double-iterated implementation:
    We know that XORO32 has near-perfect 1D 16-bit equidistribution, with two missing zeros (since two periods of underlying xoroshiro).
    XOROACC (in any discussed form) will always have perfect 1D 16-bit and near-perfect 2D 16-bit equidistribution.
    However, without removing the +1 from (A), and/or using rotl(,1) on (B), the second iteration, we cannot achieve the near-perfect 1D 32-bit equidistribution.
    That is because about half of the 2D 16-bit distribution is not aligned on same 32-bit boundary as the 32-bit output DWord.
    Doing A and/or B forces the 32-bit alignment, with both A and B together likely producing the highest degree of randomness.
    This is obviously a nice addition, since a true improvement on XORO32 would benefit from a concert of simultaneous improvements, if practical.
    Since the 32-bit 1D is not exactly perfect (65536 32-bit values occur once less often), a value judgement had to be made.
    Tony and I agreed in favor of, based on the same logic that allows for the missing 0s in XORO32...
    Not ideal, but close enough to be disregarded as being a problem, since it would be very difficult to detect that any given 32-bit value did not occur exactly the same number of times as any other value.
    That, as opposed to allowing only a complete normal 32-bit distribution (which was an incomplete normal with XORO32), where it is much easier to prove that all 32-bit values do/will not occur equally.

    Regarding +1:
    We know that +1 and ^1 are interchangeable (or any odd value, likely prime if larger values were used), and can be substituted with the same effects discussed above.
    The only requirement is that if xor is used, it be applied to PRN prior to addition with ACC.
    However, it is possible that the randomness will be affected, sometimes in a positive way, thus can be checked for comparison with +, if desired.
  • evanhevanh Posts: 16,007
    TonyB_ wrote: »
    Evan, I'm not quite sure what you've been testing but I think you should be doing the following:

    Amend the pfreq/zfreq/nzfreq C program used to create the binary data files that showed [13,5,10,9] was the best overall candidate when combining PractRand scores and frequency distributions. You will need to add code for acc as shown above and use the 32-bit acc value as the index to the 4GB array instead of prn. There must be 2^32-1 double iterations.

    Yep all that is automatic with each posted algorithm. I just plug it in now and the rest is sorted for me.

    I make a new algorithm file then name it matching a new entry in the config. Here's my plugging in config lines (bottom one is newest):
    REFERENCE	ENGINE TITLE		scrambler TITLE		COMMENT
    scro		Scro-rarns		
    xoo		Xoroshiro		+
    xoop		Xoroshiro		+p
    xo		Xoroshiro		++		scrambler = rotl( s1 + s0, CONSTANT_D ) + s0
    xop		Xoroshiro		+p+
    xogr		Xorograyro		++
    xon		Xoroshiron		++
    xofb		Xoroshirofb		++
    xom		Xoroshiro		+x1+
    xom0		Xoroshiro		+x0+		scrambler = rotl( s1 + s0, CONSTANT_D ) + s0 * 5
    xofbm		Xoroshirofb		+x+
    
    xo01		Xoroshiro		s1+s0-		scrambler = rotl( s1 + s0, CONSTANT_D ) - s0
    xo02		Xoroshiro		s1-s0+		scrambler = rotl( s1 - s0, CONSTANT_D ) + s0
    xo03		Xoroshiro		s1-s0-		scrambler = rotl( s1 - s0, CONSTANT_D ) - s0
    xo10		Xoroshiro		-s1+s0+		scrambler = rotl( -s1 + s0, CONSTANT_D ) + s0
    xo11		Xoroshiro		-s1+s0-		scrambler = rotl( -s1 + s0, CONSTANT_D ) - s0
    xo12		Xoroshiro		-s1-s0+		scrambler = rotl( -s1 - s0, CONSTANT_D ) + s0
    xo13		Xoroshiro		-s1-s0-		scrambler = rotl( -s1 - s0, CONSTANT_D ) - s0
    
    xom01		Xoroshiro		s1+s0-x 	scrambler = rotl( s1 + s0, CONSTANT_D ) - s0 * 5
    xom02		Xoroshiro		s1-s0+x 	scrambler = rotl( s1 - s0, CONSTANT_D ) + s0 * 5
    xom03		Xoroshiro		s1-s0-x 	scrambler = rotl( s1 - s0, CONSTANT_D ) - s0 * 5
    xom10		Xoroshiro		-s1+s0+x	scrambler = rotl( -s1 + s0, CONSTANT_D ) + s0 * 5
    xom11		Xoroshiro		-s1+s0-x	scrambler = rotl( -s1 + s0, CONSTANT_D ) - s0 * 5
    xom12		Xoroshiro		-s1-s0+x	scrambler = rotl( -s1 - s0, CONSTANT_D ) + s0 * 5
    xom13		Xoroshiro		-s1-s0-x	scrambler = rotl( -s1 - s0, CONSTANT_D ) - s0 * 5
    
    xoa1		Xoroshiro		+1a		scrambler += rotl( s1 + s0, CONSTANT_D ) + 1
    xoa00		Xoroshiro		++1a		scrambler += rotl( s1 + s0, CONSTANT_D ) + s0 + 1
    
  • evanhevanh Posts: 16,007
    So just completed full distribution run of ++1a. CSV file attached
Sign In or Register to comment.