#!/usr/bin/env python3
"""
base4096_hunter.py — PATH 3: Arbitrary-Precision Mersenne Candidate Scorer

DESIGN PATH SEPARATION:
  Path 1 (resonance)  = mersenne_sieve.c + pln_final.c  — generalize
  Path 2 (overfit)    = bigg_overfit.c                   — memorize
  Path 3 (this file)  = precision-unlimited scoring      — scale

  This tool uses Python's arbitrary-precision integers and the base4096
  encoding layer for compact I/O at any precision. It re-implements the
  PLN coordinate system and D_n(r) operator with full Python mpmath
  precision (unlimited digits), enabling:

    1. Extended-precision PLN coordinates for huge exponents where
       C's long double loses significant bits.
    2. Base4096 compact storage/retrieval of Mersenne candidate lists
       and fingerprints.
    3. Future expansion: GPU-friendly base4096-encoded arithmetic
       via the OpenGL-compatible alphabet.

Usage:
  python base4096_hunter.py --pln <p1> [p2 ...]       PLN coords at full precision
  python base4096_hunter.py --encode-candidates <file> Encode candidate list to base4096
  python base4096_hunter.py --decode-candidates <file> Decode base4096 candidate list
  python base4096_hunter.py --dn <n> <r> <omega>      Compute D_n(r) at full precision
  python base4096_hunter.py --fingerprint <p1> [p2 ..] Full HDGL fingerprint
  python base4096_hunter.py --help
"""

import sys
import os
import struct
import math

# Add base4096 library to path
SCRIPT_DIR = os.path.dirname(os.path.abspath(__file__))
BASE4096_DIR = os.path.join(SCRIPT_DIR, "base4096-2.0.1")
sys.path.insert(0, BASE4096_DIR)

try:
    import base4096
except ImportError:
    print("ERROR: base4096 library not found in base4096-2.0.1/")
    print("       Make sure base4096.py and frozen_base4096_alphabet.txt are present.")
    sys.exit(1)

# Try mpmath for arbitrary precision; fall back to math
try:
    from mpmath import mp, mpf, log, sqrt, cos, sin, fmod, floor, fabs, exp, pi, power
    mp.dps = 50  # 50 decimal digits of precision by default
    HAVE_MPMATH = True
except ImportError:
    HAVE_MPMATH = False
    # Shim: use stdlib math (limited to float64)
    from math import log, sqrt, cos, sin, fmod, floor, fabs, exp, pi
    mpf = float
    def power(x, y): return x ** y

# ================================================================
# CONSTANTS (identical to hdgl_analog_v30.c / bigg_overfit.c)
# ================================================================
PHI = mpf("1.6180339887498948482045868343656381177203091798057628")
LN_PHI = mpf("0.48121182505960344749775891342436842076978601214440582")
LN2 = mpf("0.69314718055994530941723212145817656807550013436025525")
LOG10_2 = mpf("0.30102999566398119521373889472449302676818596588662250")

FIB_TABLE = [1, 1, 2, 3, 5, 8, 13, 21]
PRIME_TABLE = [2, 3, 5, 7, 11, 13, 17, 19]

