Integers

"""Why is a negative times a negative positive? Three kinds of evidence.

1. Experimental: check (-a)(-b) = ab on a grid of integers.
2. Symbolic: let SymPy manipulate the letters a and b themselves.
3. Explanation: the lesson's proof from the distributive law, whose steps we
   display (SymPy cannot tell you *why* - it already has the rules built in).

Runs in the browser (Pyodide) with the ``sympy`` package loaded.
"""

import textwrap

import sympy

def grid(radius):
    """Return all integer pairs (a, b) with -radius <= a, b <= radius."""
    if radius < 0:
        raise ValueError("radius must be non-negative")
    values = range(-radius, radius + 1)
    return [(a, b) for a in values for b in values]

def check_on_grid(rule, pairs):
    """Apply ``rule(a, b) -> bool`` to every pair; return (cases, list of failing pairs)."""
    failures = [(a, b) for a, b in pairs if not rule(a, b)]
    return len(pairs), failures

def neg_times_neg(a, b):
    """The rule (-a)(-b) = ab for one pair of numbers."""
    return (-a) * (-b) == a * b

def neg_times_pos(a, b):
    """The rule (-a)b = -(ab) for one pair of numbers."""
    return (-a) * b == -(a * b)

def symbolic_difference(lhs, rhs):
    """Return ``sympy.expand(lhs - rhs)``; it is 0 exactly when SymPy proves lhs = rhs as polynomials."""
    return sympy.expand(lhs - rhs)

def symbolic_examples():
    """Return a list of (description, SymPy result) rows used in the printout."""
    a, b = sympy.symbols("a b")
    A, B = sympy.symbols("A B", commutative=False)
    return [
        ("(-a)*(-b)", (-a) * (-b)),
        ("(-a)*(-b) - a*b, expanded", symbolic_difference((-a) * (-b), a * b)),
        ("(-a)*b + a*b, expanded", symbolic_difference((-a) * b, -(a * b))),
        ("(a - b)*(a + b), expanded", sympy.expand((a - b) * (a + b))),
        ("how SymPy stores -a", sympy.srepr(-a)),
        ("(-A)*(-B), non-commuting A, B", (-A) * (-B)),
        ("A*B - B*A, non-commuting A, B", sympy.expand(A * B - B * A)),
    ]

#: Steps of the proof in the lesson. Each uses one ring law.
PROOF_STEPS = [
    ("0 = (-a) * 0", "x * 0 = 0, itself proved from distributivity"),
    ("  = (-a) * (b + (-b))", "b + (-b) = 0"),
    ("  = (-a)b + (-a)(-b)", "distributive law"),
    ("  = -(ab) + (-a)(-b)", "(-a)b = -(ab), proved the same way"),
    ("so (-a)(-b) = ab", "add ab to both sides"),
]

def main():
    print(f"1. Experiment: exact integer arithmetic on a grid (SymPy {sympy.__version__} used below)")
    pairs = grid(10)
    for label, rule in [("(-a)(-b) = ab   ", neg_times_neg), ("(-a)b  = -(ab)  ", neg_times_pos)]:
        cases, failures = check_on_grid(rule, pairs)
        print(f"   {label} for all a, b in -10..10: {cases} cases, {len(failures)} failures")
    print("   A grid can never reach every integer: this is evidence, not proof.")
    print()
    print("2. Symbolic algebra: SymPy works with the letters themselves")
    for description, result in symbolic_examples():
        print(f"   {description:<34} -> {result}")
    print(textwrap.fill(
        "SymPy turns (-a)*(-b) into a*b the moment it is typed: it stores -a as "
        "(-1)*a and multiplies the coefficients (-1)*(-1) = 1. A zero difference "
        "means the identity holds for all values in any commutative ring - "
        "provided SymPy's built-in rules are implemented correctly. The last two "
        "rows show the sign rule does not even need commutativity.",
        width=78, initial_indent="   ", subsequent_indent="   ",
    ))
    print()
    print("3. Why it is true (the lesson's proof, using only the ring laws):")
    for step, reason in PROOF_STEPS:
        print(f"   {step:<24} [{reason}]")
    print()
    print(textwrap.fill(
        "Epistemic status: part 1 is experimental evidence for 441 cases; part 2 "
        "is a symbolic verification that is only as trustworthy as SymPy's "
        "implementation of the ring axioms; only the proof in part 3 explains "
        "why (-a)(-b) = ab in every ring.",
        width=78,
    ))

if __name__ == "__main__":
    main()

"""Euler's prime-generating polynomial: forty successes, then a failure.

p(n) = n^2 + n + 41 is prime for n = 0, 1, 2, ..., 39 - forty cases in a row.
It is tempting to conclude "always prime". n = 40 breaks the pattern:
p(40) = 1681 = 41 * 41. Evidence, however plentiful, is not proof.

Runs unchanged in the browser (Pyodide) using only the standard library.
"""

import textwrap

def euler(n):
    """Euler's polynomial n^2 + n + 41."""
    return n * n + n + 41

def is_prime(n):
    """Exact primality test by trial division (fine for the small numbers used here)."""
    if n < 2:
        return False
    if n % 2 == 0:
        return n == 2
    d = 3
    while d * d <= n:
        if n % d == 0:
            return False
        d += 2
    return True

def smallest_factor(n):
    """The smallest prime factor of an integer n >= 2."""
    if n < 2:
        raise ValueError("n must be at least 2")
    if n % 2 == 0:
        return 2
    d = 3
    while d * d <= n:
        if n % d == 0:
            return d
        d += 2
    return n

def first_failure(poly, start=0, stop=10_000):
    """Return the first n in start..stop-1 where poly(n) is not prime, or None."""
    for n in range(start, stop):
        if not is_prime(poly(n)):
            return n
    return None

