from collections.abc import Generator
from typing import NamedTuple, Mapping

#define RN (((seed = seed*11109+13849) >> 16) & 0xffff)

a = 11109
c = 13849
m = 2**31

from math import gcd
from toy_crypto.nt import factor

n = 18
f_n = factor(n)
print(f'Factorization of {n} is {f_n}.')
print(f'Factorization as a list of (prime, exponent) tuples: {list(f_n)}')

# The radical_value of a factorization is the product of each prime factor
# each used only once.
print(f'Radical value of {n} is {f_n.radical_value()}')

Factorization of 18 is 2 * 3^2.
Factorization as a list of (prime, exponent) tuples: [(2, 1), (3, 2)]
Radical value of 18 is 6

f_m = factor(m)
print(f'factorization of "m" ({m}) is {f_m}')

factorization of "m" (2147483648) is 2^31

f_c = factor(c)
print(f'factorization of "c" ({c}) is {f_c}')

factorization of "c" (13849) is 11 * 1259

a_dim = a - 1
f_adim = factor(a_dim)
print(f'Factorization of {a} - 1 is {f_adim}.')

Factorization of 11109 - 1 is 2^2 * 2777.

class Conditions:
    descriptons: Mapping[str,str] = {
        'SIMPLE': "simple value conditions",
        'COPRIME': "rel prime",
        'EACH_PRIME': "each factor",
        'RULE4': "rule of 4",
        'NO_EXTRA': "no extra",
        'SQUARE_LADEN': "not square-free"
    }
    full_set: set[str] = set(descriptons.keys())

    def __init__(self) -> None:
        self.value = set()
        

    def add(self, c) -> None:
        if c not in self.full_set:
            raise ValueError(f'Unknown condition "{c}"')
        self.value |= {c}

    def contains(self, c) -> bool:
        return c in self.value

    @classmethod
    def all_conds(cls) -> set[str]:
        return cls.full_set

    def all(self) -> bool:
        return self.value == self.full_set

    def items(self) -> list[tuple[str, bool]]:
        r: list[tuple[str,bool]] = []
        for c in self.full_set:
            r.append((c, c in self.value))
        return r

# Hull-Dobell theorem conditions (and a few others)
hd_conditions = Conditions()

# simple conditions
if m > 0 and (c % m) != 0:
     hd_conditions.add('SIMPLE')
print("c and m "
     + ("do " if hd_conditions.contains('SIMPLE') else "do not ")
     +  "meet simple value conditions")

# m and c are relatively prime
if gcd(m, c) == 1:
     hd_conditions.add('COPRIME')
print("c and m "
    + ("are " if hd_conditions.contains('COPRIME') else "are not ")
    + "relatively prime")

c and m do meet simple value conditions
c and m are relatively prime

# $a - 1$ is divisible by a single instance of each prime factors of $m$
if (a - 1) % f_m.radical_value() == 0:
        hd_conditions.add('EACH_PRIME')
print("a - 1 "
        + ("is " if hd_conditions.contains('EACH_PRIME') else "is not " )
        + "divisible by all prime factors of m")

a - 1 is divisible by all prime factors of m

if  m % 4 != 0 or (a - 1) % 4 == 0:
        hd_conditions.add('RULE4')
print("a - 1 "
        + ("is " if hd_conditions.contains('RULE4') else "is not " )
        + "divisible by 4 if m is divisible by 4")

a - 1 is divisible by 4 if m is divisible by 4

# already computed above
# f_m = factor(m)
# f_adim = factor(a-1)

print(f'factorization of m: {f_m}')
print(f'factorization of a - 1: {f_adim}')

# We will produce a reduced a - 1 by dividing each by the unique
# prime factors of m. This operation depends on condition 2 being met
reduced_a = (a - 1) // f_m.radical_value()

if m % 4 == 0: # this step assumes that condition 3 is met
    reduced_a //= 2

f_reduced_a = factor(reduced_a)
print(f'reduced a factorization: {f_reduced_a}')

if gcd(reduced_a, m) == 1:

    hd_conditions.add('NO_EXTRA')
print('a - 1 and m '
        + ('do not ' if hd_conditions.contains('NO_EXTRA') else 'do ')
        + 'share "extra" factors')

factorization of m: 2^31
factorization of a - 1: 2^2 * 2777
reduced a factorization: 2777
a - 1 and m do not share "extra" factors

#hd_conditions[COND_SQUARE_LADEN] = not is_squarefree(m)
if a % m != 1:
    hd_conditions.add('SQUARE_LADEN')
