Chapter 32 — Cycles and Remainders

# --- Visualization: Multiplicative order cycles as circular graphs ---
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches

plt.style.use('seaborn-v0_8-whitegrid')

def multiplicative_order(a, n):
    """Find smallest k > 0 such that a^k ≡ 1 (mod n). Returns None if gcd(a,n)!=1."""
    import math
    if math.gcd(a, n) != 1:
        return None
    power = a % n
    for k in range(1, n + 1):
        if power == 1:
            return k
        power = (power * a) % n
    return None

PRIME = 13   # <-- try other primes: 7, 11, 17, 19
bases = range(1, PRIME)

fig, axes = plt.subplots(3, 4, figsize=(14, 10))
axes = axes.flatten()

for idx, a in enumerate(bases):
    ax = axes[idx]
    order = multiplicative_order(a, PRIME)
    sequence = []
    val = 1
    for _ in range(order):
        val = (val * a) % PRIME
        sequence.append(val)

    # Place sequence elements on a circle
    n_pts = len(sequence)
    angles = [2 * np.pi * i / n_pts for i in range(n_pts)]
    xs = [np.cos(th) for th in angles]
    ys = [np.sin(th) for th in angles]

    color = 'crimson' if order == PRIME - 1 else 'steelblue'  # red = primitive root
    ax.scatter(xs, ys, s=60, c=color, zorder=3)
    for i in range(n_pts):
        ax.annotate(str(sequence[i]), (xs[i], ys[i]),
                    ha='center', va='center', fontsize=7)
        j = (i + 1) % n_pts
        ax.annotate('', xy=(xs[j], ys[j]), xytext=(xs[i], ys[i]),
                    arrowprops=dict(arrowstyle='->', color=color, lw=1))

    ax.set_xlim(-1.5, 1.5)
    ax.set_ylim(-1.5, 1.5)
    ax.set_aspect('equal')
    ax.axis('off')
    label = f'a={a}, ord={order}'
    if order == PRIME - 1:
        label += ' [gen]'
    ax.set_title(label, fontsize=9, color=color)

# Hide unused subplots
for idx in range(len(bases), len(axes)):
    axes[idx].axis('off')

fig.suptitle(f'Multiplicative cycles mod {PRIME}  (red = primitive roots / generators)',
             fontsize=13)
plt.tight_layout()
plt.show()

# --- Implementation: Euler's totient, primitive roots, discrete log ---
import math
import numpy as np

def euler_totient(n):
    """
    Compute Euler's totient φ(n): count of integers in [1,n] coprime to n.

    Args:
        n: positive integer

    Returns:
        φ(n)
    """
    result = n
    p = 2
    temp = n
    while p * p <= temp:
        if temp % p == 0:
            while temp % p == 0:
                temp //= p
            result -= result // p  # result *= (1 - 1/p)
        p += 1
    if temp > 1:
        result -= result // temp
    return result


def multiplicative_order(a, n):
    """
    Compute the multiplicative order of a modulo n.
    Returns None if gcd(a, n) != 1.

    Args:
        a: integer
        n: modulus

    Returns:
        smallest k > 0 with a^k ≡ 1 (mod n), or None
    """
    if math.gcd(a, n) != 1:
        return None
    order = 1
    power = a % n
    while power != 1:
        power = (power * a) % n
        order += 1
    return order


def find_primitive_roots(n):
    """
    Find all primitive roots modulo n.

    Args:
        n: modulus

    Returns:
        list of primitive roots
    """
    phi = euler_totient(n)
    roots = []
    for a in range(2, n):
        if math.gcd(a, n) == 1 and multiplicative_order(a, n) == phi:
            roots.append(a)
    return roots


def discrete_log_brute(g, b, n):
    """
    Brute-force discrete logarithm: find k such that g^k ≡ b (mod n).
    Only feasible for small n.

    Args:
        g: generator (primitive root)
        b: target value
        n: modulus

    Returns:
        k, or None if not found
    """
    power = 1
    for k in range(euler_totient(n)):
        if power == b % n:
            return k
        power = (power * g) % n
    return None


# --- Tests ---
for n in [7, 12, 13, 15, 17]:
    phi = euler_totient(n)
    roots = find_primitive_roots(n)
    print(f"n={n:2d}: φ(n)={phi:2d}, primitive roots: {roots if roots else 'none'}")

print("\n--- Discrete log examples (mod 7, generator g=3) ---")
g, p = 3, 7
for b in range(1, p):
    k = discrete_log_brute(g, b, p)
    print(f"  {g}^{k} ≡ {b} (mod {p})   →  log_{g}({b}) mod {p} = {k}")

# --- Experiment 1: Order distribution across elements ---
# Hypothesis: for prime p, every divisor of p-1 appears as some element's order.
# Try changing: PRIME to different values.

import matplotlib.pyplot as plt
from collections import Counter

PRIME = 19  # <-- modify: try 7, 11, 13, 17, 23

orders = [multiplicative_order(a, PRIME) for a in range(1, PRIME)]
order_counts = Counter(orders)

divisors_of_phi = [d for d in range(1, PRIME) if (PRIME - 1) % d == 0]

fig, ax = plt.subplots(figsize=(9, 5))
bars = ax.bar([str(o) for o in sorted(order_counts)],
              [order_counts[o] for o in sorted(order_counts)],
              color='steelblue', edgecolor='white')
