import numpy as np
from sympy import primerange

# Constants
phi = (1 + np.sqrt(5)) / 2

# Fibonacci generator (simple)
def fib(n):
    a, b = 0, 1
    for _ in range(n):
        a, b = b, a + b
    return a if n > 0 else 1

# Get prime list (first 50 primes)
prime_list = list(primerange(2, 230))  # Enough primes

# Symbolic reconstruction function
def symbolic_param(p_phi, f_fib, q_prime):
    f_fib = max(0, int(round(f_fib)))
    q_prime = max(0, int(round(q_prime)))
    fib_val = fib(f_fib)
    prime_val = prime_list[q_prime] if q_prime < len(prime_list) else 1
    return (phi ** p_phi) * fib_val * prime_val

# Find best symbolic approx for a single parameter float
def fit_symbolic_param(value, max_fib=20, max_prime=30):
    best_error = float('inf')
    best_params = (0, 0, 0)

    for f in range(max_fib + 1):
        fib_val = fib(f)
        for q in range(max_prime + 1):
            prime_val = prime_list[q]
            # Avoid zero division or tiny numbers:
            base = fib_val * prime_val
            if base == 0:
                continue
            # Solve for p_phi: value = phi^p_phi * base => p_phi = log(value/base)/log(phi)
            if value / base <= 0:
                continue
            p_phi = np.log(value / base) / np.log(phi)
            approx = symbolic_param(p_phi, f, q)
            error = abs(approx - value)
            if error < best_error:
                best_error = error
                best_params = (p_phi, f, q)

    return best_params, best_error

# Your best-fit parameters (example from your printout)
best_fit_params = {
    'k': 1.049342,
    'r0': 1.049676,
    'Omega0': 1.049675,
    's0': 0.994533,
    'alpha': 0.340052,
    'beta': 0.360942,
    'gamma': 0.993975,
}

print("Symbolic decomposition of parameters:")
for name, val in best_fit_params.items():
    (p_phi, f_fib, q_prime), err = fit_symbolic_param(val)
    fib_val = fib(f_fib)
    prime_val = prime_list[q_prime]
    approx = symbolic_param(p_phi, f_fib, q_prime)
    print(f"{name:6}: {val:.6f} ≈ φ^{p_phi:.3f} * Fib({f_fib})={fib_val} * Prime({q_prime})={prime_val} = {approx:.6f} (err={err:.6e})")