print('m ' + ('is not ' if hd_conditions.contains('SQUARE_LADEN') else 'is ') + 'square free')

m is not square free

if hd_conditions.all():
    print('all conditions are met')

all conditions are met

def dh_lcg_check(m: int, a: int, c: int) -> Conditions:
    # We are going to check simple value conditions and treat
    # them as value errors instead of a condition that fails.
    for name, value in {"m": m, "a": a, "c": c}.items():
        if value == 0:
            raise ValueError(f'"{name}" cannot be zero')

    if m < 0:
        raise ValueError(f'a ({a}) must be positive')

    # If we knew that m would either be prime or a power of primes,
    # we could perform all of these tests more simply (though perhaps)
    # more abstractly) in GF(m). But we can't assume that. So these tests
    # are messy.

    f_m = factor(m) 
    
    conds = Conditions()
    if m > 0 and (c % m) != 0:
        conds.add('SIMPLE')
    if gcd(m, c) == 1:
        conds.add('COPRIME')
    if  a % m != 1:
        conds.add('SQUARE_LADEN')
    if (a - 1) % f_m.radical_value() == 0:
        conds.add('EACH_PRIME')
    if m % 4 != 0 or (a - 1) % 4 == 0:
        conds.add('RULE4')


    # We can only compute the "no extra common factors" test if 
    # if we have passed the "each factor" test and the "rule of 4" test.
    #
    # We will consider the test a pass if the conditions for the test
    # are not met
    if not ( conds.contains('EACH_PRIME') and  conds.contains('RULE4')):
        conds.add('NO_EXTRA')
    else:
        # We test for no extra factors by first dividing a - 1 by
        # the things needed for the other conditions. If that reduced a
        # shares no factors with m there are no extra common factors. 
        reduced_a: int = (a - 1) // f_m.radical_value()
        if m % 4 == 0: # this step assumes that condition 3 is met
            reduced_a //= 2
        if gcd(reduced_a, m) == 1:
            conds.add('NO_EXTRA')

    return conds

class LCGSetting(NamedTuple):
    m: int # modulus
    a: int # multipier
    c: int # increment
    d: int = 0 # discard bits
    u: int = 0 # upper bits, to discard

lcgs: dict[str, LCGSetting] = {}

lcgs["rogue"] = LCGSetting(m, a, c, 16, 0)
lcgs["Knuth & Lewis"] = LCGSetting(2**32, 1664525, 1013904223)
lcgs["ANSI C"] = LCGSetting(2**31, 1103515245, 12345, 16, 0)
# The Sinclair ZX81. The UK answer to the Commodore64
lcgs["ZX81 (Bad)"] = LCGSetting(2**16 + 1, 75, 74)
lcgs["VMS"] = LCGSetting(2**32, 69069, 1)
lcgs["SV32"] = LCGSetting(2**32, 0x915f77f5, 1)