ax.set_xlabel('Multiplicative order', fontsize=11)
ax.set_ylabel('Number of elements with this order', fontsize=11)
ax.set_title(f'Order distribution in ℤ/{PRIME}ℤ*  (φ({PRIME}) = {PRIME-1})', fontsize=12)
plt.tight_layout()
plt.show()

print(f"Divisors of φ({PRIME}) = {PRIME-1}: {divisors_of_phi}")
print(f"Orders that appear: {sorted(order_counts.keys())}")
print(f"φ(d) elements should have order d. Check: ", end="")
for d in divisors_of_phi:
    phi_d = euler_totient(d)
    print(f"order {d}: expected {phi_d}, got {order_counts.get(d, 0)}", end='  |  ')

# --- Experiment 2: Linear Congruential Generator cycle length ---
# Hypothesis: a LCG x_(n+1) = (a*x_n + c) mod m achieves full period m
#             iff Hull-Dobell conditions are met.
# Try changing: a, c, m to see full vs partial cycles.

def lcg_cycle(a, c, m, seed=1):
    """Generate LCG sequence until it repeats. Returns cycle length."""
    seen = {}
    x = seed % m
    for i in range(m + 1):
        if x in seen:
            return i - seen[x], i  # (cycle_length, steps_to_enter_cycle)
        seen[x] = i
        x = (a * x + c) % m
    return m, 0

# Hull-Dobell conditions for full period:
#   1. gcd(c, m) = 1
#   2. a ≡ 1 (mod p) for every prime p dividing m
#   3. if 4 | m, then a ≡ 1 (mod 4)

M = 16   # <-- modify
test_cases = [
    (5, 3, M),   # full period
    (3, 0, M),   # no addend — partial
    (2, 1, M),   # a not ≡ 1 mod 2? let's check
    (5, 0, M),   # c=0 breaks Hull-Dobell condition 1
]

print(f"LCG cycle analysis (m={M}):")
print(f"{'a':>3}  {'c':>3}  {'cycle len':>10}  {'full period?':>14}")
print("-" * 38)
for a, c, m in test_cases:
    clen, _ = lcg_cycle(a, c, m)
    print(f"{a:>3}  {c:>3}  {clen:>10}  {'YES' if clen == m else 'NO':>14}")

# --- Mini Project: LCG Random Number Generator Analysis ---
# Problem: Build an LCG, verify its period, test uniformity, reveal lattice structure.
# Task: complete TODOs and interpret the plots.

import numpy as np
import matplotlib.pyplot as plt

class LCG:
    """
    Linear Congruential Generator: x_{n+1} = (a * x_n + c) mod m
    """
    def __init__(self, a, c, m, seed=42):
        self.a = a
        self.c = c
        self.m = m
        self.state = seed % m

    def next(self):
        """Advance and return next value in [0, m-1]."""
        self.state = (self.a * self.state + self.c) % self.m
        return self.state

    def generate(self, n):
        """Generate n values, normalized to [0, 1)."""
        return np.array([self.next() / self.m for _ in range(n)])


# Parameters from Numerical Recipes (historically widely used)
NR_LCG   = LCG(a=1664525,  c=1013904223, m=2**32, seed=12345)
# Intentionally bad LCG for comparison
BAD_LCG  = LCG(a=65539,    c=0,          m=2**31, seed=12345)  # RANDU

N = 100_000
good_samples = NR_LCG.generate(N)
bad_samples  = BAD_LCG.generate(N)

fig, axes = plt.subplots(2, 2, figsize=(13, 10))

# Uniformity histogram — good LCG
axes[0, 0].hist(good_samples, bins=50, color='steelblue', edgecolor='white')
axes[0, 0].set_title('Good LCG: Distribution', fontsize=12)
axes[0, 0].set_xlabel('Value (normalized)', fontsize=10)
axes[0, 0].set_ylabel('Count', fontsize=10)
axes[0, 0].axhline(N / 50, color='red', linestyle='--', label='Expected uniform')
axes[0, 0].legend()

# Uniformity histogram — bad LCG
axes[0, 1].hist(bad_samples, bins=50, color='salmon', edgecolor='white')
axes[0, 1].set_title('RANDU LCG: Distribution', fontsize=12)
axes[0, 1].set_xlabel('Value (normalized)', fontsize=10)
axes[0, 1].set_ylabel('Count', fontsize=10)

# 2D scatter — good LCG: no visible structure
N_PLOT = 5000
axes[1, 0].scatter(good_samples[:N_PLOT], good_samples[1:N_PLOT+1],
                   s=1, alpha=0.3, color='steelblue')
axes[1, 0].set_title('Good LCG: Consecutive pairs (x_n, x_{n+1})', fontsize=12)
axes[1, 0].set_xlabel('x_n', fontsize=10)
axes[1, 0].set_ylabel('x_{n+1}', fontsize=10)

# 2D scatter — RANDU: visible lattice structure!
axes[1, 1].scatter(bad_samples[:N_PLOT], bad_samples[1:N_PLOT+1],
                   s=1, alpha=0.3, color='salmon')
axes[1, 1].set_title('RANDU: Consecutive pairs — lattice structure!', fontsize=12)
axes[1, 1].set_xlabel('x_n', fontsize=10)
axes[1, 1].set_ylabel('x_{n+1}', fontsize=10)

plt.tight_layout()
plt.show()

print("The RANDU LCG (a=65539, m=2^31) has all output triples (x_n, x_{n+1}, x_{n+2})")
print("lying on just 15 parallel planes in 3D — catastrophic for Monte Carlo methods.")