Add an arbitrary-precision generator for the angular grid tables - #110
Open
susilehtola wants to merge 17 commits into
Open
Add an arbitrary-precision generator for the angular grid tables#110susilehtola wants to merge 17 commits into
susilehtola wants to merge 17 commits into
Conversation
The tabulated Lebedev-Laikov and Delley grids carry about 16 significant
digits, so LebedevLaikov<T> cannot be more accurate than double however
wide T is. The format makes this worse than it looks: the tables are
unsuffixed floating literals inside std::array<T,N>, and an unsuffixed
literal has type double whatever T is, so writing more digits into the
existing format would change nothing at all.
Add a generator that refines the existing rules to any requested
precision.
Finding a spherical quadrature rule and refining one are different
problems, and only the second is needed here. Construction is a global
optimisation with many local minima and wants second-order globalisation
(a trust-region solver), but only at double precision, since the goal is
just to land in the right basin. Refinement starts from a rule already
correct to 16 digits, is purely local, and Newton converges quadratically
-- two or three iterations reach any precision, with no globalisation.
For the existing families construction is already done: the tables are
the result.
What makes it cheap is symmetry. Both families are octahedrally
invariant, so points fall into orbits carrying one weight and zero, one
or two angular parameters. The 5810-point Lebedev rule is a
385-parameter dense Newton solve rather than a 17430-parameter one.
The exactness conditions are imposed on even monomials in x and y only:
on the unit sphere z^2 = 1 - x^2 - y^2, so every even monomial reduces to
a fixed linear combination of those, and imposing them implies the rest
at about a seventh of the cost. Results are always verified against the
full redundant set, never the reduced one.
Emitted headers carry decimal strings rather than literals, converted
through detail::grid_scalar<T>::parse at construction. That works for the
built-in types as well as Boost.Multiprecision and MPFR, and costs one
parse per grid.
Measured on the 110-point rule integrating x^4 y^2 z^2:
current table regenerated
double 1.91e-15 1.74e-15
long double 1.26e-15 2.55e-19
cpp_bin_float_50 unreachable 1.94e-40
This commit adds the tool and the scalar-parsing helper only; no shipped
table is regenerated yet. Ahrens-Beylkin needs icosahedral orbit algebra
and Womersley needs a matrix-free least-squares step; both are described
in the README.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01FDTFYJMQ76iujDFNHzZyXF
Adds the Womersley path to the generator. These grids have no symmetry to
exploit -- every point is its own orbit -- so the unknowns are the 2N
tangential degrees of freedom of N points and the Jacobian is never
formed: CGLS needs only J*v and J^T*u.
Points are carried as unit 3-vectors with steps taken in each point's
tangent plane, not as (theta, phi). The spherical chart is singular at
the poles and real designs put points there: the 32-point design has one
at exactly z = 1, where d(p)/d(phi) vanishes. In that chart the Jacobian
acquires a structurally zero column that is a coordinate artifact rather
than a symmetry, and rotations about x and y are not expressible as
finite tangent vectors at all. Measured on the 32-point design:
zero cols near-null rotation generators
(theta, phi) 1 3 3.4e-02, 1.3e-16, 2.0e-16
tangent frame 0 3 2.0e-16, 1.4e-16, 2.9e-16
The remaining three-dimensional null space is the genuine one: any
rotation of a design is a design. It needs no special handling, because
CGLS started from a zero step keeps every iterate in range(J^T), which is
orthogonal to the null space, and so returns the minimum-norm -- already
rotation-projected -- step. Design.rotation_generators() exposes them
anyway, since a trust-region method used to *construct* a design does
need to project them out.
Two defects found by measurement rather than by reasoning, both recorded
in the source so they are not reintroduced:
* The tabulated points are unit vectors only to double precision, and
carrying 3-vectors inherits that inconsistency with |p| = 1. It
corrupts the Jacobian at the same 1e-16 level the refinement is trying
to work below -- 8.0e-16 relative error against finite differences,
1.7e-21 after renormalising -- and stalled the first Newton step. The
(theta, phi) chart had normalised implicitly.
* CG terminates in rank-many steps only in exact arithmetic; its rate
goes as the square root of the condition number. A cap of 2*2N starves
it, and the outer iteration then converges convincingly to a spurious
floor. On the 50-point design, capped at 200 iterations the residual
stalls at 2.45e-18, while at 600 the same first step reaches 2.16e-28.
The cap is now 12*2N with the residual tolerance doing the real
stopping and a stagnation guard for early exit.
With both fixed, against the full monomial set:
womersley_32 (order 7): 7.39e-16 -> 3.82e-30 -> 9.18e-41
womersley_50 (order 9): 5.75e-16 -> 2.16e-28 -> 4.96e-41
Conditions use monomials with c in {0,1}: on the sphere z^2 reduces,
leaving exactly (order+1)^2 of them -- the dimension of the polynomials
of degree <= order on S^2 -- rather than C(order+3,3). For order 7 that
is 64 rather than 120.
read_grid() also learns the second weight format: the Womersley headers
synthesise their equal weights with create_array<N,T>(4*M_PI/N) instead
of listing them.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01FDTFYJMQ76iujDFNHzZyXF
The AB grids are invariant under the icosahedral rotation group I (order
60) rather than the octahedral group, and without inversion -- which is
why they are not antipodally symmetric and why 15012 points divide as
15012/60 = 250.2 against 251 distinct weights.
The group is constructed numerically instead of from memorised
generators: it is exactly the set of rotations carrying the
icosahedron's vertices to themselves, and each is fixed by naming where
one vertex and one of its neighbours go. Verified: 60 elements, all
proper rotations, max |R^T R - I| = 3.9e-31.
Every tabulated AB grid decomposes cleanly:
AB-552 10 orbits 12 + 9x60 28 free parameters
AB-612 11 orbits 12 + 10x60 31 free parameters
AB-15012 251 orbits 12 + 250x60 751 free parameters
so even the largest is a 751-parameter problem rather than a
45036-parameter one.
Two tolerance defects fixed along the way, both the same shape -- a
tolerance tighter than the accuracy of the data it is applied to:
* read_grid()'s numeric pattern required an exponent. ahrens_beylkin_552
mixes plain decimals with exponent form inside one array, so the
parser silently read 48 points where there are 552. Any file in that
format was being misread.
* orbit() de-duplicated at 1e-18 while its input is tabulated to double
precision, so images that should coincide differ by ~1e-16 and a
12-point orbit came back as 60 near-duplicates. De-duplication is by
distance rather than by formatting coordinates, since a coordinate
that should be zero comes out as +-1e-31 and two such values have
different decimal representations.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01FDTFYJMQ76iujDFNHzZyXF
Adds the icosahedral refinement path, mirroring the octahedral one. Each
orbit carries one weight; a generic 60-point orbit carries two more
parameters as the position of its representative, stepped in that point's
tangent plane rather than through spherical angles. Because the orbit is
generated by the group action, the derivative of an orbit point with
respect to a tangent step of its representative is just the group element
applied to that step.
AB-72 order 14 2 orbits 4 params 1.75e-15 -> 1.58e-30
AB-132 order 19 3 orbits 7 params 4.34e-16 -> 7.89e-31
AB-192 order 23 4 orbits 10 params 7.66e-16 -> 1.58e-30
two Newton steps each, from the tabulated double-precision data.
Two defects found and fixed while getting there:
* decompose() used the rounded weight as both the grouping key and the
weight value, truncating it to the key's 12 digits. It showed up as
1.7e-11 in the sum of weights while the points matched to 4.7e-33.
* A 12-, 20- or 30-point orbit has no free parameter -- its position is
fixed by the symmetry -- but taking the representative from tabulated
double-precision data makes it a *generic* point whose 60 images
cluster in fives rather than coinciding. De-duplication then keeps an
arbitrary member of each cluster and injects ~1e-16 of asymmetry into a
grid that should be exactly symmetric. Special orbits are now snapped
onto the exact vertex, face-centre or edge-midpoint position, rotated
onto the tabulated one.
The second was masked by the residual: summing a monomial over all 60
group images annihilates odd moments identically, whatever the
representative, so the odd-moment conditions carry no information about
the fit and could not see the asymmetry. Verification against the
materialised grid does see it, and now agrees with the residual.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01FDTFYJMQ76iujDFNHzZyXF
Three changes, all needed before the larger Ahrens-Beylkin grids are reachable. Condition selection. The full set is (order+1)^2 monomials -- 1600 at order 39 -- against a few dozen unknowns, since the icosahedral symmetry makes most of them dependent. Candidates are now visited stratified across degree and kept when their Jacobian row increases the column rank, stopping early once the rank is full. AB-192 drops from 576 conditions to 30. Results are always verified against the full set, so a bad selection surfaces as a failed verification rather than a silently wrong grid. Image cache hoisted out of residual_and_jacobian. Selection asks for many single-row Jacobians, and rebuilding the 60 group images per row dominated: AB-192 selection went from 27s to 9s, and its refinement from 45s to 1s. Rank-aware step. The Jacobian is not always full rank -- AB-312 comes out 13 of 16 -- and both a normal-equations LU and mpmath's qr_solve refuse a singular matrix outright. The step is now solved only in an independent set of columns, found by modified Gram-Schmidt with a relative threshold, leaving the unresolved directions at zero. Same lesson as the Womersley path, where CGLS gave the minimum-norm step for free. Also records a negative result: refining AB-552 from its tabulated data diverges, 0.4619 -> 0.03831 -> 145.9. That is expected rather than surprising -- its points are wrong, not merely imprecise, and no weighting of them integrates to order 39 -- but it settles that the corrupt table cannot be recovered by Newton refinement. It needs globalisation. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01FDTFYJMQ76iujDFNHzZyXF
generate.py now dispatches on the family: octahedral refinement for
Lebedev-Laikov and Delley, icosahedral for Ahrens-Beylkin, and the
equal-weight design path for Womersley. Each verifies against the full
condition set before emitting, and refuses to write a table that did not
reach the requested precision.
ahrens_beylkin_72: order 14, 2 icosahedral orbits, 4 free parameters
selected 12 of 225 conditions, rank 4/4
1.754e-15 -> 2.091e-33
verified 2.431e-32 against the full condition set
The emitted load() is now templated on its containers, so it fills the
library's std::vector<cartesian_pt_t<T>> and std::vector<T> directly --
the same role detail::copy_grid plays for the literal tables, which is
what the eventual switch needs. Verified by filling those vectors from a
generated header.
Two smaller fixes: the emit namespace and title tables only knew the two
octahedral families, so an Ahrens-Beylkin header came out in namespace
`ahrens_beylkin` rather than `AhrensBeylkinGrids`; and condition
selection now caps its walk, since a genuinely rank-deficient Jacobian
means the rank never fills and the search would otherwise visit all 1600
candidates. An incomplete basis is harmless -- the step solver handles
deficiency, and verification against the full set is the real guard.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01FDTFYJMQ76iujDFNHzZyXF
batch.py regenerates every tabulated size of a family, reading the sizes
from the family's own algebraic_order_by_npts table so the two cannot
drift apart. Sizes are independent, so they run in separate processes. A
size that fails to reach the requested precision is reported and skipped
rather than written, so a partial run never emits an unverified table.
delley: 6 sizes, 30 digits, 4 jobs
ok delley_14 err 2.552e-51
ok delley_26 err 4.786e-51
ok delley_50 err 5.264e-51
ok delley_110 err 1.246e-49
ok delley_194 err 6.018e-49
ok delley_302 err 2.163e-48
switch_to_load.py points a family's dispatch header at the regenerated
tables:
detail::copy_grid<delley_50<RealType>>(points, weights)
-> delley_50<RealType>::load(points, weights)
and drops the `weights[i] *= 4.0*M_PI` loop for the two families whose
shipped tables are normalised to one -- a regenerated table already
carries the 4*pi convention. The loop index is `auto` in delley.hpp and
`size_t` in lebedev_laikov.hpp, so the pattern does not pin the type.
Dry runs cover all four families: 31 calls in lebedev_laikov, 18 in
delley, 56 in ahrens_beylkin, 125 in womersley.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01FDTFYJMQ76iujDFNHzZyXF
Records what is known about the one table that cannot be regenerated by
refinement, with AB-612 as a control throughout.
The decisive result is that the AB-552 points admit no valid weighting.
The exactness conditions are linear in the weights, so the best possible
weights for a fixed point set follow from one least-squares solve:
AB-612 (known good) 8.3e-14, all positive
AB-552 2.4e-02, five negative weights
The positions are wrong, not just the weights, and not merely imprecise.
Refinement behaves accordingly: Gauss-Newton diverges (0.4619 -> 0.03831
-> 145.9), and Levenberg-Marquardt globalises the step and makes real
progress (0.4619 -> 9.99e-4 -> 7.75e-6 -> 1.19e-7) but plateaus five
orders short, on a condition subset that does not determine all 28
parameters.
Also records what is *not* evidence, since I nearly reported it as such:
the Jacobian's rank over low-degree conditions is 6 against 28
parameters, which looks alarming until AB-612 is measured the same way
and gives rank 6 against 31. Few icosahedral invariants exist at low
degree, so that number says nothing about either grid.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01FDTFYJMQ76iujDFNHzZyXF
Pointing --out at the family's own data directory destroys the input. Generation reads the shipped literal table for each size, so a partial run replaces the tables the remaining sizes still need, and every later size then fails with "substring not found" -- which is what happened, on a scratch worktree, when I did exactly this. Generate into a separate directory and install the results deliberately. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01FDTFYJMQ76iujDFNHzZyXF
The octahedral refinement still solved its normal equations with a plain lu_solve, and mpmath refuses a singular matrix outright. That surfaced as delley_1454 failing with "matrix is numerically singular" at 20 digits while succeeding at 40 -- so it is a robustness property of the solver rather than of any particular grid, and it would have bitten whichever size happened to be marginal at the chosen precision. Move independent_columns() and the step solver into linalg.py and use them from both paths. delley_1454 now regenerates at 20 digits, verified to 1.287e-24 against the full condition set for order 65. Regression suite 11/11. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01FDTFYJMQ76iujDFNHzZyXF
The rank-aware step I added earlier used a fixed 1e-16 threshold on the
normal equations, and that discarded columns which were merely
ill-conditioned rather than null. delley_974 lost its convergence as a
result -- 3.11e-29 instead of 6.93e-51, failing the generator's own
precision check. The batch caught it; the unit tests did not.
Chasing the threshold does not work, because the two failure modes trade
off against each other:
threshold delley_974 @ 40 delley_1454 @ 20
fixed 1e-16 stalls 3.1e-29 solves
precision-tied 6.93e-51 "numerically singular"
QR, 1e-14 cut 11 iters, 2e-33 diverges, 1.7e+146
QR, precision-tied 6.93e-51, 2 iters needs more precision
The tension is not in the threshold. It is that J^T J squares the
condition number. The step is now a rank-revealing modified Gram-Schmidt
QR of J itself, back-substituted, never forming the normal equations.
Checked against a synthetic rank-deficient system where it finds rank 2
and leaves a 1e-31 residual.
The threshold that remains cuts only what is unresolvable at the working
precision, so it discards genuinely null directions and keeps everything
else. delley_1454 at 20 digits is then simply a case that needs more
precision than it was given -- 102 parameters at that conditioning are
not resolvable in 20 digits by any of these methods, and it refines
cleanly at 40. The generator's final verification catches it and refuses
to write the table, which is the right outcome.
Regression suite 11/11.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01FDTFYJMQ76iujDFNHzZyXF
Records why the step is a rank-revealing QR of the Jacobian rather than a solve of the normal equations, and that a size needing more working precision than it is given fails verification instead of being written. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01FDTFYJMQ76iujDFNHzZyXF
The octahedral orbit types carry implicit constraints. bk is
(l, l, sqrt(1-2 l^2)) and is only real for l <= 1/sqrt(2); ck and dk have
their own. A Newton step can overshoot one, and mpmath then returns a
complex square root:
l = 0.30000000 sqrt(1-2 l^2) = 0.90553851
l = 0.70710678 sqrt(1-2 l^2) = 5.7931539e-5
l = 0.72000000 sqrt(1-2 l^2) = (0.0 + 0.19183326j)
The complex value then propagates silently until something tries to order
two of them, which is where it finally surfaced -- as
"TypeError: '<=' not supported between instances of 'mpc' and 'mpc'",
far from the cause. It took out delley_2354, 2702, 3074 and 3470, the
four largest sizes in the family.
Halve the step until every orbit base is real again. Backtracking is
cheaper than reparameterising the orbits, and it leaves the well-behaved
sizes untouched.
Found by the batch, not by the test suite: the eleven regression cases
use small grids whose parameters sit well inside their domains, and only
the 2354-and-larger sizes push one to the boundary.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01FDTFYJMQ76iujDFNHzZyXF
Backtracking only on the orbits' domain was not enough. It stops the complex-square-root failure, but on the larger Delley grids the undamped step is simply bad: delley_2354 (order 83, 64 orbits, 161 parameters) walks away from a starting residual of 1.1e-16 to 6.4, 6.2, 5.6 over successive iterations. Halve the step until it both keeps every orbit real and reduces the residual, and give up rather than accept an increase. This is ordinary damped Newton; the quadratic convergence on the well-behaved sizes is unaffected, since the full step is accepted on the first try there. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01FDTFYJMQ76iujDFNHzZyXF
A weight column and an angular column differ by orders of magnitude, and that disparity alone wrecks the conditioning of the least-squares solve. On delley_2354 (order 83, 161 parameters) the undamped step is unusable, and once damped the iteration crawls -- a factor of two per step instead of squaring. Scale every column to unit norm before the QR and unscale the step afterwards. It costs one pass over the Jacobian and is the difference between converging and not. Regression suite 11/11. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01FDTFYJMQ76iujDFNHzZyXF
The guard digits were a fixed margin, and that is wrong because the conditioning varies enormously with grid size. Measured after column equilibration, delley_1730's Jacobian has a condition number of 2.0e+24 -- so 24 of the working digits are gone before the step is even computed, and the larger sizes are worse. That is what has been defeating the four biggest Delley grids: a bad step from an under-resourced solve, which then overshot an orbit's domain and surfaced as a complex square root. Rather than predict the conditioning, try a modest guard and retry with more when the final verification falls short. The well-conditioned sizes pay nothing, since they pass on the first attempt. This is the root cause behind the last three commits. Each of those was a correct fix for a real symptom -- complex roots, divergence, crawling convergence -- but all three were downstream of a step computed with insufficient precision for the conditioning. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01FDTFYJMQ76iujDFNHzZyXF
Retrying on failure works but wastes a long attempt first, and on the
biggest Delley grids that attempt takes many minutes. The conditioning is
predictable enough to start from a sensible guard: measured on the
equilibrated Jacobian it is 2.0e+24 at 120 parameters (delley_1730) and
8.19e+26 at 140 (delley_2030), so log10(cond) is roughly 0.15 * params.
The extra guard this gives:
delley_110 0 delley_2030 24
delley_974 10 delley_2354 28
delley_1730 20 delley_3470 42
Small grids are unaffected and pay nothing. The retry ladder stays as the
backstop for anything the estimate underestimates.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01FDTFYJMQ76iujDFNHzZyXF
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Adds
tools/angular_grids, an offline generator that regenerates thetabulated angular grids at arbitrary precision, plus the small runtime
helper the regenerated tables need.
Why
The shipped tables carry about 16 significant digits, so
LebedevLaikov<T>cannot be more accurate than
doublehowever wideTis. The format makesthis worse than it looks: the tables are unsuffixed floating literals inside
std::array<T, N>, and an unsuffixed literal has typedoublewhateverTis. Writing more digits into the existing format would change nothing.
Measured on the 110-point rule, integrating
x^4 y^2 z^2:doublelong doublecpp_bin_float_50Approach
Finding a spherical quadrature rule and refining one are different problems,
and only the second is needed here. Construction is a global optimisation
with many local minima and wants a trust-region method, but only at double
precision. Refinement starts from a rule already correct to 16 digits, is
purely local, and Newton converges quadratically. For every family already
in the library, construction is done: the tables are its output.
Tractability comes from symmetry. Lebedev-Laikov and Delley are octahedrally
invariant; Ahrens-Beylkin is invariant under the icosahedral rotation group
(order 60, no inversion -- which is why those grids are not antipodally
symmetric, and why 15012 points divide as 15012/60 = 250.2 against 251
distinct weights). Womersley grids are equal-weight spherical designs with no
symmetry at all, handled matrix-free.
So the largest Lebedev rule is a 385-parameter solve rather than a
17430-parameter one.
Output format
Regenerated headers store the decimal digits as strings and convert them
through
detail::grid_scalar<T>::parseat grid construction, which works forthe built-in types as well as Boost.Multiprecision and MPFR, and costs one
parse per grid.
load()is templated on its containers, so it fills thelibrary's
std::vector<cartesian_pt_t<T>>andstd::vector<T>in exactlythe role
detail::copy_gridplays today.Verified
switch_to_loadapplied to the Delley dispatch: 18 calls rewritten, the4.0*M_PInormalisation dropped, 9/9 library tests pass.the full redundant condition set rather than the reduced one used to fit.
only library-side addition is
util/grid_scalar.hpp.Nothing shipped is changed by this PR -- the generator is offline and no
table is replaced here. Regenerating and switching a family is a deliberate
follow-up, one family at a time.
Numerics worth reviewing
The larger grids are badly conditioned, and getting them to converge took
several passes:
Jacobian, never the normal equations.
J^T Jsquares the condition number,and with a threshold on the normal equations the failures trade off against
each other: delley_974 stalls at 3.1e-29 instead of 6.9e-51, while
tightening the threshold puts delley_1454 back to "numerically singular".
differ by orders of magnitude, and that disparity alone wrecks the solve.
equilibrated Jacobian, the condition number is 2.0e+24 at 120 parameters
(delley_1730) and 8.19e+26 at 140 (delley_2030). A fixed guard is simply the
wrong model; the retry ladder is the backstop.
bkis(l, l, sqrt(1-2 l^2))and only real forl <= 1/sqrt(2); an undamped stepovershoots, mpmath returns a complex root, and it surfaces far away as a
TypeErrorabout ordering complex numbers.Every one of those was found by running the full batch, not by the regression
suite, which uses small well-conditioned grids. If this lands, the suite would
benefit from at least one large size per family even at the cost of runtime.
Not covered
The 552-point Ahrens-Beylkin table (#108) cannot be regenerated by refinement.
Its points admit no valid weighting at all: the conditions are linear in the
weights, and a least-squares solve over the full order-39 set gives 2.4e-02
with five negative weights, against 8.3e-14 for AB-612 as a control. So the
positions are wrong rather than imprecise. Gauss-Newton diverges from them and
Levenberg-Marquardt plateaus five orders short. That size needs the original
Ahrens-Beylkin data, or to stay withdrawn as #101 leaves it.
🤖 Generated with Claude Code