vvs/sim-numerics
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Sign tests by product underflowed: FIXED everywhere (branch, not master - the paper is not yet updated). FA * fnew < 0.0 with FA ~ -1e-200, fnew ~ 1e-150 is -0.0 and '-0.0 < 0.0' is False, so on flat functions (x^19) the bracket update moved the wrong end, the root left the bracket and the run still reported convergence - the proposed rule returned 4.9e-9, bisection 3.5e-9, ITP 1.6e-8 with a claimed 1e-12 bracket; the 1e-6 harness scored them solved. Same defect in our bisection, false-position family, Ridders, Brent, ITP and twelve copies in stage scripts, simroot/sim.py and one test; all replaced by _opposite / _same_sign, pinned by test_sign_underflow.

Browse Source 
Impact, every sweep regenerated: the paper's key configurations are unchanged to the evaluation (proposed 2806/3220, credit 3568/3020, Chandrupatla, TOMS 748, Brent, bisection); 46 stated numbers move (ITP CHA 2857 -> 2896, ITP oracle 2611 -> 2647, record CHA figures by 14-18, the grid means by ~1%) and one finding weakens: the measured-order exponent of section 7.1 was cheaper in 29 of 60 far cells and 7 of 7 safeguards; with the fix 12 of 61 and 6 of 7, and on the wide-bracket set it now loses in every safeguard and drops 5 problems.