UPPER_FIELD = [
    mpf("170.6180339887"), mpf("150.9442719100"), mpf("12.6180339887"),
    mpf("8.8541019662"), mpf("4.2360679775"), mpf("3.6180339887"), mpf("1.6180339887")
]
LOWER_FIELD = [
    mpf("0.0000000001"), mpf("0.0344465435"), mpf("0.0557280900"), mpf("0.0901699437"),
    mpf("0.1458980338"), mpf("0.2360679775"), mpf("0.3819660113"), mpf("0.6180339887")
]
SIBLING_HARMONICS = [
    mpf("0.0901699437"), mpf("0.1458980338"), mpf("0.2360679775"), mpf("0.3090169944"),
    mpf("0.3819660113"), mpf("0.4721359549"), mpf("0.6545084972"), mpf("0.8729833462")
]
BASE_INF_SEEDS = [
    0.6180339887, 1.6180339887, 2.6180339887, 3.6180339887, 4.8541019662,
    5.6180339887, 6.4721359549, 7.8541019662, 8.3141592654, 0.0901699437,
    0.1458980338, 0.2360679775, 0.3090169944, 0.3819660113, 0.4721359549,
    0.6545084972, 0.8729833462, 1.0000000000, 1.2360679775, 1.6180339887,
    2.2360679775, 2.6180339887, 3.1415926535, 3.6180339887, 4.2360679775,
    4.8541019662, 5.6180339887, 6.4721359549, 7.2360679775, 7.8541019662,
    8.6180339887, 9.2360679775, 9.8541019662, 10.6180339887, 11.0901699437,
    11.9442719100, 12.6180339887, 13.6180339887, 14.2360679775, 14.8541019662,
    15.6180339887, 16.4721359549, 17.2360679775, 17.9442719100, 18.6180339887,
    19.2360679775, 19.8541019662, 20.6180339887, 21.0901699437, 21.9442719100,
    22.6180339887, 23.6180339887, 24.2360679775, 24.8541019662, 25.6180339887,
    26.4721359549, 27.2360679775, 27.9442719100, 28.6180339887, 29.0344465435,
    29.6180339887, 30.2360679775, 30.8541019662, 31.6180339887
]

ALL_MERSENNE_P = [
    2, 3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127,
    521, 607, 1279, 2203, 2281, 3217, 4253, 4423,
    9689, 9941, 11213, 19937, 21701, 23209, 44497,
    86243, 110503, 132049, 216091, 756839, 859433,
    1257787, 1398269, 2976221, 3021377, 6972593,
    13466917, 20996011, 24036583, 25964951, 30402457,
    32582657, 37156667, 42643801, 43112609,
    57885161, 74207281, 77232917, 82589933, 136279841
]


# ================================================================
# D_n(r) — Exact formula from hdgl_analog_v30.c at arbitrary precision
# ================================================================
def compute_Dn_r(n, r, omega):
    """D_n(r) = sqrt(phi * F_n * 2^n * P_n * omega) * r^k"""
    if n < 1 or n > 8:
        return mpf(0)
    idx = n - 1
    F_n = mpf(FIB_TABLE[idx])
    two_n = power(mpf(2), n)
    P_n = mpf(PRIME_TABLE[idx])
    k = mpf(n + 1) / mpf(8)
    base_val = sqrt(PHI * F_n * two_n * P_n * omega)
    r_power = power(fabs(r), k)
    return base_val * r_power


# ================================================================
# PLN Coordinates at arbitrary precision
# ================================================================
def orbit_r(omega):
    frac_part = omega - floor(omega)
    theta = frac_part * mpf(2) * pi
    D_cont = fmod(omega, mpf(8)) + mpf(1)
    lat = PHI * PHI * (mpf(1) - power(PHI, -D_cont)) + mpf(1) / D_cont
    lat = max(mpf(1), min(lat, PHI * PHI))
    Rre = cos(theta) - mpf(1) + lat
    Rim = sin(theta)
    return sqrt(Rre * Rre + Rim * Rim)


def compute_pln(p):
    """Full PLN coordinate at arbitrary precision."""
    p_mpf = mpf(p)
    lnx = p_mpf * LN2
    arg = lnx / LN_PHI
    omega = log(arg) / LN_PHI if arg > 1 else mpf(0)
    D = int(fmod(floor(omega), mpf(8))) + 1
    theta = fmod(omega, mpf(1))
    r = orbit_r(omega)
    digits = int(floor(p_mpf * LOG10_2)) + 1
    return {
        'p': p,
        'omega': omega,
        'D': D,
        'theta': theta,
        'r': r,
        'digits': digits
    }


# ================================================================
# Base4096 Candidate Encoding/Decoding
#
# Format: each exponent is encoded as a big-endian 8-byte uint64,
# then the entire blob is base4096-encoded for compact storage.
# ================================================================
def encode_candidates(exponents):
    """Encode a list of uint64 exponents to base4096 string.
    Prepend 0xFF sentinel + 4-byte count to prevent leading-zero loss."""
    data = b'\xff' + struct.pack('>I', len(exponents))
    data += b''.join(struct.pack('>Q', p) for p in exponents)
    return base4096.encode(data)


