erdos_solver_db.py

"""
erdos_solver_db.py — Auto-generating Erdős solver

Loads atlas operators. For each problem in the registry, applies the
appropriate pipeline. Stores results in SQLite. Generates new entries
on demand.

Author: Elias Oulad Brahim (Brahim Framework)
ORCID: 0009-0009-3302-9532
"""
import math
import sqlite3
import os
import json
import time
import sys
from math import gcd
from collections import defaultdict
from functools import lru_cache

try:
    from sympy import isprime, primerange, divisors, factorint, totient, mobius
    from math import comb as sym_comb
except ImportError:
    print("ERROR: sympy required. pip install sympy", file=sys.stderr)
    sys.exit(1)


# ============================================================
# FRAMEWORK CONSTANTS
# ============================================================
PHI = (1 + math.sqrt(5)) / 2
PSI = (1 - math.sqrt(5)) / 2
SIGMA = math.log(PHI)
K = 107
S_ATLAS = 214
N_C = 3
N_ST = 4
ICAS = 29
M_840 = 840
GF_841 = 841

BRAHIM = [27, 42, 60, 75, 97, 117, 139, 154, 172, 187]

MORDELL = {1, 121, 169, 289, 361, 529}        # = QR((Z/840)*)
PI_HEEG = {311, 479, 551, 671, 719, 839}      # = -QR((Z/840)*)
UNIT_PRIMES = [1, 11, 13, 17, 19, 23]


# ============================================================
# UTILITIES
# ============================================================
def L(n):
    """Lucas number L_n. L_0=2, L_1=1."""
    if n == 0: return 2
    if n == 1: return 1
    a, b = 2, 1
    for _ in range(n - 1):
        a, b = b, a + b
    return b


def F(n):
    """Fibonacci F_n. F_0=0, F_1=1."""
    if n == 0: return 0
    if n == 1: return 1
    a, b = 0, 1
    for _ in range(n - 1):
        a, b = b, a + b
    return b


def D_pitch(x):
    """Phi-graded depth: D(x) = -log(x)/log(phi)."""
    if x <= 0: return None
    return -math.log(x) / SIGMA


# ============================================================
# ATLAS OPERATORS
# ============================================================
class Atlas:
    """Atlas operator suite. Each method is a named operator."""
    
    @staticmethod
    def R_hat(n, mod=M_840):
        """Residue operator."""
        return n % mod
    
    @staticmethod
    def W_hat(n, mod=M_840):
        """W-mirror: x -> mod - x."""
        return (mod - n) % mod
    
    @staticmethod
    def M_K(p):
        """K-radial offset: p - K."""
        return p - K
    
    @staticmethod
    def Phi_Gal(x, p=ICAS, mod=M_840):
        """Frobenius x -> x^p mod N."""
        return pow(x, p, mod)
    
    @staticmethod
    def CRT_840(p):
        """CRT signature mod (8, 3, 5, 7)."""
        return (p % 8, p % 3, p % 5, p % 7)
    
    @staticmethod
    def is_QR_mod(a, n):
        """Is a a quadratic residue mod n?"""
        a_mod = a % n
        for x in range(n):
            if (x * x) % n == a_mod:
                return True
        return False
    
    @staticmethod
    def Mordell_root(p):
        """If p mod 840 is in Mordell, return q ∈ {1,11,13,17,19,23} with q² ≡ p mod 840."""
        r = p % M_840
        if r not in MORDELL: return None
        for q in UNIT_PRIMES:
            if (q * q) % M_840 == r: return q
        return None
    
    @staticmethod
    def Y_Z_decomp(p, x):
        """Given p and x, compute all (y, z) with 4/p = 1/x + 1/y + 1/z, x ≤ y ≤ z.
        Closed-form via divisor enumeration of b² where b = denom in lowest terms.
        Returns list of (y, z) pairs."""
        if 4 * x <= p: return []
        rem_num = 4 * x - p
        rem_den = p * x
        g = gcd(rem_num, rem_den)
        a = rem_num // g
        b = rem_den // g
        b2 = b * b
        pairs = []
        for u in divisors(b2):
            if u * u > b2: break
            v = b2 // u
            if (u + b) % a != 0: continue
            if (v + b) % a != 0: continue
            y = (u + b) // a
            z = (v + b) // a
            if y < x or y > z: continue
            if (y * z + x * z + x * y) * p == 4 * x * y * z:
                pairs.append((y, z))
        return pairs