PREREG G1-G3 outcome on the fixed code: G2 (repeating step pattern, period from the data) PASSES Y1-Y4 and is ACCEPTED by the stop rule - proposed rule APS 2806 -> 2782 (below Chandrupatla's 2787), CHA 3220 -> 3169, both 6026 -> 5951, per problem vs Chandrupatla on APS 95/20/38 p = 8e-7. G1 fails Y1/Y2 (its Y3 miss was the sign bug), G3 fails Y2, A fails Y2. Regression 887 pass; the two red tests are the paper-number check (awaiting the paper update) and test_blend (fixed since, copy in stage16).

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
This commit is contained in:
vvs
2026-09-21 12:49:01 +01:00
co-authored by Claude Opus 5
parent 29ab49599c
commit 68df1d25db
34 changed files with 680 additions and 64 deletions
Download Patch File Download Diff File Expand all files Collapse all files
simroot/baselines.py
+21 -6
View File
@@ -47,6 +47,21 @@ class _Budget(Exception):
pass
def _opposite(u, v):
"""True when u and v have strictly opposite signs. NOT u * v < 0: for the tiny
values of a flat function near its root (x^19 there is ~1e-200) the product
underflows to -0.0, '-0.0 < 0.0' is False, and a bracket update then moves the
wrong end and loses the root while still reporting convergence (found 2026-09-21
on fun6.5, x^19 over [-10, 100])."""
return (u < 0.0 < v) or (v < 0.0 < u)
def _same_sign(u, v):
"""True when u and v are both strictly positive or both strictly negative; see
_opposite for why this is not u * v > 0."""
return (u > 0.0 and v > 0.0) or (u < 0.0 and v < 0.0)
def _converged(a, b, xtol):
return abs(b - a) <= xtol
@@ -68,7 +83,7 @@ def bisection(f, a, b, xtol, maxfev):
fx = g(x)
if fx == 0.0:
return x, g.n, True
if fa * fx < 0.0:
if _opposite(fa, fx):
b, fb = x, fx
else:
a, fa = x, fx
@@ -104,7 +119,7 @@ def _false_position_family(f, a, b, xtol, maxfev, correction):
fx = g(x)
if fx == 0.0:
return x, g.n, True
if fa * fx < 0.0:
if _opposite(fa, fx):
# b replaced; a became stale -> correct fa
fa = correction(fa, fb, fx)
b, fb = x, fx
@@ -170,9 +185,9 @@ def ridders(f, a, b, xtol, maxfev):
if fx == 0.0:
return x, g.n, True
# rebuild the tightest bracket among a, c, x, b
if fc * fx < 0.0:
if _opposite(fc, fx):
a, fa, b, fb = (c, fc, x, fx) if c < x else (x, fx, c, fc)
elif fa * fx < 0.0:
elif _opposite(fa, fx):
b, fb = x, fx
else:
a, fa = x, fx
@@ -202,7 +217,7 @@ def brent(f, a, b, xtol, maxfev, stats=None):
c, fc = a, fa
d = e = b - a
while True:
if fb * fc > 0.0:
if _same_sign(fb, fc):
c, fc = a, fa
d = e = b - a
if abs(fc) < abs(fb):
@@ -305,7 +320,7 @@ def itp(f, a, b, xtol, maxfev, kappa1=None, kappa2=2.0, n0=1):
fx = g(x)
if fx == 0.0:
return x, g.n, True
if fx * fa < 0.0:
if _opposite(fx, fa):
b, fb = x, fx
else:
a, fa = x, fx
simroot/pure.py
+79 -24
View File
@@ -93,7 +93,7 @@ from __future__ import annotations
import math
from .baselines import Counter, _Budget
from .baselines import Counter, _Budget, _opposite, _same_sign
# --------------------------------------------------------------------------
# [0,1] similarity functions of a normalised residual z: S(0) = 1 (identical
@@ -126,6 +126,63 @@ DIRECTIONS = ("distance", "dissimilar")
EXPONENTS = ("none", "order")
EXTRAPOLATE_MODES = ("none", "steps", "confirmed", "pattern", "steps+results")
def _ratio_sim(a, b):
"""The SIM ratio form for two magnitudes: min/max, 1 when they are equal."""
hi = max(a, b)
return min(a, b) / hi if hi > 0.0 else 0.0
def _extrapolation_s(mode, step, x, x_prev, x_prev2, fx, f_prev, f_prev2, moves):
"""How much of the geometric-series lengthening the history CONFIRMS: returns s in
[0, 1); the step becomes step/(1 - s). 0 means leave the step alone.
"steps" A: similarity of this step to the previous displacement.
"confirmed" G1: three moves the same way, both length ratios q1, q2 < 1; the
latest ratio q2, scaled by how similar q1 and q2 are.
"pattern" G2: A's s, applied only when the pattern of successive length ratios
repeats itself - period = the lag (2 .. half the history, so that two
repeats fit) of highest mean ratio similarity; applied iff that
similarity exceeds 1/2, the bisection constant.
"steps+results" G3: the smaller of A's s and the similarity of the last two RESULTS
f_n/f_(n-1) and f_(n-1)/f_(n-2); a result of either sign change
(a crossing) makes it 0.
No constant beyond bisection's 1/2 appears in any of them (PREREG 2026-09-21)."""
d1 = x - x_prev
if step == 0.0 or d1 == 0.0 or (step > 0.0) != (d1 > 0.0):
return 0.0
s_step = _ratio_sim(abs(step), abs(d1))
if mode == "steps":
return s_step
if mode == "confirmed":
if x_prev2 is None:
return 0.0
d2 = x_prev - x_prev2
if d2 == 0.0 or (d2 > 0.0) != (d1 > 0.0):
return 0.0
q1, q2 = abs(d1) / abs(d2), abs(step) / abs(d1)
if not (q1 < 1.0 and q2 < 1.0):
return 0.0
return q2 * _ratio_sim(q1, q2)
if mode == "pattern":
qs = [abs(moves[j]) / abs(moves[j - 1]) for j in range(1, len(moves))
if moves[j - 1] != 0.0 and moves[j] != 0.0]
k, best = len(qs), 0.0
for p in range(2, k // 2 + 1):
sims = [_ratio_sim(qs[j], qs[j - p]) for j in range(p, k)]
best = max(best, sum(sims) / len(sims))
return s_step if best > 0.5 else 0.0
if mode == "steps+results":
if f_prev2 is None or f_prev == 0.0 or f_prev2 == 0.0:
return 0.0
r_n, r_p = fx / f_prev, f_prev / f_prev2
if not (r_n > 0.0 and r_p > 0.0):
return 0.0
return min(s_step, _ratio_sim(r_n, r_p))
return 0.0
def _read_z(resid, sc, exponent, m):
"""The normalised residual the similarity function is applied to.
@@ -595,8 +652,9 @@ def pure_root(f, a, b, xtol, maxfev, similarity="ratio", scale="derived",
if scale not in ("derived", "derived3", "derived_rho", "derived3_rho",
"derived3_m", "derived3_cap", "start", "bracket", "far"):
raise ValueError(f"unknown scale {scale!r}")
if extrapolate not in ("none", "steps"):
if extrapolate not in EXTRAPOLATE_MODES:
raise ValueError(f"unknown extrapolate {extrapolate!r}")
moves = [] # displacement history, kept only for extrapolate="pattern"
if exponent not in EXPONENTS:
raise ValueError(f"unknown exponent {exponent!r}")
if exponent == "order" and scale not in ("derived3", "derived3_cap", "bracket", "far"):
@@ -641,7 +699,7 @@ def pure_root(f, a, b, xtol, maxfev, similarity="ratio", scale="derived",
return a, g.n, True
if fb == 0.0:
return b, g.n, True
if fa * fb > 0.0:
if _same_sign(fa, fb):
raise ValueError("f(a) and f(b) must straddle zero")
A, FA, B, FB = a, fa, b, fb # the proven bracket
@@ -866,26 +924,23 @@ def pure_root(f, a, b, xtol, maxfev, similarity="ratio", scale="derived",
step = _log_space_step(S, x, fx, far, slope, c_near, xtol)
if step is None:
step = (1.0 - S(z_now)) * (far - x)
if extrapolate == "steps":
# PREREG 2026-09-21 (A in the method). The similarity of THIS
# step to the previous displacement: two successive moves the
# same way, of similar length, mean the iterate creeps along a
# geometric series whose remainder is step*S/(1-S) - so the
# step is lengthened by 1/(1-S). On x^m under a Newton step
# S = 1 - 1/m and the factor is m, which lands on the root.
# At a simple root successive steps shrink superlinearly, S ->
# 0 and the step is untouched. Capped at the far end, the same
# bound the step itself has; no constant.
d_prev = x - x_prev
if (step != 0.0 and d_prev != 0.0
and (step > 0.0) == (d_prev > 0.0)):
s_steps = (min(abs(step), abs(d_prev))
/ max(abs(step), abs(d_prev)))
if s_steps < 1.0:
longer = step / (1.0 - s_steps)
reach = abs(far - x)
step = (longer if abs(longer) <= reach
else math.copysign(reach, step))
if extrapolate != "none":
# PREREG 2026-09-21 (A in the method) and PREREG 2026-09-21
# (gates G1-G3). Two successive moves the same way, of similar
# length, mean the iterate creeps along a geometric series
# whose remainder is step*S/(1-S): the step is lengthened by
# 1/(1-S). The modes differ only in how much of S the history
# CONFIRMS - see _extrapolation_s. Capped at the far end, the
# same bound the step itself has; no constant.
if extrapolate == "pattern":
moves.append(x - x_prev)
s_eff = _extrapolation_s(extrapolate, step, x, x_prev,
x_prev2, fx, f_prev, f_prev2, moves)
if 0.0 < s_eff < 1.0:
longer = step / (1.0 - s_eff)
reach = abs(far - x)
step = (longer if abs(longer) <= reach
else math.copysign(reach, step))
if stats is not None and "trace" in stats:
# DIAGNOSTIC HOOK, behaviour-neutral: records the state each
# step was taken from, so another model can be evaluated on the
@@ -934,7 +989,7 @@ def pure_root(f, a, b, xtol, maxfev, similarity="ratio", scale="derived",
if stats is not None:
stats['forced'] = forced; stats['pre_refused'] = pre_refused; stats['triple_changed'] = triple_changed; stats['dir_changed'] = dir_changed
return x_new, g.n, True
if FA * fnew < 0.0:
if _opposite(FA, fnew):
B, FB = x_new, fnew
else:
A, FA = x_new, fnew
simroot/sim.py
+4 -2
View File
@@ -40,6 +40,8 @@ no sigma, no temperature, no threshold, no tuned coefficient anywhere below.
"""
from __future__ import annotations
from .baselines import _opposite, _same_sign # 2026-09-21 fix
import math
from .baselines import Counter, _Budget
@@ -155,7 +157,7 @@ def sim_root(f, a, b, xtol, maxfev, similarity="neg_sq", prototypes="pair",
pending = False
if fx == 0.0:
return x, g.n, True
if fa * fx < 0.0:
if _opposite(fa, fx):
b, fb = x, fx
else:
a, fa = x, fx
@@ -228,7 +230,7 @@ def sim_root(f, a, b, xtol, maxfev, similarity="neg_sq", prototypes="pair",
mem.append((x, fx))
if fx == 0.0:
return x, g.n, True
if fa * fx < 0.0:
if _opposite(fa, fx):
b, fb = x, fx
else:
a, fa = x, fx