''' 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