# ============================================================
# ES SOLVER (the working pipeline)
# ============================================================
class ES_Solver:
    """Erdős-Straus 4/n = 1/x + 1/y + 1/z solver via 2-recipe pipeline."""
    
    KNOWN_EXCEPTIONS = {2521, 7681, 33601}
    
    @staticmethod
    def solve(p):
        """Return (x, y, z, recipe) or None."""
        if not isprime(p): return None
        
        x_min = (p + 3) // 4
        x_max = 3 * p // 4
        
        # Recipe 1: x ≡ 0 mod 8
        x_start = x_min + (8 - x_min % 8) % 8
        for x in range(x_start, x_max + 1, 8):
            pairs = Atlas.Y_Z_decomp(p, x)
            if pairs:
                y, z = pairs[0]
                return (x, y, z, "R1_mod8")
        
        # Recipe 2: boundary scan
        upper = min(x_min + 300, x_max)
        for x in range(x_min, upper + 1):
            pairs = Atlas.Y_Z_decomp(p, x)
            if pairs:
                y, z = pairs[0]
                return (x, y, z, "R2_boundary")
        
        # Recipe 3: exhaustive
        for x in range(x_min, x_max + 1):
            pairs = Atlas.Y_Z_decomp(p, x)
            if pairs:
                y, z = pairs[0]
                return (x, y, z, "R3_exhaustive")
        
        return None
    
    @staticmethod
    def classify(p):
        """Return structural classification dict."""
        if not isprime(p):
            return {"prime": False}
        r = p % M_840
        return {
            "prime": True,
            "p": p,
            "residue_mod_840": r,
            "is_hard": r in MORDELL,
            "Mordell_root": Atlas.Mordell_root(p),
            "W_mirror_residue": Atlas.W_hat(r),
            "W_mirror_in_Pi_Heeg": Atlas.W_hat(r) in PI_HEEG,
            "CRT_signature": Atlas.CRT_840(p),
            "M_K_offset": Atlas.M_K(p),
            "D_pitch": D_pitch(p),
            "S_p_atlas_bound": p ** (1 / N_C),
            "search_range": (x_min := (p + 3) // 4, 3 * p // 4),
        }


# ============================================================
# OTHER PROBLEM SOLVERS
# ============================================================
class GoldbachSolver:
    """G(2c) ≥ 1: find Goldbach pair."""
    
    @staticmethod
    def solve(c):
        n = 2 * c
        if n < 4: return None
        for p in primerange(2, c + 1):
            q = n - p
            if isprime(q):
                return (p, q)
        return None
    
    @staticmethod
    def count(c):
        n = 2 * c
        return sum(1 for p in primerange(2, c + 1) if isprime(n - p))


class TwinPrimeCounter:
    """π_2(N) = #{p ≤ N : p, p+2 both prime}."""
    
    @staticmethod
    def count(N):
        return sum(1 for p in primerange(2, N) if isprime(p + 2))
    
    @staticmethod
    def HL_estimate(N):
        """Hardy-Littlewood: π_2(N) ~ 2 C_2 N/(log N)²."""
        if N < 3: return 0
        # C_2 hardcoded for speed; computed as ∏_{p>2}(1 - 1/(p-1)²) over primes > 2
        C2 = 0.6601682965
        return 2 * C2 * N / (math.log(N)) ** 2


class BrocardSolver:
    """n!+1 = m². Verify or check."""
    
    @staticmethod
    def check(n):
        f = math.factorial(n)
        val = f + 1
        m = int(math.isqrt(val))
        return (m, m * m == val)
    
    @staticmethod
    def find_solutions(max_n=25):
        return [(n, m) for n in range(1, max_n + 1)
                for m, ok in [BrocardSolver.check(n)] if ok]


class BealSolver:
    """x^a + y^b = z^c, a,b,c≥3. Conjecture: gcd(x,y,z) > 1 always."""
    
    @staticmethod
    def search_counterexamples(max_base=20, max_exp=5):
        ces = []
        for x in range(1, max_base):
            for y in range(1, max_base):
                for a in range(3, max_exp + 1):
                    for b in range(3, max_exp + 1):
                        val = x ** a + y ** b
                        for c in range(3, max_exp + 1):
                            z = round(val ** (1 / c))
                            if z > 0 and z ** c == val:
                                g = gcd(gcd(x, y), z)
                                if g == 1:
                                    ces.append((x, a, y, b, z, c))
        return ces


class CatalanSolver:
    """x^p - y^q = 1. Mihailescu: only (3,2,2,3)."""
    
    @staticmethod
    def search(max_base=50, max_exp=7):
        sols = []
        for x in range(2, max_base):
            for p in range(2, max_exp + 1):
                for y in range(2, max_base):
                    for q in range(2, max_exp + 1):
                        if x ** p - y ** q == 1:
                            sols.append((x, p, y, q))
        return sols


class RamseyR3kSolver:
    """R(3, k) — known values + framework anchors."""
    
    KNOWN = {3: 6, 4: 9, 5: 14, 6: 18, 7: 23, 8: 28, 9: 36}
    
    @staticmethod
    def value(k):
        return RamseyR3kSolver.KNOWN.get(k)
    
    @staticmethod
    def lower_bound(k):
        if k < 2: return None
        return k * k / (math.log(k)) ** 2
    
    @staticmethod
    def upper_bound(k):
        return sym_comb(k - 1 + 2, 2)


class SidonSolver:
    """Sidon set: pairwise sums distinct. Perfect difference set search."""
    
    @staticmethod
    def is_sidon(s):
        sums = set()
        for i in range(len(s)):
            for j in range(i, len(s)):
                v = s[i] + s[j]
                if v in sums: return False
                sums.add(v)
        return True


# ============================================================
# DATABASE
# ============================================================
class ErdosSolverDB:
    """SQLite-backed database of solver results. Auto-populates."""
    
    SCHEMA = """
    CREATE TABLE IF NOT EXISTS problems (
        problem_id TEXT PRIMARY KEY,
        name TEXT NOT NULL,
        family TEXT NOT NULL,
        equation TEXT,
        atlas_operators TEXT,
        status TEXT,
        notes TEXT
    );
    
    CREATE TABLE IF NOT EXISTS es_solutions (
        prime INTEGER PRIMARY KEY,
        residue INTEGER,
        is_hard INTEGER,
        x INTEGER, y TEXT, z TEXT,
        recipe TEXT,
        elapsed_ms REAL,
        verified INTEGER
    );
    
    CREATE TABLE IF NOT EXISTS goldbach_pairs (
        c INTEGER PRIMARY KEY,
        n_2c INTEGER,
        p1 INTEGER, p2 INTEGER,
        n_pairs INTEGER
    );
    
    CREATE TABLE IF NOT EXISTS twin_primes (
        N INTEGER PRIMARY KEY,
        count INTEGER,
        HL_estimate REAL,
        ratio REAL
    );
    
    CREATE TABLE IF NOT EXISTS brocard (
        n INTEGER PRIMARY KEY,
        n_factorial_plus_1 TEXT,
        is_square INTEGER,
        m INTEGER
    );
    
    CREATE TABLE IF NOT EXISTS atlas_anchors (
        anchor_id TEXT PRIMARY KEY,
        problem_family TEXT,
        anchor_input TEXT,
        anchor_output TEXT,
        framework_form TEXT,
        verified INTEGER
    );
    """
    
    def __init__(self, path):
        self.path = path
        self.conn = sqlite3.connect(path)
        self.conn.executescript(self.SCHEMA)
        self.conn.commit()
        self._populate_problems()
        self._populate_anchors()
    
    def _populate_problems(self):
        problems = [
            ("ES_4n", "Erdős-Straus 4/n", "diophantine_unit_fractions",
             "4/n = 1/x + 1/y + 1/z", "Atlas.R_hat, Atlas.W_hat, Atlas.M_K, Atlas.Y_Z_decomp",
             "OPEN; 100% pipeline coverage to 5e7", "QR((Z/840)*) hard residues"),
            ("Goldbach_binary", "Goldbach binary", "additive_prime",
             "G(2c) ≥ 1, p+q=2c both prime", "GP_hat, M_K antisymmetric",
             "OPEN; verified to 4e18", "G(214)=8=2^N_c anchor"),
            ("twin_primes", "Twin primes", "additive_prime",
             "infinitely many p with p+2 prime", "GP_hat gap-2",
             "OPEN; bounded gaps ≤ 246 (Maynard)", "no closed-form anchor"),
            ("Brocard", "Brocard n!+1=m²", "diophantine",
             "n! + 1 = m²", "P_phi, KB_hat",
             "OPEN beyond n=7", "7! = 70·72 anchor"),
            ("Beal", "Beal conjecture", "diophantine_power",
             "x^a + y^b = z^c, a,b,c≥3, gcd=1: no solution",
             "P_phi, B_hat", "OPEN", "no counterexample known"),
            ("Catalan", "Catalan-Mihailescu", "diophantine_power",
             "x^p - y^q = 1: only (3,2,2,3)", "N_c anchor",
             "SOLVED 2002", "N_c² - 2^N_c = 1 atlas reading"),
            ("Ramsey_R3k", "Ramsey R(3,k)", "combinatorial_extremal",
             "R(3,k) growth", "Sigma, Casimir",
             "OPEN; CJMS 2025 lower bound", "R(3,3)=2N_c, R(3,4)=N_c²"),
            ("EKR", "Erdős-Ko-Rado", "combinatorial_extremal",
             "|F| ≤ C(n-1, k-1) intersecting", "Sigma, P_+",
             "SOLVED (Wilson 1984)", "n=29=|I|, k=5 framework anchor"),
            ("Turan", "Turán's theorem", "combinatorial_extremal",
             "ex(n; K_{r+1}) = (1-1/r) n²/2", "Sigma",
             "SOLVED 1941", "density 1-1/N_c=2/3"),
            ("Sidon", "Sidon sets", "combinatorial",
             "pairwise distinct sums", "Sigma, GP_hat",
             "various OPEN extensions", "framework integers not Sidon"),
        ]
        cur = self.conn.cursor()
        for p in problems:
            cur.execute(
                "INSERT OR IGNORE INTO problems VALUES (?,?,?,?,?,?,?)", p)
        self.conn.commit()
    
    def _populate_anchors(self):
        anchors = [
            ("ES_QR_840", "ES_4n", "H_840",
             "{1,121,169,289,361,529}",
             "QR((Z/840)*) — exactly 6 = 1·1·2·3 from CRT", 1),
            ("ES_W_PiHeeg", "ES_4n", "W(H_840)",
             "{311,479,551,671,719,839}",
             "Pi_Heeg = -QR((Z/840)*)", 1),
            ("Goldbach_S", "Goldbach_binary", "G(214)",
             "8", "G(S) = 2^N_c", 1),
            ("Ramsey_Nc", "Ramsey_R3k", "R(3,3)",
             "6", "R(3, N_c) = 2·N_c", 1),
            ("Ramsey_Nst", "Ramsey_R3k", "R(3,4)",
             "9", "R(3, N_st) = N_c²", 1),
            ("Catalan_Nc", "Catalan", "(3,2,2,3)",
             "3² - 2³ = 1", "N_c² - 2^N_c = 1", 1),
            ("Turan_Nc", "Turan", "density at r=N_c",
             "2/3", "1 - 1/N_c", 1),
            ("Brocard_7fact", "Brocard", "7!", "5040 = 70·72",
             "(71-1)(71+1) = 7! at largest known Brocard m", 1),
            ("EKR_29_5", "EKR", "(n,k)=(29,5)",
             f"{int(sym_comb(28, 4))}", "n=|I|, k=Sigma 5-fold", 1),
        ]
        cur = self.conn.cursor()
        for a in anchors:
            cur.execute(
                "INSERT OR IGNORE INTO atlas_anchors VALUES (?,?,?,?,?,?)", a)
        self.conn.commit()
    
    # ============================================================
    # AUTO-GENERATING METHODS
    # ============================================================
    def populate_es(self, p_min=2, p_max=10000, hard_only=True):
        """Compute ES solutions for primes in [p_min, p_max]."""
        cur = self.conn.cursor()
        count = 0
        for p in primerange(p_min, p_max + 1):
            r = p % M_840
            is_hard = r in MORDELL
            if hard_only and not is_hard: continue
            
            # Skip if already in DB
            existing = cur.execute(
                "SELECT prime FROM es_solutions WHERE prime=?", (p,)).fetchone()
            if existing: continue
            
            t0 = time.time()
            sol = ES_Solver.solve(p)
            elapsed_ms = (time.time() - t0) * 1000
            
            if sol is None:
                cur.execute(
                    "INSERT INTO es_solutions VALUES (?,?,?,?,?,?,?,?,?)",
                    (p, r, int(is_hard), None, None, None, "FAIL", elapsed_ms, 0))
            else:
                x, y, z, recipe = sol
                # Verify
                verified = (y*z + x*z + x*y) * p == 4 * x * y * z
                cur.execute(
                    "INSERT INTO es_solutions VALUES (?,?,?,?,?,?,?,?,?)",
                    (p, r, int(is_hard), x, str(y), str(z), recipe, elapsed_ms, int(verified)))
            count += 1
        self.conn.commit()
        return count
    
    def populate_goldbach(self, c_min=2, c_max=1000):
        cur = self.conn.cursor()
        count = 0
        for c in range(c_min, c_max + 1):
            existing = cur.execute(
                "SELECT c FROM goldbach_pairs WHERE c=?", (c,)).fetchone()
            if existing: continue
            pair = GoldbachSolver.solve(c)
            n_pairs = GoldbachSolver.count(c)
            if pair:
                p1, p2 = pair
                cur.execute(
                    "INSERT INTO goldbach_pairs VALUES (?,?,?,?,?)",
                    (c, 2*c, p1, p2, n_pairs))
            else:
                cur.execute(
                    "INSERT INTO goldbach_pairs VALUES (?,?,?,?,?)",
                    (c, 2*c, None, None, n_pairs))
            count += 1
        self.conn.commit()
        return count
    
    def populate_twin_primes(self, scales=None):
        if scales is None:
            scales = [100, 1000, 10000, 100000, 1000000]
        cur = self.conn.cursor()
        count = 0
        for N in scales:
            existing = cur.execute(
                "SELECT N FROM twin_primes WHERE N=?", (N,)).fetchone()
            if existing: continue
            actual = TwinPrimeCounter.count(N)
            est = TwinPrimeCounter.HL_estimate(N)
            ratio = actual / est if est > 0 else 0
            cur.execute(
                "INSERT INTO twin_primes VALUES (?,?,?,?)",
                (N, actual, est, ratio))
            count += 1
        self.conn.commit()
        return count
    
    def populate_brocard(self, n_max=25):
        cur = self.conn.cursor()
        count = 0
        for n in range(1, n_max + 1):
            existing = cur.execute(
                "SELECT n FROM brocard WHERE n=?", (n,)).fetchone()
            if existing: continue
            f1 = math.factorial(n) + 1
            m, is_sq = BrocardSolver.check(n)
            cur.execute(
                "INSERT INTO brocard VALUES (?,?,?,?)",
                (n, str(f1), int(is_sq), m if is_sq else None))
            count += 1
        self.conn.commit()
        return count
    
    def autogenerate_all(self, p_max=10000, c_max=1000, brocard_max=25):
        """Run all populators."""
        report = {}
        print(f"Auto-generating ES solutions for primes ≤ {p_max} (hard only)...")
        report['es'] = self.populate_es(p_max=p_max)
        print(f"  Added {report['es']} ES entries")
        
        print(f"Auto-generating Goldbach pairs for c ≤ {c_max}...")
        report['goldbach'] = self.populate_goldbach(c_max=c_max)
        print(f"  Added {report['goldbach']} Goldbach entries")
        
        print(f"Auto-generating twin prime counts...")
        report['twin'] = self.populate_twin_primes()
        print(f"  Added {report['twin']} twin prime entries")
        
        print(f"Auto-generating Brocard checks for n ≤ {brocard_max}...")
        report['brocard'] = self.populate_brocard(n_max=brocard_max)
        print(f"  Added {report['brocard']} Brocard entries")
        
        return report
    
    def stats(self):
        cur = self.conn.cursor()
        result = {}
        result['problems'] = cur.execute(
            "SELECT COUNT(*) FROM problems").fetchone()[0]
        result['atlas_anchors'] = cur.execute(
            "SELECT COUNT(*) FROM atlas_anchors").fetchone()[0]
        result['es_solutions'] = cur.execute(
            "SELECT COUNT(*) FROM es_solutions").fetchone()[0]
        result['es_hard_solved'] = cur.execute(
            "SELECT COUNT(*) FROM es_solutions WHERE is_hard=1 AND x IS NOT NULL").fetchone()[0]
        result['es_failures'] = cur.execute(
            "SELECT COUNT(*) FROM es_solutions WHERE recipe='FAIL'").fetchone()[0]
        result['goldbach_pairs'] = cur.execute(
            "SELECT COUNT(*) FROM goldbach_pairs").fetchone()[0]
        result['twin_prime_scales'] = cur.execute(
            "SELECT COUNT(*) FROM twin_primes").fetchone()[0]
        result['brocard_solutions'] = cur.execute(
            "SELECT COUNT(*) FROM brocard WHERE is_square=1").fetchone()[0]
        return result
    
    def close(self):
        self.conn.close()


# ============================================================
# CLI
# ============================================================
if __name__ == "__main__":
    import argparse
    parser = argparse.ArgumentParser(description="Erdős solver auto-DB")
    parser.add_argument("--db", default="erdos_solver.sqlite",
                        help="SQLite path")
    parser.add_argument("--p_max", type=int, default=10000,
                        help="Max prime for ES population")
    parser.add_argument("--c_max", type=int, default=1000,
                        help="Max c for Goldbach")
    parser.add_argument("--brocard_max", type=int, default=25,
                        help="Max n for Brocard")
    parser.add_argument("--autogen", action="store_true",
                        help="Run auto-generation")
    parser.add_argument("--stats", action="store_true",
                        help="Print DB stats")
    parser.add_argument("--solve", type=int, metavar="P",
                        help="Solve ES for specific prime P")
    parser.add_argument("--classify", type=int, metavar="P",
                        help="Classify prime P via atlas")
    
    args = parser.parse_args()
    
    db = ErdosSolverDB(args.db)
    
    if args.solve is not None:
        sol = ES_Solver.solve(args.solve)
        if sol:
            x, y, z, rec = sol
            print(f"4/{args.solve} = 1/{x} + 1/{y} + 1/{z}  (recipe {rec})")
        else:
            print(f"No solution found for p={args.solve}")
    
    if args.classify is not None:
        cls = ES_Solver.classify(args.classify)
        print(json.dumps(cls, indent=2, default=str))
    
    if args.autogen:
        report = db.autogenerate_all(
            p_max=args.p_max, c_max=args.c_max, brocard_max=args.brocard_max)
        print(f"\nAuto-generation complete: {report}")
    
    if args.stats:
        s = db.stats()
        print(f"\nDatabase stats: {args.db}")
        for k, v in s.items():
            print(f"  {k}: {v}")
    
    db.close()