Upload problem_solvers/ramanujan_machine.py

problem_solvers/ramanujan_machine.py

@@ -0,0 +1,241 @@
+'''
+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