Repository navigation
4.0.0: every function rounded once, with special functions, enclosures and directed rounding - #16
Merged
Merged
Conversation
The accuracy tool gains a second sweep. It runs every function on 104 layouts from I8F0 to I3F125: every value of each 8- and 16-bit one, and 40,000 seeded samples of each wider one, crowded round each function's hard points. Inputs reach MPFR exactly, from raw bits, and the error is counted in ulps of the type rather than through f64. The same inputs check odd and even symmetry bit for bit, and monotonicity over sorted inputs and their one- and two-ulp neighbours. worst.txt holds every row, the README a table for I16F16, I32F32 and I64F64, and CI fails a row that gets worse than on main. On today's code the sweep finds: - sin, cos and tan wrong by up to 0.76 near ±π on types with three integer bits, where 2π overflows (debug builds panic there instead); - argument reduction through a rounded 2π, 4,300 ulps off for I16F16 at 32,766 and 2.2e17 for I64F64 at 2^63; - asin and acos losing about √ulp near ±1, 93,000 ulps for I32F32; - symmetry failing on most types for 14 odd or even functions, since two's-complement shifts and multiplies round toward −∞; - atan, ln, log2, log10, tanh and the inverse hyperbolics stepping backwards on most types, about one adjacent pair in 500; - errors near 2^62 ulps on every type with more than 63 fractional bits, from the I1F63 tables. tests/golden.txt pins a hash of each function's outputs on four types, and a second test wants CHANGELOG.md to name every function whose pinned bits change, under the version that changed them.
Each function now converts its arguments exactly into a binary float with a 64-, 128- or 192-bit mantissa, twice its type's width or more, evaluates there and rounds once to the type, so every type from I8F0 to I3F125 gets the same treatment. Sign and magnitude are apart, so odd and even functions are exactly so. sin and cos reduce through 384 bits of 1/(2π), exactly at any magnitude, and every series runs a fixed number of terms around 384-bit tables, summed in 62-, 126- or 190-bit fixed point. asin and acos take √((1 − a)(1 + a)), and the hyperbolics and their inverses go through expm1 and log1p, so none cancels. CordicNumber is sealed, since each type now carries its working arithmetic. On the accuracy tool's 104 layouts every function is within half an ulp of MPFR wherever it looked: every value of the 8- and 16-bit types and 40,000 samples of each wider one. Worst cases were up to 2e21 ulps. Symmetry and monotonicity now hold everywhere, saturation at MIN aside, and the tool fails any function at or past ERROR_ULPS, one ulp, or that breaks either. New functions round an exact result once: mul_add and mul_div through a 256-bit product, dot through a 384-bit sum, div_or_saturate, which never panics, and sqrt_or_zero. The price is speed: on the benchmarks I16F16 functions take 2.4 to 11 times as long as before and I64F64 ones 1.3 to 9 times, sqrt aside. Every function's mean error against MPFR falls, by up to 20 times for sin on I16F16, and its worst case falls to half an ulp on every layout. The golden pins gain I3F5 and I3F61, which 3.1.1 could not run in debug builds, and the arithmetic functions; CHANGELOG.md names each function whose bits moved, which is all but sqrt. The version moves to 4.0.0, which the rest of this work joins unreleased.
enclose turns any result into the interval of one ulp each way that holds its exact value, by ERROR_ULPS. sin_cos_range finds the extremes inside an interval from the exact turn fraction of its start, counting one within the working error as inside, and exp_range, atan_range and tanh_range take their rising ends; all round outward. A property test samples 34,000 points of random intervals and finds each inside its range. hypot, norm and normalize form their squares in the working arithmetic, so a vector of a few ulps keeps its exact length and direction. scale_exponent and scale_by scale by powers of two exactly. sqrt_product, rsqrt and asin_from_gap, which takes 1 − |x| in place of x, round once. On all 104 layouts each is within half an ulp of MPFR, and no other output moves.
…ions expm1, exprel, exprel2 and exp_moments come from φₖ(z) = (e^z − Σ zⁿ/n!)/z^k, summed by its series within ½ of 0 so it keeps relative precision there; Eₖ(x) = (k − 1)!·e^(−x)·φₖ(x). log1p_residual writes ln(1 + v) − v as two terms of order v² that do not cancel, and log1p_ratio and log1p_residual_ratio divide them through. stumpff gives c₀ to c₃ of either sign without a branch at 0: the series within ¼ of it, past that the closed forms in √|z|, with c₂ through the half angle. tan_ratio, atan_ratio, atanh_ratio and tanh_ratio keep their bits near 0, and sin and cos of a working value now reduce exactly from its own mantissa. On all 104 layouts each is within half an ulp of MPFR at 960 bits, and no other output moves.
A Modulus holds m or its complement k' exact as given, and the working values of m, k' and k'² all come from it, so no pair can disagree. Carlson's R_F, R_J, R_C and R_D run by duplication to a fixed depth per width, closed by DLMF 19.36.1 or Carlson's 1995 series, deep enough for arguments hundreds of octaves apart. R_C and the R_J terms take 1 + e = 2√p(p + λ)/d with e, so neither cancels near −1, and R_J trades a p far above the others for a small q by DLMF 19.21.12, where duplication would need over a hundred steps. K comes from the AGM, E and Π through Carlson. sn, cn, dn and am follow the AGM down from φ_N with each arcsine taken from its exact gap to 1, reflect past K/2, and take dn = √(cn² + k'²sn²); the separatrix has closed forms. jacobi_inverse and ellip_pi take (sn, cn, dn) together and extend past π/2. Against independent 512-bit references every function is within half an ulp on all 104 layouts, sampled 3,000 times each, and matches mpmath to 1e-15 on I64F64.
sin_cos_turns reduces a fraction of a turn by its own bits, so no rounded π enters at any magnitude; an I64F64 value holds the same bits as a U64F64 fraction below one half. atan2_turns and div_by_tau scale by 1/(2π) before their one rounding. wrap_radians and wrap_principal take the exact turn fraction, and where the nearest value of the type falls outside [0, 2π) or (−π, π] they take the end of the range nearer modulo 2π. On all 104 layouts each is within its bound against MPFR, the wraps' error taken modulo 2π: half an ulp, but 0.98 for wrap_principal at ±π, where the nearest value to π can lie outside the range.
…ction Complex<T> holds a value of any of the types in each part. Multiplication and division form each part from exact sums of products, through a new dot2 in the working arithmetic, so a part rounds once however far its products cancel, and no division by a non-zero divisor overflows. exp, exprel and exprel2 follow, then the Faddeeva function, the complex erf, E₁ and Γ(a, z) for a in (0, 128], on the principal branch; erf, erfc, fresnel and gamma take a real argument. e^(−z²) takes y² − x² and the phase 2xy from exact products, the phase reduced by 576 bits of 1/(2π), so 128-bit arguments keep their phase. E₁ and Γ(a, z) choose between series, asymptotic expansion and Lentz's fraction by r + Re z, the power of e by which the series' terms outgrow their sum. Γ(a, z) takes its 1/a parts off exactly for a ≤ ½, z^a reduces e^(iπa) exactly left of the imaginary axis, and Re erf near that axis has its own Taylor series, so small parts beside huge ones keep their precision. beyond() rises to 2^8192, past e^HUGE; no output changes. On all 104 layouts each part is within half an ulp of MPFR and MPC references, which the harness gains along with a --quick mode that CI runs and holds to the bounds alone. Comments across the crate are trimmed so comment density is back below main's: 5.78% of lines against 5.79%, and 10.50% of characters against 11.48%.
…E₄, E₅ exp_moments now returns E₁ to E₅ with e^(−x), all from one exponential: near 0 the series of φ₅ and φₖ = 1/k! + z·φₖ₊₁, past it expm1 and the forward recurrence, which loses at most 11 bits. Complex::exprel_weighted gives exprel and exprel2 under the weight ((1 − εt)/(1 − εs))^p/(1 − εs), the kernels of a burn whose mass falls by ε against linear drag. Past |z| = 4(p + 1 + PREC)·δ, δ = ε/(1 − ε), they are sums by parts with the weight's derivatives in closed form; else the integral runs from each end toward the least of its log-convex modulus, in steps that take e^(ct) whole through exact moments ∫₀¹ sⁿ·e^(ws) ds and only the weight's curvature as a series, so neither |z| nor the drag sets the step count. Each pass stops where its integrand times what is left no longer moves the sum, which ε a few ulps from 1 needs. Complex::ln and Complex::pow, a real power whose phase's multiple of π/4 is reduced exactly; the operators by name as saturating_add and its kin; mul_by_tau and sin_cos_range_turns; root_in, bracketed by the ITP method within a fixed budget; and kepler_universal, the universal anomaly by Newton's method kept to a geometric bracket by rtsafe's rule. The harness gains references for each: Gauss–Legendre panels and sums by parts for the kernels, a 960-bit Kepler solve, and powers with exact phases on the axes, plus sweep_some to sweep named functions alone. Every new function is within half an ulp on all 104 layouts. Comment density stays below main's: 5.78% of lines against 5.79%, and 10.76% of characters against 11.48%. A z one ulp off a diagonal with a huge modulus, as MIN + i·MAX on I104F24 squared, has a real part 2^−127 of |z|², which η measured from the real axis held only to its absolute precision. arg z is now j·π/4 + ψ about the nearest of the eight lines, ψ from (|y| − |x|)/(|x| + |y|) near a diagonal, so |ψ| ≤ π/8 keeps its relative precision; Complex::pow is within half an ulp on all 104 layouts.
…ernels The working arithmetic adds and multiplies two and three limbs in 256 bits, correctly rounded and so bit for bit as before, two-limb sums on a 192-bit path; division and square roots start from table seeds refined in Q62, and two-limb division from the exact 128/64 quotient. Each type also has a working type as wide as itself: a function keeps its result there when every output leaves 24 bits past the type's ulp, plus what it may cancel inside (48 for the Faddeeva function and complex erf, 24 for E₁ and Γ(a, z), 16 for the weighted kernels), else works again in the wider type, and the Jacobi functions check the amplitude too. Carlson's duplication stops once its closing series is exact, R_J and R_D closing with DLMF 19.36.2 through degree 7; R_C takes its series below 2^−4. The Jacobi functions take Bulirsch's descent, sharing one AGM with K, and the arctangent reduces (16y − j·x)/(16x + j·y) with one division. On 128-bit types with at most 96 fractional bits, sin_cos, exp, ln and atan2 run in Q126 integers after their exact reductions, within 2^−120 before the one rounding. On I64F64, against the start of this work: carlson_rj_roots 48.4 → 6.7 µs, jacobi_sn_cn_dn 15.4 → 2.5 µs, sin_cos 255 → 69 ns, exp 163 → 66 ns, ln 259 → 89 ns and atan2 370 → 116 ns, against 3.1.1's 28, 70, 140 and 99; on I32F32 the elementary functions run 2–5× faster and the elliptic ones 7–12×. Outputs move only for exprel, exp_moments, log1p_ratio, stumpff and Complex::exprel, all new in 4.0.0; every function stays within its bound on all 104 layouts, within half an ulp plus at most 2^−30 on 32-bit fractions and 2^−31 from 80 fractional bits up.
…st the edges CI's no-panic step counted lines starting `pub fn` in src/ops against the no_panic attributes there, so the indented methods of Complex, already verified by the link check, made the counts differ, 83 against 98, and the step fail. It now requires every non-const pub fn, methods included, to sit under its attribute, and Complex's conj, norm and arg and Modulus's four methods carry one, so the link check covers them as well. Tests reach the wraps whose exact value rounds onto an end of the range, the solvers' and elliptic functions' edge cases, complex zeros, and the coefficient tables built again at run time; the working sum's generic path past three limbs, which no type reaches, goes. No output changes.
Each kernel trajectory's closed forms are built from now has a range over an argument interval, rounded outward, so its solvers can be re-evaluated in intervals that provably hold their rounded outputs. exprel, exprel2, exp_moments, log1p and its ratios, tanh_ratio, atan_ratio, atanh_ratio, sqrt and hypot are monotone in x or |x| and take their ends; tan_ratio's range is everything when cos changes sign or the ends are π apart; atan2's is its corners' angles unless the box holds the origin or crosses the negative x-axis. The Stumpff functions fall up to z = 4, and c₃ everywhere; past it |cₖ''| ≤ 2/(k + 4)!, from their integral forms, rules an extreme out from the slopes at the ends or bounds it by that times w²/8 and by the extremes' own sizes, so the range over [20.18, 20.2] holds c₁'s first minimum within 1.7·10⁻⁶. Rounding::Down and Rounding::Up direct mul_rounded, div_rounded, mul_div_rounded, mul_add_rounded and dot_rounded from the exact products and sums their nearest forms use; mul_div and the rest are now their Rounding::Nearest, and no output moves. ellip_pi_increment(n, u, step, m) is ∫ dv/(1 − n·sn²v) from u to u + step for n < 1: the ends reduced by 2K with 2Π(n|m) a period, 1 − n·sn² as cn² + (1 − n)·sn² so it does not cancel, and a step far inside the nearest zero of 1 − n·sn² by five-point Gauss–Legendre, where a difference of the ends' integrals would cancel. It matches mpmath to 1e-15 on I64F64, a 2^−30 step to the last bit. The harness sweeps the directed forms against MPFR, failing any rounded the wrong way or an ulp off, and the increment against 512-bit references. Unit tests sample every range densely, crowded round the Stumpff extremes and the poles of tan. Comment density stays below main's: 5.76% of lines against 5.79%, and 11.01% of characters against 11.48%.
The full sweep of all 104 layouts adds the directed forms and ellip_pi_increment to worst.txt and the README. The increment is within half an ulp everywhere. The directed forms are below an ulp and on their own side at every point, checked exactly, though their maxima read 1.00 once rounded, which the README now says. No other function's worst case moves, and the mean errors against MPFR match the baseline's, so baseline.json stands.
|
Codecov Report❌ Patch coverage is
📢 Thoughts on this report? Let us know! |
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
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.
Releases 4.0.0, the first release since 3.1.1. This version is breaking:
CordicNumberis now sealed.What changes
ERROR_ULPS) on every type from I8F0 to I3F125, and within half an ulp wherever the sweep looked. Odd and even functions are exact, and monotone ones stay monotone. The 23 existing functions other thansqrtchange output, all as corrections; the changelog names each one.Roundingdown or up;ellip_pi_increment;Complexwith the Faddeeva function, erf, E₁ and Γ(a, z);root_inandkepler_universal;lnandatan2faster;sqrtis unchanged.src/opsto carry the annotation, methods included. The old line count failed on indentedComplexmethods.Commits: one per step, each passing its own tests, ending with the accuracy baseline update.
Checks, run locally:
check.shpasses, and comment density stays below main's.