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

## No comments:

Post a Comment