"""Reversible finite-state wave automata with exact rational force evaluation. All arithmetic used by an automaton step is performed with Python integers. Floating point is used only to sample prescribed smooth initial data and to analyse already computed trajectories. Rounding ties are away from zero. """ from __future__ import annotations from dataclasses import dataclass from fractions import Fraction from functools import reduce from math import gcd import numpy as np def as_fraction(value): return value if isinstance(value, Fraction) else Fraction(str(value)) def lcm(a, b): return abs(a * b) // gcd(a, b) def round_ratio(numerator: int, denominator: int) -> int: """Odd nearest-integer rounding of an exact rational; ties away from zero.""" if denominator <= 0: raise ValueError('positive denominator required') numerator = int(numerator) sign = 1 if numerator >= 0 else -1 return sign * ((2 * abs(numerator) + denominator) // (2 * denominator)) def _round_array(numerator, denominator): numerator = np.asarray(numerator, dtype=object) return np.sign(numerator) * ((2 * np.abs(numerator) + denominator) // (2 * denominator)) def centered(values, modulus): half = modulus // 2 return (np.asarray(values, dtype=object) + half) % modulus - half def _common_integer_coefficients(coefficients): denominator = reduce(lcm, (x.denominator for x in coefficients), 1) return tuple(int(x * denominator) for x in coefficients), denominator @dataclass class Model: sites: int length: Fraction courant: Fraction = Fraction(1, 2) delta_factor: int = 64 mass_phi_sq: Fraction = Fraction(0) mass_chi_sq: Fraction = Fraction(0) eta: Fraction = Fraction(0) lambda_phi: Fraction = Fraction(0) lambda_chi: Fraction = Fraction(0) g: Fraction = Fraction(0) species: int = 1 modulus: int = (1 << 61) - 1 delta_override: Fraction | None = None def __post_init__(self): for name in ('length', 'courant', 'mass_phi_sq', 'mass_chi_sq', 'eta', 'lambda_phi', 'lambda_chi', 'g'): setattr(self, name, as_fraction(getattr(self, name))) if self.sites < 3 or self.species not in (1, 2): raise ValueError('sites >= 3 and one or two species required') if self.modulus < 3 or self.modulus % 2 == 0: raise ValueError('odd modulus >= 3 required') self.h = self.length / self.sites self.tau = self.courant * self.h self.delta = as_fraction(self.delta_override) if self.delta_override is not None else self.h**4 / self.delta_factor self.nu_sq = self.courant**2 t2 = self.tau**2 d2 = self.delta**2 masses = (self.mass_phi_sq, self.mass_chi_sq) lambdas = (self.lambda_phi, self.lambda_chi) self.onsite = [] for a in range(self.species): coeffs = (t2*masses[a], t2*self.eta if self.species == 2 else Fraction(0), t2*lambdas[a]*d2, t2*self.g*d2 if self.species == 2 else Fraction(0)) self.onsite.append(_common_integer_coefficients(coeffs)) def to_dict(self): result = {name: str(getattr(self, name)) for name in ( 'length','courant','h','tau','delta','mass_phi_sq','mass_chi_sq','eta', 'lambda_phi','lambda_chi','g')} result.update(sites=self.sites, species=self.species, modulus=self.modulus, delta_factor=self.delta_factor) mass_matrix=np.array([[float(self.mass_phi_sq+self.eta),-float(self.eta)], [-float(self.eta),float(self.mass_chi_sq+self.eta)]]) maximum_mass=(float(np.linalg.eigvalsh(mass_matrix).max()) if self.species==2 else float(self.mass_phi_sq)) result['linear_cfl_upper_bound']=4*float(self.courant)**2+float(self.tau)**2*maximum_mass result['linear_cfl_strict_margin']=4-result['linear_cfl_upper_bound'] return result def force(self, fields): q = np.asarray(fields, dtype=object) if q.shape != (self.species, self.sites): raise ValueError(f'expected {(self.species, self.sites)}, got {q.shape}') result = np.empty(q.shape, dtype=object) numerator, denominator = self.nu_sq.numerator, self.nu_sq.denominator for a in range(self.species): current = q[a] other = q[1-a] if self.species == 2 else current * 0 spatial = _round_array(numerator*(np.roll(current, -1)-current), denominator) spatial += _round_array(numerator*(np.roll(current, 1)-current), denominator) (cm, ce, cl, cg), den = self.onsite[a] onsite = cm*current + ce*(current-other) if cl: onsite = onsite + cl*current*current*current if cg: onsite = onsite + cg*current*other*other result[a] = spatial - _round_array(onsite, den) return result def step(self, fields, momenta, modular=True): q = np.asarray(fields, dtype=object) p = np.asarray(momenta, dtype=object) if modular: q = centered(q, self.modulus) p = centered(p, self.modulus) p_new = p + self.force(q) q_new = q + p_new if modular: p_new = centered(p_new, self.modulus) q_new = centered(q_new, self.modulus) return q_new, p_new def inverse(self, fields_new, momenta_new, modular=True): qp = np.asarray(fields_new, dtype=object) pp = np.asarray(momenta_new, dtype=object) q = qp - pp if modular: q = centered(q, self.modulus) p = pp - self.force(q) if modular: p = centered(p, self.modulus) return q, p def continuum_force(self, fields): q = np.asarray(fields, dtype=float) out = np.zeros_like(q) out[0] = float(self.mass_phi_sq)*q[0] + float(self.lambda_phi)*q[0]**3 if self.species == 2: out[0] += float(self.eta)*(q[0]-q[1]) + float(self.g)*q[0]*q[1]**2 out[1] = (float(self.mass_chi_sq)*q[1] + float(self.eta)*(q[1]-q[0]) + float(self.lambda_chi)*q[1]**3 + float(self.g)*q[1]*q[0]**2) return out def _sample_round(array, delta): scaled = np.asarray(array, dtype=float)/float(delta) if not np.all(np.isfinite(scaled)) or np.max(np.abs(scaled)) >= 2**52: raise ValueError('initial sampling requires finite scaled values below 2**52') rounded = np.sign(scaled)*np.floor(np.abs(scaled)+0.5) return rounded.astype(np.int64).astype(object) def initialize(model, displacement, velocity, laplacian): """Second-order backward sample, then exact memory-register encoding.""" displacement = np.asarray(displacement, dtype=float).reshape(model.species, model.sites) velocity = np.asarray(velocity, dtype=float).reshape(model.species, model.sites) laplacian = np.asarray(laplacian, dtype=float).reshape(model.species, model.sites) tau = float(model.tau) acceleration = laplacian-model.continuum_force(displacement) q = _sample_round(displacement, model.delta) q_prev = _sample_round(displacement-tau*velocity+0.5*tau*tau*acceleration, model.delta) return q, q-q_prev def _safe_int64(array): lo, hi = -(1 << 63), (1 << 63)-1 if any(int(x) < lo or int(x) > hi for x in np.asarray(array, dtype=object).flat): raise OverflowError('saved register does not fit int64; no silent conversion permitted') return np.asarray(array, dtype=np.int64) def evolve(model, fields, momenta, steps, paired=True): """Compute and retain every lifted step and an independent modular update.""" q, p = np.asarray(fields,dtype=object).copy(), np.asarray(momenta,dtype=object).copy() qm, pm = centered(q,model.modulus), centered(p,model.modulus) shape = (steps+1,model.species,model.sites) q_history, p_history = np.empty(shape,dtype=np.int64), np.empty(shape,dtype=np.int64) maximum = 0 consistent, nowrap = True, True for n in range(steps+1): q_history[n], p_history[n] = _safe_int64(q), _safe_int64(p) maximum = max(maximum, max(abs(int(x)) for x in q.flat), max(abs(int(x)) for x in p.flat)) nowrap = nowrap and 2*maximum < model.modulus if paired: consistent = consistent and np.array_equal(centered(q,model.modulus),qm) and np.array_equal(centered(p,model.modulus),pm) if n == steps: break q,p=model.step(q,p,modular=False) if paired: qm,pm=model.step(qm,pm,modular=True) if not consistent: raise AssertionError('modular and lifted trajectories disagree') if not nowrap: raise AssertionError('trajectory left its centered no-wrap interval') physical = float(model.delta)*q_history.astype(float) return dict(Q=q_history,P=p_history,fields=physical, times=np.arange(steps+1)*float(model.tau)), dict( paired_consistent=bool(consistent),nowrap_pass=bool(nowrap),steps=int(steps), max_abs_register=int(maximum),modulus=int(model.modulus), physical_register_bound=float(maximum*model.delta), physical_half_modulus=float(model.modulus*model.delta/2), integer_margin=int(model.modulus//2-maximum),all_steps_saved=True, arithmetic='Python arbitrary-precision integer; exact rational nearest-integer force')