def decode_candidates(encoded_str):
    """Decode base4096 string back to list of uint64 exponents."""
    data = base4096.decode(encoded_str)
    if data[0] != 0xff:
        raise ValueError("Invalid sentinel byte — data may be corrupted")
    count = struct.unpack('>I', data[1:5])[0]
    exponents = []
    for i in range(count):
        p = struct.unpack('>Q', data[5 + i*8:5 + (i+1)*8])[0]
        exponents.append(p)
    return exponents


# ================================================================
# Fingerprint (full HDGL, matches bigg_overfit.c but at arb precision)
# ================================================================
def compute_fingerprint(p):
    pln = compute_pln(p)
    features = []

    # D_n(r) for all 8 dimensions
    for n in range(1, 9):
        features.append(float(compute_Dn_r(n, pln['r'], pln['omega'])))

    # omega, theta, orbit_r, D
    features.append(float(pln['omega']))
    features.append(float(pln['theta']))
    features.append(float(pln['r']))
    features.append(float(pln['D']))

    # Lattice projections
    features.append(float(fmod(pln['omega'] * UPPER_FIELD[0], mpf(1))))
    features.append(float(fmod(pln['omega'] * UPPER_FIELD[6], mpf(1))))
    features.append(float(fmod(pln['r'] * UPPER_FIELD[3], mpf(1))))

    # Lower field correlations
    features.append(float(fmod(pln['theta'] + LOWER_FIELD[7], mpf(1))))
    features.append(float(fmod(pln['omega'] * LOWER_FIELD[5], mpf(1))))
    features.append(float(fabs(sin(pln['omega'] * LOWER_FIELD[6] * pi))))

    # Sibling harmonic projections
    for i in range(4):
        features.append(float(fabs(cos(pln['omega'] * SIBLING_HARMONICS[i*2] * pi))))

    # Analog dimension phase locks (using first 4)
    ANALOG_DIMS = [mpf("8.3141592654"), mpf("7.8541019662"),
                   mpf("6.4721359549"), mpf("5.6180339887")]
    for i in range(4):
        phase = fmod(pln['omega'], ANALOG_DIMS[i])
        features.append(float(phase / ANALOG_DIMS[i]))

    # Base-inf seed distances
    omega_mod32 = float(fmod(pln['omega'], mpf(32)))
    dists = sorted(abs(omega_mod32 - s) for s in BASE_INF_SEEDS)
    features.extend(dists[:4])

    # D_n cross ratio
    d1 = float(compute_Dn_r(1, pln['r'], pln['omega']))
    d4 = float(compute_Dn_r(4, pln['r'], pln['omega']))
    d5 = float(compute_Dn_r(5, pln['r'], pln['omega']))
    d8 = float(compute_Dn_r(8, pln['r'], pln['omega']))
    features.append((d1 * d8) / (d4 * d5) if d4 * d5 > 1e-15 else 0.0)

    # Golden depth
    golden_dist = float(pln['r']) - 1.0 / float(PHI)
    features.append(math.exp(-golden_dist * golden_dist * 10.0))

    return features


# ================================================================
# CLI
# ================================================================
def cmd_pln(exponents):
    prec_label = f"mpmath {mp.dps} digits" if HAVE_MPMATH else "float64 (install mpmath for more)"
    print(f"=== PLN COORDINATES — {prec_label} ===\n")
    print(f"  {'p':>12}  D   {'omega':>20}  {'theta':>20}  {'|R|':>20}  {'digits':>12}")
    print(f"  {'─'*12}  ──  {'─'*20}  {'─'*20}  {'─'*20}  {'─'*12}")

    for p in exponents:
        pln = compute_pln(p)
        known = "KNOWN" if p in ALL_MERSENNE_P else ""
        if HAVE_MPMATH:
            print(f"  {p:>12}  D{pln['D']}  {mp.nstr(pln['omega'], 18):>20}  "
                  f"{mp.nstr(pln['theta'], 18):>20}  {mp.nstr(pln['r'], 18):>20}  "
                  f"{pln['digits']:>12}  {known}")
        else:
            print(f"  {p:>12}  D{pln['D']}  {float(pln['omega']):>20.15f}  "
                  f"{float(pln['theta']):>20.15f}  {float(pln['r']):>20.15f}  "
                  f"{pln['digits']:>12}  {known}")


