Keep the argument precision in expint and gamma(a, x) - #548
Conversation
`expint(ν, z)` widened `Float16`, `Float32` and `ComplexF32` arguments to
`Float64` on several code paths, because `Float64` values leaked into the
computation:
* the near-pole correction in `En_expand_origin_general` evaluates the
`Float64` constants `π^2` and `π^4`,
* `gamma(1-ν)` there and `loggamma(n+1)` in `En_safe_expfact` are computed in
`Float64` for an integer order,
* so are `cospi`, `sinpi` and `logabsgamma` in `En_imagbranchcut`.
The results are converted back to the type that promoting the arguments
gives, which keeps the `Float64` precision of the intermediate computations
and only corrects the type of the result. `expint(2.000001f0, 0.5f0)` for
instance returned a `Float64` and now returns a `Float32`.
In `En_expand_origin_general` the guard value is additionally bound to a
variable so that the type of `n` in that branch, which is passed on to
`polygamma`, can be inferred at all.
`gamma(a, x)` and `loggamma(a, x)` are computed via `expint` for
non-positive-integer `a` and inherit this: `gamma(a::Float32, x::Float32)` is
now inferred as `Float32` instead of `Union{Float32,Float64}`.
All `Float64` results are unchanged, bit for bit.
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## master #548 +/- ##
==========================================
+ Coverage 94.49% 94.71% +0.22%
==========================================
Files 14 14
Lines 3016 3030 +14
==========================================
+ Hits 2850 2870 +20
+ Misses 166 160 -6
Flags with carried forward coverage won't be shown. Click here to find out more. ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
| gammaterm = gamma(1-ν)*z^(ν-1) | ||
| # `gamma(1-ν)` is computed in `Float64` for e.g. integer `ν`, so it is converted | ||
| # to the type of `z` (which `_expint` promoted with `ν`) to keep the result type | ||
| gammaterm = oftype(z, gamma(1-ν))*z^(ν-1) |
There was a problem hiding this comment.
Based on the function signature it's not clear that this is possible at all - it allows e.g. for z::Int.
There was a problem hiding this comment.
You are right, and it is worse than unclear — it actually breaks. On the version you reviewed:
julia> SpecialFunctions.En_expand_origin_general(0.5, 1, 1000)
ERROR: InexactError: Int64(1.772453850905516)where master happily returned a Float64. The public entry points are unaffected, because _expint does z, = promote(float(z), ν) first, but the signature admitted more than the body can handle.
Fixed in fbf8ca7 by putting that contract into the signatures of the four internal helpers that rely on it (En_expand_origin, En_expand_origin_general, En_expand_origin_posint, En_safe_expfact, En_imagbranchcut): z::Union{AbstractFloat,Complex{<:AbstractFloat}}. The call above is now a MethodError instead of an InexactError, and expint(0.5, 1) still works. ComplexOrReal{<:AbstractFloat} would have read better, but it is defined in gamma.jl, which is included after expint.jl.
(reply written by Claude Code on behalf of @andreasnoack)
| # `reνz` is bound to a variable such that the type of the values in this branch | ||
| # (in particular of `n`, which is passed on to `polygamma`) can be inferred |
There was a problem hiding this comment.
Even without this variable it should have been possible to infer the type of n.
There was a problem hiding this comment.
Correct, and my comment was simply wrong — real(round(ν)) - 1 is inferrable on its own:
julia> Base.infer_return_type(ν -> real(round(ν)) - 1, Tuple{BigInt})
BigIntWhat the binding actually buys is the narrowing of real(ν + z) to Union{Float64,Float32} inside the branch: the isa refinement applies to the bound value, so a recomputed real(ν + z) would still be Any. That is what allows n to be converted to the working precision, which keeps the correction — and hence the polygamma calls — out of BigInt/BigFloat, and incidentally is what removed the _polygamma(::Int64, ::BigFloat) finding from report_package.
Comment rewritten accordingly in fbf8ca7.
(reply written by Claude Code on behalf of @andreasnoack)
| # `loggamma(n+1)` is computed in `Float64` for integer `n`, so the result | ||
| # has to be converted back to the type of `z` |
There was a problem hiding this comment.
If z is a BigFloat, then also loggamma(n + 1) should be computed as BigFloat in the same precision, otherwise the results will be silently inaccurate. But maybe such cases can't occur?
There was a problem hiding this comment.
Such cases do occur, and you are right that the conversion was papering over a real loss. At the default 256 bits:
julia> setprecision(BigFloat, 256) do
v = SpecialFunctions.En_safe_expfact(199, BigFloat(1))
ref = -1/gamma(BigFloat(200))
-log2(abs(v - ref)/abs(ref))
end
43.4 # of 256 bitsloggamma(200) in Float64 is off by 8.6e-14 in absolute terms, and it goes straight into an exponent. The branch is reached from e.g. expint(200, big(1.0)) through En_expand_origin_posint.
That particular expint call is saved only by the term being about 1e-373 next to a sumterm of 5e-3, so the end result measured 255 of 256 correct bits against an independent reference — but the helper itself was wrong, and for abs(z) < 3 (this branch) the term is always negligible, so I could not construct a case where it shows up in expint.
fbf8ca7 evaluates the exponent in promote_type(Float64, real(typeof(z))): 43.4 → 248.7 bits, unchanged for the hardware floats.
(reply written by Claude Code on behalf of @andreasnoack)
| # `cospi`/`sinpi`/`logabsgamma` are computed in `Float64` for e.g. integer `ν`, | ||
| # so the result has to be converted back to the type of `z` |
There was a problem hiding this comment.
This one is observable, so thanks for pushing on it. At 256 bits, against a reference built from MPFR's expinti via E₁(-x + 0im) = -Ei(x) - πi and E_{n+1}(z) = (exp(-z) - z·E_n(z))/n — independent of the code under test:
| before | after fbf8ca7 | |
|---|---|---|
expint(2, big(-4)+0im) |
247.4 | 247.4 |
expint(5, big(-4)+0im) |
52.8 | 244.4 |
expint(9, big(-4)+0im) |
50.3 | 245.1 |
(of 256 bits). ν = 2 escaped only because logabsgamma(2) == 0 exactly; for ν = 5 the Float64 value 3.178... enters the exponent and caps the result at Float64 accuracy. cospi/sinpi are exact for an integer ν but have the same problem for a lower-precision float order, so they are converted too.
What remains on this path is the Float32 inaccuracy in #545, which is a different and much larger problem (up to 85% relative error), and the last ~12 bits at 256, which are the fixed 101-term limit in En_taylor from #546.
(reply written by Claude Code on behalf of @andreasnoack)
Follow-up to the review of #548. `En_safe_expfact` and `En_imagbranchcut` converted a `Float64` result back to the type of `z`, which is right when that type is narrower but hides a loss of accuracy when it is wider: `loggamma(n+1)` and `logabsgamma(ν)` enter an exponent, where their absolute error becomes a relative error of the result. Both now evaluate those terms in `promote_type(Float64, real(typeof(z)))`, i.e. unchanged for the hardware floats and in its own precision for a `BigFloat`. At the default 256 bits, `En_safe_expfact(199, big(1.0))` goes from 43 to 249 correct bits and `expint(5, big(-4)+0im)` from 53 to 244. The internal functions that rely on `z` having been promoted with `ν` by `_expint` now say so in their signatures instead of assuming it, so that `oftype(z, ...)` cannot be reached with something like `z::Int`. Also corrects the comment on the bound guard value: the type of `n` was always inferrable; what the binding buys is narrowing `real(ν + z)` to `Union{Float64,Float32}` inside the branch, which is what keeps the correction and its `polygamma` calls in the working precision. All values of the previous version are reproduced bit for bit.
|
@devmotion don't spend more time on this for now. Let me take a closer look at what Claude proposed before continuing with the review. |
|
Closing in favour of a different approach. @devmotion's review here was the trigger: the conversions in this PR let a kernel receive arguments of two different precisions and then converted the result, which hides the mixing rather than preventing it — and, as his third and fourth comments predicted, it papered over a real loss of accuracy for The replacement restricts the internal helpers to arguments of a single precision, #549 (the type stability sweep) no longer depends on this branch — it is now based on Comment written by Claude Code on behalf of @andreasnoack. |
2 participants
expint(ν, z)widensFloat16,Float32andComplexF32arguments toFloat64on several code paths, becauseFloat64values leak into the computation. On master:and correspondingly
Base.infer_return_type(gamma, Tuple{Float32,Float32})isUnion{Float32, Float64}rather thanFloat32, which is what #520 ran into.There are four separate leaks:
En_expand_origin_generalevaluates theFloat64constantsπ^2andπ^4,gamma(1-ν)in the same function is computed inFloat64for an integer orRationalorder,loggamma(n+1)in then >= 100branch ofEn_safe_expfactlikewise,cospi,sinpiandlogabsgammainEn_imagbranchcut, which handles a complex argument sitting exactly on the negative real axis.Each is fixed by converting the result back to the type that promoting the arguments gives, which keeps the
Float64precision of the intermediate computations and only corrects the type. InEn_expand_origin_generalthe guard value is additionally bound to a variable, so that the type ofnin that branch — which is passed on topolygamma— can be inferred at all.gamma(a, x)andloggamma(a, x)are computed viaexpintfor non-positive-integeraand inherit all of this.Effect
Every combination of order type × argument type I checked now infers concretely:
Int,Int32,Rational,Float16/32/64,ComplexF32/64andBigFloatorders againstFloat16/32/64,ComplexF32/64,BigFloatandComplex{BigFloat}arguments, forexpint,expintxandgamma(a, x).Float64results are unchanged bit for bit — I compared 144 values acrossexpint,expintx,expinti,gamma(a, x),zeta,digamma,polygamma,eta,invdigammaandloggamma. The only values that change are the ones that had the wrong type, and they agree with the old ones toFloat32rounding.Tests
test/expint.jlcovering the series about the origin, the near-pole correction, the continued fraction and the negative-real-axis procedure forFloat16,Float32andFloat64, real and complex;gamma(a, x)/loggamma(a, x)type stability and allocation tests in the existingGPU compatibilitytestset oftest/gamma_inc.jl, in the shape #520 asked for.Full suite passes on 1.10, 1.11, 1.12 and 1.13-rc3.
While verifying this I found that
expintis separately inaccurate forFloat32arguments near the negative real axis (up to 85% relative error) — that is pre-existing, unrelated to the types, and filed as #545.Prepared by Claude Code on behalf of @andreasnoack; the investigation and the text above are Claude's.