# Some specifically bad ones
lcgs["Bad 4"] = LCGSetting(m, ((a-1)//2) +1, c) 
lcgs["Bad rel prime"] = LCGSetting(m , a, c * 2)
lcgs["Bad extra"] = LCGSetting(m, 1 + (a-1)*2, c)
lcgs["Bad each prime"] = LCGSetting(m*3, a, c)

# The only way I can think of to have something fail at "square free"
# while "each primne" is to make a = m + 1.
square_free_m = 2 * 3 * 5 * 7 * 11 * 13 * 17 * 19 * 23
lcgs["Bad square free"] = LCGSetting(square_free_m, square_free_m + 1, 29)

for (name, v) in lcgs.items():
    conditions = dh_lcg_check(v.m, v.a, v.c)
    if conditions.all():
        print(f'All conditions hold for {name}')
    else:
        fails = [k for (k, v) in conditions.items() if not v ]
        print(f'Conditions {fails} fail for {name}')

All conditions hold for rogue
All conditions hold for Knuth & Lewis
All conditions hold for ANSI C
Conditions ['EACH_PRIME'] fail for ZX81 (Bad)
All conditions hold for VMS
All conditions hold for SV32
Conditions ['RULE4'] fail for Bad 4
Conditions ['COPRIME'] fail for Bad rel prime
Conditions ['NO_EXTRA'] fail for Bad extra
Conditions ['EACH_PRIME'] fail for Bad each prime
Conditions ['SQUARE_LADEN'] fail for Bad square free

def unshifted_lcg(
                modulus: int,
                multiplier: int,
                increment: int,
                seed: int) -> Generator[int, None, None]:
    
    state = seed
    a = multiplier % modulus
    c = increment % modulus
    while True:
        state = (a * state + c) % modulus
        yield state

seed = 0xABAD5EED # well, not a particularly bad seed.
rogue_unshifted = unshifted_lcg(m, a, c, seed)

for _ in range(10):
    print(next(rogue_unshifted))

1515747482
19507419
1959566720
1932450201
1342751350
203341991
1920877820
1634189701
1534125714
172340147

for _ in range(5):
    r = next(rogue_unshifted)
    lsb = r & 1
    print(lsb)

0
1
0
1
0

for _ in range(16):
    r = next(rogue_unshifted)
    lsb3 = r & 6
    lsb3 = lsb3 >> 1
    print(f'bits 3 and 2: {lsb3:02b}')

bits 3 and 2: 10
bits 3 and 2: 01
bits 3 and 2: 01
bits 3 and 2: 00
bits 3 and 2: 00
bits 3 and 2: 11
bits 3 and 2: 11
bits 3 and 2: 10
bits 3 and 2: 10
bits 3 and 2: 01
bits 3 and 2: 01
bits 3 and 2: 00
bits 3 and 2: 00
bits 3 and 2: 11
bits 3 and 2: 11
bits 3 and 2: 10

def lcg(modulus: int,
        multiplier: int,
        increment: int,
        seed: int,
        discard_bits: int) -> Generator[int, bool, None]:

    state = seed
    m = modulus
    a = multiplier % m
    c = increment % m

    # Python has int.bit_length(), but Sage doesn't. Sage has Integer.nbits()
    if not discard_bits < m.bit_length():
        raise ValueError("can't discard so many bits")

    while True:
        state = (a * state + c) % m
        reseed = yield state >> discard_bits
        if reseed:
            state = seed

# And now we can constuct the (intended) rogue RNG
rogue_lcg = lcg(m, a, c, seed, 16)

for _ in range(10):
    print(f'{next(rogue_lcg):8}')

   23128
     297
   29900
   29486
   20488
    3102
   29310
   24935
   23408
    2629

int seed;
seed = /* initial seed set from time and PID */
#define RN (((seed = seed*11109+13849) >> 16) & 0xffff)

def acg(modulus: int,
        multiplier: int,
        increment: int,
        seed: int,
        discard_bits: int,
        discard_high) -> Generator[int, bool, None]:

    state = seed
    m = modulus
    a = multiplier % m
    c = increment % m

    # Python has int.bit_length(), but Sage doesn't. Sage has Integer.nbits()
    if not discard_bits < m.bit_length():
        raise ValueError("can't discard so many bits")

    if discard_high > discard_bits:
        raise ValueError("Can't discard more high bits than all discarded")
    
    high_bit_modulus = m - (2**discard_high)

    while True:
        state = (a * state + c) % m
        reseed = yield (state % high_bit_modulus) >> (discard_bits - discard_high)
        if reseed:
            state = seed

sv32_settings = lcgs["SV32"]
sv32h_acg = acg(sv32_settings[0], sv32_settings[1], sv32_settings[2], seed=1, discard_bits=8, discard_high=6)

print([next(sv32h_acg) for _ in range(10)])

[807164125, 527615617, 307547762, 773737748, 340240924, 890032332, 643728697, 708426375, 454352427, 724802599, 63664208, 960981130, 623692810, 141421458, 725167799, 940147229, 794794681, 991040013, 584724334, 928858688]

Rogue background¶

The PRNG¶

Sageless¶

m and c relatively prime¶

$a - 1$ divisible by all unique prime factors of m¶

Rule of 4¶

Simple conditions on m and c¶

No extra common factors¶

a is not 1 mod m¶

checking more generally¶

Dealing with low order bits¶

The rogue LCG modulus¶

Unsigned overflow¶

Unsigned bit shift¶

`sizeof(int)`¶

15 bit multiplier and increment¶

Discarding some high bits¶

Rogue background¶

The PRNG¶

Sageless¶

m and c relatively prime¶

$a - 1$ divisible by all unique prime factors of m¶

Rule of 4¶

Simple conditions on m and c¶

No extra common factors¶

a is not 1 mod m¶

checking more generally¶

Dealing with low order bits¶

The rogue LCG modulus¶

Unsigned overflow¶

Unsigned bit shift¶

sizeof(int)¶

15 bit multiplier and increment¶

Discarding some high bits¶

`sizeof(int)`¶