For the case of CMWCs, it's beneficial for us to use use w-1. As noted by Marsaglia, this is because we can have w-1 as a primitive root. Briefly, with g as w-1 and n as a*g+1, g^{phi(n)/f} != 1 for all factors of phi(n). This is made easier when n is prime, because phi(n)=n-1. When n is composite, more work is required for establising the order: we need to find the least common multiple across the order of all the factors.
With w-1=65535, let's consider an a of 65518 (as used in [1]), and then arbitrarily its neighbour 65517.
# 65518*65535+1
65535^{4293722130/2} mod 4293722131 == 4293722130
65535^{4293722130/3} mod 4293722131 == 1070428133
65535^{4293722130/5} mod 4293722131 == 3696002153
65535^{4293722130/17} mod 4293722131 == 2477700139
65535^{4293722130/41} mod 4293722131 == 2876866897
65535^{4293722130/47} mod 4293722131 == 1432659653
65535^{4293722130/257} mod 4293722131 == 1579676435
... which indeed has a period of 4293722130, and we needn't look further because none of the results are 1. However, for the composite 65517*65535+1 we need to perform an exhaustive search.
# 65517*65535+1
65535^{1036399840/2} mod 4293656596 == 1
65535^{1036399840/5} mod 4293656596 == 1
65535^{1036399840/7} mod 4293656596 == 3849485225
65535^{1036399840/19} mod 4293656596 == 318890497
65535^{1036399840/113} mod 4293656596 == 3330338425
65535^{1036399840/431} mod 4293656596 == 2842235365
Walking the generator gives the period `12954998`. Now, to prove it we need to do something like:
n = 4293656596
# 4293656596: 2 2 29 37014281
# ...
# 1:
# 1:
# 28: 2 2 7
# 37014280: 2 2 2 5 19 113 431
# product of 28 and 37014280's factors
phi_n = 1036399840
for f in [2,5,7,19,113,431]:
print '# factor',f
for x in xrange(1,phi_n/f):
if pow(65535,f*x,n) == 1:
print(f, x)
break
# (2, 6477499)
# (5, 12954998)
# (7, 1850714)
# (19, 681842)
# (113, 114646)
# (431, 30058)
>>> reduce(lcm,(6477499,12954998,1850714,681842,114646,30058))
12954998
Bingo. Now, returning to Knuth's quote. There is indeed some more work to be done with w+1 when we consider c != 0, but not a lot more than handling w-1. The catch with using w+1, is that we'll always lose at least a factor of 2 for our order because of word sizes being 2^n. So, adjusting for w+1 doesn't actually give us a CMWC of longer period than MWC, but CMWCs do have other nice properties such as always being a multiple of the base for instance. While CMWCs in base 2^n will at least lose a factor of 2^4, 2^n+1 can optimally lose 2^1.
With 2^n+1=65537 here, we can use the following code to generate the expansion 65514*65537+1 with a period of `2146795509`. We handle the case of a real zero by stepping to the next state, ensuring the 0=256 mapping.
rand proc near
push dx
mov dx,65514
db 0xb8 ; mov ax,imm16
x dw 1
mul dx
db 0x05 ; add ax,imm16
c dw 1
adc dx,0
; handle wrap
sub ax,dx
ja done
jnz wrap
; skip an actual 0
xchg ax,dx
mov dx,65514
sub ax,dx
wrap: inc ax
dec dx
done: neg ax ; (b-1-x) == 0-x == neg ax
mov cs:[x],ax
mov cs:[c],dx
pop dx
ret
rand end
I wrote it for density, not performance. Modern processors will invalidate their I$ if modified, resulting in a performance hit. Changing the self-modifying code to point x and c into a data section is left as an exercise for the reader.
[1] - https://b2d-f9r.blogspot.com/2010/10/lcg-mod-65537.html
[2] - https://b2d-f9r.blogspot.com/2010/11/multiply-with-carry-generator-in-16-bit.html
