problem_solvers/ramanujan_machine.py

'''
Ramanujan Machine-style Continued Fraction Search
====================================================
Inspired by arXiv:2308.11829 (Elimelech et al.)
Key insight: polynomial continued fractions converging to mathematical constants
exhibit "factorial reduction" — gcd(p_n, q_n) grows super-exponentially.
This property acts as a FILTER that dramatically reduces the search space.

We adapt this to search for continued fractions converging to zeta-related values:
- ζ(2) = π²/6
- ζ(3) (Apéry's constant)
- Values at zero ordinates: ζ(1/2 + iγ_n)
- Li's criterion values λ_n

The "conservative matrix field" structure from the paper suggests that formulas
for related constants share matrix structures. We search for polynomial continued
fractions and test if they converge to known zeta values.
'''

import numpy as np
from typing import Dict, List, Tuple
from dataclasses import dataclass


@dataclass
class ContinuedFractionFormula:
    '''A discovered continued fraction formula.'''
    target_constant: str
    a_coeffs: List[int]  # partial denominators (polynomial in n)
    b_coeffs: List[int]  # partial numerators (polynomial in n)
    convergence_value: float
    error: float
    convergence_rate: float
    factorial_reduction_score: float
    depth: int


class RamanujanMachineSearcher:
    '''
    Search for polynomial continued fractions converging to zeta values.
    
    A polynomial continued fraction has the form:
        a_0 + b_1/(a_1 + b_2/(a_2 + ...))
    where a_n = polynomial in n, b_n = polynomial in n.
    
    The "factorial reduction" heuristic: for formulas converging to constants,
    gcd(p_n, q_n) is unusually large, causing the reduced convergents to grow
    much slower than (n!)^d.
    '''

    def __init__(self, max_degree: int = 2, max_coeff: int = 5, max_depth: int = 20):
        self.max_degree = max_degree
        self.max_coeff = max_coeff
        self.max_depth = max_depth
        self.results = []

    def _eval_polynomial(self, coeffs: List[int], n: int) -> int:
        '''Evaluate polynomial Σ coeffs[k] * n^k.'''
        return sum(c * (n ** k) for k, c in enumerate(coeffs))

    def _compute_convergents(self, a_coeffs: List[int], b_coeffs: List[int],
                             depth: int) -> Tuple[List[int], List[int], List[int]]:
        '''
        Compute convergents p_n/q_n for polynomial continued fraction.
        Returns (p_list, q_list, gcd_list).
        '''
        p = [0] * (depth + 2)
        q = [0] * (depth + 2)
        g = [0] * (depth + 2)

        # Initial conditions
        p[-1] = 1  # p_{-1} = 1
        p[0] = self._eval_polynomial(a_coeffs, 0)  # p_0 = a_0
        q[-1] = 0  # q_{-1} = 0
        q[0] = 1   # q_0 = 1

        for n in range(1, depth + 1):
            a_n = self._eval_polynomial(a_coeffs, n)
            b_n = self._eval_polynomial(b_coeffs, n)
            p[n] = a_n * p[n-1] + b_n * p[n-2]
            q[n] = a_n * q[n-1] + b_n * q[n-2]
            g[n] = self._gcd(abs(p[n]), abs(q[n]))

        return p[1:depth+1], q[1:depth+1], g[1:depth+1]

    def _gcd(self, a: int, b: int) -> int:
        while b:
            a, b = b, a % b
        return a

    def _factorial_reduction_score(self, p: List[int], q: List[int],
                                    g: List[int], depth: int) -> float:
        '''
        Measure how much gcd reduces the growth of convergents.
        For "good" formulas, reduced p_n/g_n and q_n/g_n grow like s^n
        (exponential), not like (n!)^d (super-exponential).
        
        Score = average log-growth-rate of reduced convergents.
        Lower = stronger factorial reduction.
        '''
        if depth < 5:
            return float('inf')

        growth_rates = []
        for n in range(3, depth):
            if g[n] > 1 and p[n] != 0 and q[n] != 0:
                reduced_p = abs(p[n]) // g[n]
                reduced_q = abs(q[n]) // g[n]
                if reduced_p > 0 and reduced_q > 0:
                    rate = (np.log(float(reduced_p)) + np.log(float(reduced_q))) / (2 * n)
                    growth_rates.append(rate)

        if len(growth_rates) < 3:
            return float('inf')

        return float(np.mean(growth_rates))

    def _test_convergence(self, a_coeffs: List[int], b_coeffs: List[int],
                          targets: Dict[str, float], depth: int = 20) -> List[ContinuedFractionFormula]:
        '''Test if the continued fraction converges to any target constant.'''
        p, q, g = self._compute_convergents(a_coeffs, b_coeffs, depth)

        if q[-1] == 0 or p[-1] == 0:
            return []

        frac_value = p[-1] / q[-1]
        frac_score = self._factorial_reduction_score(p, q, g, depth)

        # Compute convergence rate: how fast do convergents stabilize?
        if depth >= 10:
            late_values = [p[n] / q[n] for n in range(depth-5, depth) if q[n] != 0]
            if len(late_values) >= 3:
                convergence_rate = float(np.std(late_values))
            else:
                convergence_rate = 1.0
        else:
            convergence_rate = 1.0

        found = []
        for name, target in targets.items():
            if target == 0:
                continue
            error = abs(frac_value - target) / abs(target)
            if error < 0.01:  # Within 1%
                found.append(ContinuedFractionFormula(
                    target_constant=name,
                    a_coeffs=a_coeffs,
                    b_coeffs=b_coeffs,
                    convergence_value=frac_value,
                    error=error,
                    convergence_rate=convergence_rate,
                    factorial_reduction_score=frac_score,
                    depth=depth
                ))

        return found

    def search(self, targets: Dict[str, float] = None, n_candidates: int = 5000) -> Dict:
        '''
        Search for polynomial continued fractions converging to targets.
        
        Default targets include zeta values and related constants.
        '''
        if targets is None:
            targets = {
                'zeta_2': np.pi**2 / 6,
                'zeta_3': 1.202056903159594,
                'zeta_4': np.pi**4 / 90,
                'catalan': 0.915965594177219,
                'euler_mascheroni': 0.5772156649015329,
                'golden_ratio': (1 + np.sqrt(5)) / 2,
            }

        print(f"  [RamanujanMachine] Searching {n_candidates} candidates...")
        found_formulas = []
        best_scores = []

        # Generate candidates: simple polynomial forms
        # a_n = a0 + a1*n + a2*n^2, b_n = b0 + b1*n + b2*n^2
        count = 0
        for a0 in range(-self.max_coeff, self.max_coeff + 1):
            for a1 in range(-self.max_coeff, self.max_coeff + 1):
                for b0 in range(1, self.max_coeff + 1):  # b0 > 0 for convergence
                    for b1 in range(-self.max_coeff, self.max_coeff + 1):
                        if count >= n_candidates:
                            break
                        a_coeffs = [a0, a1]
                        b_coeffs = [b0, b1]

                        formulas = self._test_convergence(a_coeffs, b_coeffs, targets)
                        if formulas:
                            found_formulas.extend(formulas)

                        # Also compute score for all candidates
                        p, q, g = self._compute_convergents(a_coeffs, b_coeffs, 15)
                        score = self._factorial_reduction_score(p, q, g, 15)
                        if score < 5.0:  # Interesting reduction
                            best_scores.append((score, a_coeffs, b_coeffs))

                        count += 1
                    if count >= n_candidates:
                        break
                if count >= n_candidates:
                    break
            if count >= n_candidates:
                break

        # Sort by error
        found_formulas.sort(key=lambda f: f.error)

        self.results = {
            'strategy': 'ramanujan_machine_continued_fractions',
            'n_candidates': count,
            'n_found': len(found_formulas),
            'found_formulas': [
                {
                    'target': f.target_constant,
                    'a_coeffs': f.a_coeffs,
                    'b_coeffs': f.b_coeffs,
                    'value': f.convergence_value,
                    'error': f.error,
                    'rate': f.convergence_rate,
                    'frac_score': f.factorial_reduction_score,
                }
                for f in found_formulas[:10]
            ],
            'best_factorial_reduction': sorted(best_scores, key=lambda x: x[0])[:5],
        }
        return self.results

    def summary(self) -> str:
        r = self.results
        s = f"Ramanujan Machine Search\n{'='*50}\n"
        s += f"Candidates: {r['n_candidates']:,}\n"
        s += f"Formulas found: {r['n_found']}\n"
        for f in r['found_formulas'][:5]:
            s += f"  → {f['target']}: a={f['a_coeffs']}, b={f['b_coeffs']}, "
            s += f"error={f['error']:.6f}, value={f['value']:.8f}\n"
        if r['best_factorial_reduction']:
            s += f"Best factorial reduction score: {r['best_factorial_reduction'][0][0]:.4f}\n"
        return s