Experimental mathematics: conjectures and counterexamples

"""Rediscovering the Bailey-Borwein-Plouffe formula, and using it.

The constants S_j = sum_k 16^-k / (8k + j), j = 1..8, and pi are computed to
60 digits with exact rational arithmetic (pi by Machin's formula). The PSLQ
integer-relation algorithm (Ferguson and Bailey), implemented here with
Python's decimal module, finds the relation pi = 4 S_1 - 2 S_4 - S_5 - S_6,
which is then re-checked at 100 digits. Two cautionary runs follow: with only
8 digits PSLQ "finds" a spurious relation for (pi, e, 1), and at 40 digits it
only bounds the size of any relation. Finally the formula is used to compute
hexadecimal digits of pi far after the point without the earlier ones, and
these are compared with an independent exact computation. A found relation is
numerical evidence; the lesson gives the proof.
"""

from decimal import Decimal, localcontext, ROUND_HALF_EVEN
from fractions import Fraction
import math

HEX = "0123456789ABCDEF"

def arctan_inv(x, scale):
    """floor-ish of arctan(1/x) * scale with integer arithmetic (error < terms)."""
    total = term = scale // x
    x2, k, sign = x * x, 1, -1
    while term:
        term //= x2
        total += sign * (term // (2 * k + 1))
        sign, k = -sign, k + 1
    return total

def pi_fraction(digits):
    """pi to about `digits` decimal digits (Machin: pi = 16 atan 1/5 - 4 atan 1/239)."""
    scale = 10 ** (digits + 10)
    return Fraction(4 * (4 * arctan_inv(5, scale) - arctan_inv(239, scale)), scale)

def e_fraction(digits):
    """e to well beyond `digits` digits, as an exact partial sum of 1/k!."""
    terms = 10
    while math.lgamma(terms + 1) / math.log(10) < digits + 5:
        terms += 5
    return sum(Fraction(1, math.factorial(k)) for k in range(terms))

def bbp_sum(j, digits):
    """S_j = sum_k 16^-k / (8k + j), truncated with error below 10^-digits."""
    terms = int(digits / 1.2) + 5
    return sum(Fraction(1, 16 ** k * (8 * k + j)) for k in range(terms))

def pslq(xs, digits, max_norm, max_iter=1000):
    """PSLQ on the reals xs (given as Fractions) at `digits` digits.

    Returns (relation or None, iterations, bound): if a relation is found it is
    an integer vector; otherwise every relation has norm at least `bound`.
    """
    n = len(xs)
    with localcontext() as ctx:
        ctx.prec = digits + 10
        x = [Decimal(f.numerator) / Decimal(f.denominator) for f in xs]
        eps = Decimal(10) ** (2 - digits)
        gamma = (Decimal(4) / 3).sqrt()
        A = [[int(i == j) for j in range(n)] for i in range(n)]
        B = [[int(i == j) for j in range(n)] for i in range(n)]
        s = [sum((v * v for v in x[k:]), Decimal(0)).sqrt() for k in range(n)]
        y = [v / s[0] for v in x]
        s = [v / s[0] for v in s]
        H = [[Decimal(0)] * (n - 1) for _ in range(n)]
        for i in range(n):
            if i < n - 1:
                H[i][i] = s[i + 1] / s[i]
            for j in range(min(i, n - 1)):
                H[i][j] = -y[i] * y[j] / (s[j] * s[j + 1])

        def reduce(i, j):
            t = int((H[i][j] / H[j][j]).to_integral_value(ROUND_HALF_EVEN))
            if t:
                y[j] += t * y[i]
                for k in range(j + 1):
                    H[i][k] -= t * H[j][k]
                for k in range(n):
                    A[i][k] -= t * A[j][k]
                    B[k][j] += t * B[k][i]

        for i in range(1, n):
            for j in range(i - 1, -1, -1):
                reduce(i, j)
        bound = Decimal(0)
        for it in range(1, max_iter + 1):
            m = max(range(n - 1), key=lambda i: gamma ** (i + 1) * abs(H[i][i]))
            y[m], y[m + 1] = y[m + 1], y[m]
            A[m], A[m + 1] = A[m + 1], A[m]
            H[m], H[m + 1] = H[m + 1], H[m]
            for row in B:
                row[m], row[m + 1] = row[m + 1], row[m]
            if m < n - 2:
                t0 = (H[m][m] ** 2 + H[m][m + 1] ** 2).sqrt()
                t1, t2 = H[m][m] / t0, H[m][m + 1] / t0
                for i in range(m, n):
                    t3, t4 = H[i][m], H[i][m + 1]
                    H[i][m] = t1 * t3 + t2 * t4
                    H[i][m + 1] = -t2 * t3 + t1 * t4
            for i in range(m + 1, n):
                for j in range(min(i - 1, m + 1), -1, -1):
                    reduce(i, j)
            jmin = min(range(n), key=lambda j: abs(y[j]))
            if abs(y[jmin]) < eps:
                rel = [B[k][jmin] for k in range(n)]
                sign = 1 if next(c for c in rel if c) > 0 else -1
                return [sign * c for c in rel], it, None
            bound = 1 / max(abs(H[i][i]) for i in range(n - 1))
            if bound > max_norm:
                break
        return None, it, float(bound)

def format_combination(rel):
    """Write -(rel_1 S_1 + ... + rel_8 S_8) as text, e.g. '4 S_1 - 2 S_4'."""
    text = ""
    for k, c in enumerate(rel):
        coeff = -c
        if k == 0 or coeff == 0:
            continue
        sign = "-" if coeff < 0 else "+"
        size = "" if abs(coeff) == 1 else f"{abs(coeff)} "
        text += f" {sign} {size}S_{k}"
    text = text.strip()
    return text[2:] if text.startswith("+ ") else text

def residual(rel, xs):
    """|sum rel_i * x_i| computed exactly from the Fractions xs."""
    return abs(sum(c * v for c, v in zip(rel, xs)))

def bbp_hex_digits(d, count):
    """Hex digits of pi at positions d+1, ..., d+count (BBP digit extraction)."""
    def frac_series(j):
        s = 0.0
        for k in range(d + 1):
            s = (s + pow(16, d - k, 8 * k + j) / (8 * k + j)) % 1.0
        k = d + 1
        while True:
            term = 16.0 ** (d - k) / (8 * k + j)
            if term < 1e-17:
                return s
            s, k = s + term, k + 1

    x = (4 * frac_series(1) - 2 * frac_series(4)
         - frac_series(5) - frac_series(6)) % 1.0
    out = ""
    for _ in range(count):
        x *= 16
        digit = int(x)
        out += HEX[digit]
        x -= digit
    return out

def exact_hex_digits(d, count):
    """The same digits from an exact integer computation (Machin, base 16)."""
    guard = 8
    scale = 16 ** (d + count + guard)
    pi_scaled = 4 * (4 * arctan_inv(5, scale) - arctan_inv(239, scale))
    block = (pi_scaled // 16 ** guard) % 16 ** count
    return format(block, "X").rjust(count, "0")

def main():
    digits = 60
    xs = [pi_fraction(digits)] + [bbp_sum(j, digits) for j in range(1, 9)]
    rel, its, _ = pslq(xs, digits, max_norm=10 ** 12)
    print("PSLQ on (pi, S_1, ..., S_8) at 60 digits, "
          "S_j = sum_k 16^-k/(8k+j):")
    print(f"  relation after {its} iterations: {tuple(rel)}")
    print(f"  i.e. pi = {format_combination(rel)}")
    xs100 = [pi_fraction(100)] + [bbp_sum(j, 100) for j in range(1, 9)]
    print("  residual with 100-digit values: "
          f"{float(residual(rel, xs100)):.1e}")
    print()
    low = [pi_fraction(8), e_fraction(8), Fraction(1)]
    fake, _, _ = pslq(low, 8, max_norm=10 ** 6)
    high = [pi_fraction(40), e_fraction(40), Fraction(1)]
    print(f"(pi, e, 1) at 8 digits: 'relation' {tuple(fake)}, but at 40 digits")
    print(f"  its residual is {float(residual(fake, high)):.1e}: spurious.")
    none, its, bound = pslq(high, 40, max_norm=10 ** 10)
    print(f"(pi, e, 1) at 40 digits: no relation; after {its} iterations every")
    print(f"  integer relation would need norm > {bound:.1e}.")
    print()
    d = 5000
    fast, exact = bbp_hex_digits(d, 8), exact_hex_digits(d, 8)
    print(f"Hex digits of pi after the point: {bbp_hex_digits(0, 8)} ...")
    print(f"Positions {d + 1}-{d + 8} by BBP digit extraction: {fast}")
    print(f"Same positions by exact integer arithmetic:      {exact}")
    print()
    print("Epistemic status: the relation is numerical evidence (60 digits, with")
    print("the residual rechecked at 100); its proof is in the lesson. The 'no")
    print("relation' bound is PSLQ's guarantee only if 40 digits were enough;")
    print("whether pi + e is irrational is an open problem. The hex digits agree")
    print("between two independent methods: an exact check of those 8 digits.")

if __name__ == "__main__":
    main()

"""Patterns that fail: four searches with stated bounds.

1. Euler's sum-of-powers conjecture for fifth powers: an exhaustive search
   over all a <= b <= c <= d < e <= 150 (meet in the middle: sums of two fifth
   powers, looked up with NumPy) finds Lander and Parkin's counterexample.
2. The prime race modulo 4: primes 3 mod 4 lead primes 1 mod 4 for every x
   below 26861, then fall behind (Leech, 1957).
3. pi(x) < li(x) for every real x in [2, 10^6] (and far beyond), although
   Littlewood proved that the inequality reverses infinitely often.
4. The Borwein integrals: an exact rational computation shows when the
   pattern "the integral equals pi/2" breaks, and by how little.
Searches 1-3 are exhaustive up to the stated bounds; li(x) is computed in
floating point with a convergent series.
"""

import math
from fractions import Fraction

import numpy as np

EULER_GAMMA = 0.5772156649015329

def euler_fifth_power_search(limit):
    """All (a, b, c, d, e) with a<=b<=c<=d<e<=limit and a^5+b^5+c^5+d^5 = e^5."""
    fifth = np.arange(limit + 1, dtype=np.int64) ** 5
    a, b = np.triu_indices(limit + 1, k=0)
    keep = (a >= 1)
    a, b = a[keep], b[keep]
    sums = fifth[a] + fifth[b]
    order = np.argsort(sums, kind="stable")
    sums, a, b = sums[order], a[order], b[order]
    found = set()
    for e in range(2, limit + 1):
        target = fifth[e] - sums                      # e^5 - (c^5 + d^5)
        ok = (target > 0)
        idx = np.searchsorted(sums, target[ok])
        idx = np.minimum(idx, len(sums) - 1)
        hits = np.nonzero(sums[idx] == target[ok])[0]
        cd = np.nonzero(ok)[0][hits]
        for i, j in zip(idx[hits], cd):
            quad = tuple(sorted((int(a[i]), int(b[i]), int(a[j]), int(b[j]))))
            if quad[3] < e:
                found.add(quad + (e,))
    return sorted(found)

def prime_sieve(limit):
    """Boolean array is_prime[0..limit]."""
    is_prime = np.ones(limit + 1, dtype=bool)
    is_prime[:2] = False
    for p in range(2, int(limit ** 0.5) + 1):
        if is_prime[p]:
            is_prime[p * p::p] = False
    return is_prime

def race_mod_4(limit):
    """First x where #primes = 1 mod 4 exceeds #primes = 3 mod 4, and the
    number of integers x <= limit at which it does."""
    is_prime = prime_sieve(limit)
    n = np.arange(limit + 1)
    ones = np.cumsum(is_prime & (n % 4 == 1))
    threes = np.cumsum(is_prime & (n % 4 == 3))
    ahead = np.nonzero(ones > threes)[0]
    return int(ahead[0]), int(len(ahead))

def li(x):
    """Logarithmic integral li(x) = gamma + ln ln x + sum (ln x)^k / (k k!)."""
    t = math.log(x)
    total, term, k = 0.0, 1.0, 1
    while True:
        term *= t / k
        add = term / k
        total += add
        if add < 1e-17 * total:
            break
        k += 1
    return EULER_GAMMA + math.log(t) + total

def li_array(xs):
    """li for a NumPy array of x > 1, by the same series (vectorised)."""
    t = np.log(np.asarray(xs, dtype=float))
    total, term = np.zeros_like(t), np.ones_like(t)
    for k in range(1, 200):
        term = term * t / k
        total += term / k
        if np.all(term / k < 1e-17 * total):
            break
    return EULER_GAMMA + np.log(t) + total

def li_exceeds_pi_up_to(limit):
    """Is li(x) > pi(x) for every real 2 <= x <= limit?  Since li increases and
    pi is constant between primes, it suffices to check x = each prime."""
    is_prime = prime_sieve(limit)
    primes = np.nonzero(is_prime)[0]
    margin = li_array(primes) - np.arange(1, len(primes) + 1)
    return bool(np.all(margin > 0)), float(margin.min())

def borwein_pattern(max_n):
    """For n = 1..max_n: (s_n, exact value of the n-th Borwein integral / pi).

    The n-th integral is that of prod_{k=0}^{n} sinc(x/(2k+1)) over (0, inf).
    It equals pi/2 while s_n = 1/3 + ... + 1/(2n+1) < 1; at the first n with
    s_n > 1 it equals (pi/2)(1 - (s_n - 1)^n / (2^(n-1) n! prod 1/(2k+1)))
    (D. and J. Borwein, 2001).
    """
    rows = []
    for n in range(1, max_n + 1):
        a = [Fraction(1, 2 * k + 1) for k in range(1, n + 1)]
        s = sum(a)
        if s < 1:
            rows.append((s, Fraction(1, 2)))
        else:
            prod = math.prod(a)
            deficit = (s - 1) ** n / (2 ** (n - 1) * math.factorial(n) * prod)
            rows.append((s, Fraction(1, 2) * (1 - deficit)))
            break
    return rows

def main():
    sols = euler_fifth_power_search(150)
    print("1. a^5 + b^5 + c^5 + d^5 = e^5, all 1 <= a<=b<=c<=d < e <= 150:")
    for a, b, c, d, e in sols:
        print(f"   {a}^5 + {b}^5 + {c}^5 + {d}^5 = {e}^5")
    print(f"   {len(sols)} solution(s): Euler's conjecture (1769) is false.")
    first, count = race_mod_4(10 ** 5)
    print("2. Prime race mod 4, x <= 100000: primes 1 mod 4 first lead at")
    print(f"   x = {first}; they lead at only {count} of the 100000 values of x.")
    is_prime = prime_sieve(10 ** 6)
    counts = np.cumsum(is_prime)
    print("3. pi(x) versus li(x):")
    for k in range(2, 7):
        x = 10 ** k
        print(f"   x = 10^{k}: pi(x) = {int(counts[x]):>6}, "
              f"li(x) = {li(x):>10.1f}, li(x) - pi(x) = {li(x) - counts[x]:6.1f}")
    holds, margin = li_exceeds_pi_up_to(10 ** 6)
    print(f"   li(x) > pi(x) for all real 2 <= x <= 10^6: {holds} "
          f"(smallest gap {margin:.3f})")
    print("4. Borwein integrals, exact values divided by pi:")
    for n, (s, value) in enumerate(borwein_pattern(10), start=1):
        tag = "= 1/2" if value == Fraction(1, 2) else \
            f"= 1/2 - {float(Fraction(1, 2) - value):.3e}"
        print(f"   n = {n}: 1/3 + ... + 1/{2 * n + 1} = {float(s):.4f}, "
              f"integral/pi {tag}")
    print()
    print("Epistemic status: searches 1 and 2 are exhaustive and exact up to their")
    print("stated bounds, and 4 is exact rational arithmetic on a proved formula.")
    print("Check 3 uses floating-point li(x) at every prime up to 10^6: numerical")
    print("evidence only - and Littlewood's theorem shows that this evidence,")
    print("however far it is extended, points the wrong way.")

if __name__ == "__main__":
    main()