def prime_count(poly, stop):
    """How many of poly(0), ..., poly(stop - 1) are prime."""
    return sum(is_prime(poly(n)) for n in range(stop))

def main():
    print("Euler's polynomial p(n) = n^2 + n + 41. Is p(n) always prime?")
    print()
    values = [euler(n) for n in range(40)]
    for row in range(0, 40, 8):
        cells = " ".join(f"{v:>5}" for v in values[row:row + 8])
        print(f"  n = {row:>2}..{row + 7:<2}: {cells}")
    streak = sum(is_prime(v) for v in values)
    print(f"All {streak} of these 40 values are prime (checked exactly by trial division).")
    print()
    n = first_failure(euler)
    v = euler(n)
    f = smallest_factor(v)
    print(f"First failure: n = {n}, p({n}) = {v} = {f} * {v // f}.")
    print(f"Why n = 40 had to fail: p(40) = 40*40 + 40 + 41 = 40*41 + 41 = 41*41.")
    print(f"And p(41) = 41*41 + 41 + 41 = 41 * 43 = {euler(41)} fails too.")
    print()
    for stop in (100, 1000):
        print(f"Primes among p(0), ..., p({stop - 1}): {prime_count(euler, stop)} of {stop}")
    print("Still remarkably many - but a pattern that is 'usually true' is not a theorem.")
    print()
    print(textwrap.fill(
        "Epistemic status: forty confirming cases were experimental evidence, "
        "and they were not a proof - the single counterexample at n = 40 proves "
        "that n^2+n+41 is not always prime. (No non-constant polynomial with "
        "integer coefficients gives primes for every n; that is a theorem, not "
        "something this program checked.)",
        width=78,
    ))

if __name__ == "__main__":
    main()

"""Searching for inverses among the integers, and among the rationals.

For each integer n we search a window of integers m for an *additive inverse*
(n + m = 0) and a *multiplicative inverse* (n * m = 1). Additive inverses are
always found; multiplicative ones only for n = 1 and n = -1. Then we move to
the rational numbers, using ``fractions.Fraction`` for exact arithmetic, where
every non-zero number has a multiplicative inverse.

Runs unchanged in the browser (Pyodide) using only the standard library.
"""

from fractions import Fraction
import textwrap

def window(radius):
    """Return the integers -radius..radius (inclusive) as a range."""
    if radius < 0:
        raise ValueError("radius must be non-negative")
    return range(-radius, radius + 1)

def find_additive_inverse(n, candidates):
    """Return the first m in ``candidates`` with n + m == 0, or None if none is found."""
    for m in candidates:
        if n + m == 0:
            return m
    return None

def find_multiplicative_inverse(n, candidates):
    """Return the first m in ``candidates`` with n * m == 1, or None if none is found.

    ``None`` only means "not found among the candidates": on its own a search
    cannot prove that no inverse exists further away.
    """
    for m in candidates:
        if n * m == 1:
            return m
    return None

def rational_inverse(q):
    """Return the exact multiplicative inverse 1/q of a rational q, or None when q == 0."""
    q = Fraction(q)
    if q == 0:
        return None
    return 1 / q

def search_summary(n_radius, m_radius):
    """Search every n in -n_radius..n_radius over candidates m in -m_radius..m_radius.

    Returns ``(additive_found, multiplicative_found)``: the lists of n for which
    an additive / multiplicative inverse was found.
    """
    candidates = window(m_radius)
    additive = [n for n in window(n_radius) if find_additive_inverse(n, candidates) is not None]
    multiplicative = [n for n in window(n_radius) if find_multiplicative_inverse(n, candidates) is not None]
    return additive, multiplicative

def main():
    candidates = window(50)
    print("For each n, search m in -50..50 for n + m = 0 and for n * m = 1.")
    print("Then, in the rationals Q, compute 1/n exactly with fractions.Fraction.")
    print()
    print(f"{'n':>4} | {'n + m = 0':>9} | {'n * m = 1':>9} | {'1/n in Q':>8} | n * (1/n)")
    print("-" * 52)
    for n in window(6):
        add = find_additive_inverse(n, candidates)
        mul = find_multiplicative_inverse(n, candidates)
        inv = rational_inverse(n)
        mul_s = "none" if mul is None else str(mul)
        inv_s = "none" if inv is None else str(inv)
        check = "-" if inv is None else str(n * inv)
        print(f"{n:>4} | {add:>9} | {mul_s:>9} | {inv_s:>8} | {check}")
    print()
    radius = 30
    additive, multiplicative = search_summary(radius, radius)
    total = 2 * radius + 1
    print(f"Wider search, n and m both in -{radius}..{radius}:")
    print(f"  additive inverse found for {len(additive)} of {total} integers")
    print(f"  multiplicative inverse found for {len(multiplicative)} of {total}: n = {multiplicative}")
    q = Fraction(-7, 12)
    print(f"  in Q, e.g. q = {q}: 1/q = {rational_inverse(q)} and q * (1/q) = {q * rational_inverse(q)} exactly")
    print()
    print("Why the integer search can stop: if |n| >= 2 and m != 0 then |n * m| >= 2,")
    print("and if m = 0 then n * m = 0. So n * m = 1 is impossible however far we look.")
    print()
    print(textwrap.fill(
        "Epistemic status: the searches are experiments over finite windows (a "
        "failed search alone proves nothing about numbers outside the window); "
        "the rational arithmetic is exact. The two-line size argument above, "
        "not the search, proves that only 1 and -1 have integer inverses.",
        width=78,
    ))

if __name__ == "__main__":
    main()