def cmd_dn(n, r, omega):
    result = compute_Dn_r(n, mpf(r), mpf(omega))
    prec = f"mpmath {mp.dps} digits" if HAVE_MPMATH else "float64"
    print(f"D_{n}(r={r}, omega={omega}) = {result}")
    print(f"Precision: {prec}")
    print(f"Formula: sqrt(phi * F_{n} * 2^{n} * P_{n} * omega) * |r|^((n+1)/8)")


def cmd_encode(filepath):
    """Read exponents from a text file (one per line), encode to base4096."""
    exponents = []
    with open(filepath, 'r') as f:
        for line in f:
            line = line.strip()
            if line and not line.startswith('#'):
                try:
                    exponents.append(int(line))
                except ValueError:
                    pass

    encoded = encode_candidates(exponents)
    outpath = filepath + ".b4096"
    with open(outpath, 'w', encoding='utf-8') as f:
        f.write(encoded)

    ratio = len(encoded) / max(1, os.path.getsize(filepath))
    print(f"Encoded {len(exponents)} candidates → {outpath}")
    print(f"  Raw file : {os.path.getsize(filepath):,} bytes")
    print(f"  Base4096 : {len(encoded)} chars ({len(encoded.encode('utf-8')):,} bytes UTF-8)")
    print(f"  Ratio    : {ratio:.2%}")


def cmd_decode(filepath):
    """Read base4096 file, decode back to exponent list."""
    with open(filepath, 'r', encoding='utf-8') as f:
        encoded = f.read().strip()

    exponents = decode_candidates(encoded)
    print(f"Decoded {len(exponents)} candidates from {filepath}:")
    for p in exponents:
        known = " ** KNOWN MERSENNE **" if p in ALL_MERSENNE_P else ""
        print(f"  {p}{known}")


def cmd_fingerprint(exponents):
    print(f"=== HDGL FINGERPRINTS (32 features) ===\n")
    feat_names = [
        "D1(r)","D2(r)","D3(r)","D4(r)","D5(r)","D6(r)","D7(r)","D8(r)",
        "omega","theta","orbit_r","D_dim",
        "upper_0","upper_6","upper_3",
        "lower_7","lower_5","lower_6",
        "sib_0","sib_2","sib_4","sib_6",
        "alock_0","alock_1","alock_2","alock_3",
        "inf_0","inf_1","inf_2","inf_3",
        "Dn_xratio","golden_dep"
    ]

    for p in exponents:
        fp = compute_fingerprint(p)
        known = "KNOWN" if p in ALL_MERSENNE_P else ""
        print(f"  p = {p}  {known}")
        for i, (name, val) in enumerate(zip(feat_names, fp)):
            print(f"    [{i:2d}] {name:14s} = {val:18.12f}")
        print()


def main():
    if len(sys.argv) < 2:
        print(__doc__)
        return

    cmd = sys.argv[1]

    if cmd in ('--help', '-h'):
        print(__doc__)
        return

    if cmd == '--pln' and len(sys.argv) >= 3:
        exponents = [int(x) for x in sys.argv[2:]]
        cmd_pln(exponents)
        return

    if cmd == '--dn' and len(sys.argv) >= 5:
        n = int(sys.argv[2])
        r = float(sys.argv[3])
        omega = float(sys.argv[4])
        cmd_dn(n, r, omega)
        return

    if cmd == '--encode-candidates' and len(sys.argv) >= 3:
        cmd_encode(sys.argv[2])
        return

    if cmd == '--decode-candidates' and len(sys.argv) >= 3:
        cmd_decode(sys.argv[2])
        return

    if cmd == '--fingerprint' and len(sys.argv) >= 3:
        exponents = [int(x) for x in sys.argv[2:]]
        cmd_fingerprint(exponents)
        return

    print(__doc__)


if __name__ == '__main__':
    main()
