Accuracy and domains¤

The supported input domain and measured accuracy of every public function. This page is the canonical home for these facts; other pages link here rather than restating them.

Every function is tested against a reference implementation — scipy.special where a counterpart exists, mpmath where it does not. All entries assume double precision (how to enable it).

Functions¤

Per-function limits¤

eval_gegenbauer, eval_gegenbauers : At exactly \(\alpha = 0\), spexial returns \(C_0^{(0)} = 1\) and \(C_n^{(0)} = 0\) for \(n \ge 1\), which is what the generating function gives. SciPy agrees up to 1.14; from 1.18 it returns 0.0 for every order at exactly \(\alpha = 0\), while still returning 1.0 at \(\alpha = 10^{-300}\).

eval_gegenbauer : Error is proportional to the recurrence's largest intermediate value, not to the size of the result. At \(n = 20,\ \alpha = 10\) the intermediates peak at \(9.6\times10^{7}\) on the way to a result of \(1.3\times10^{-2}\), and the worst absolute error near a root there is \(1.5\times10^{-7}\) — \(1.6\times10^{-15}\) of that peak. Expect absolute error of roughly \(10^{-15}\) times the peak intermediate, which at large \(\alpha\) can be far larger than \(10^{-15}\) times the answer. (Why.)

eval_gegenbauer : At \(x = \pm\infty\) (outside the supported \(\lvert x \rvert \le 1\)) every order returns its analytic limit. The leading coefficient is \(2^n(\alpha)_n/n!\), and for \(\alpha > -1/2\) every factor after the first is positive, so its sign is \(\operatorname{sign}(\alpha)\) — giving \(C_n(+\infty) = \operatorname{sign}(\alpha)\infty\) and \(C_n(-\infty) = \operatorname{sign}(\alpha)(-1)^n\infty\), and \(0\) for \(n \ge 1\) at \(\alpha = 0\), where the polynomial vanishes identically. SciPy differs here: it returns inf at \(+\infty\) but nan at \(-\infty\). A finite argument large enough to overflow the polynomial gives nan. The gradient at \(x = \pm\infty\) is 0 in both modes. It was nan in reverse and 0 in forward: the where that substitutes the limit transposed a zero cotangent into the \(2\alpha x\) branch, forming \(0 \cdot \infty\). Forward mode had it right — at fixed \(x = \pm\infty\) the limit is \(\operatorname{sign}(\alpha)\infty\), a step in \(\alpha\), whose derivative is \(0\) away from the jump — so the recurrence is now kept away from an infinite \(x\) altogether and the two agree. A nan \(\alpha\) still splits the modes, and is left alone as the degenerate input it is.

incomplete_beta : This is the unregularized \(B(a, b, z) = \int_0^z t^{a-1}(1-t)^{b-1}\,\mathrm{d}t\) (DLMF 8.17.1), not SciPy's betainc, which is the regularized \(I_z(a, b)\). The two differ by the complete beta function \(B(a, b)\). Reconstructing this as beta(a, b) * betainc(a, b, z) is nan for every \(b \le 0\) — \(B(a, b)\) has a pole there while the product does not — and that domain is the reason the function exists. (Why not hyp2f1.)

incomplete_beta : Two fixed-length 64-term series, switched at \(z = 1/2\), each converging like \(2^{-k}\). They meet to about \(10^{-11}\) relative rather than exactly, since they truncate differently; worst measured step across the seam is \(2.9\times10^{-11}\), at \(a = 4.5,\ b = -2.5\). That is inside the tolerance the reference tests assert, so the seam is not the accuracy-limiting feature. Above \(a \approx 16\) the error stops improving with term count at all — it is cancellation between large alternating terms, not truncation.

incomplete_beta : a and b must be scalars; z may be any shape. The \(z\)-derivative comes from a jax.custom_jvp and is exact: by Leibniz it is the integrand at the endpoint, \(z^{a-1}(1-z)^{b-1}\). The a and b tangents have no cheap closed form and fall back to differentiating the series, but only when actually requested — symbolic_zeros skips them in the common case of differentiating with respect to \(z\) alone. It is a custom_jvp and not a custom_vjp specifically so that jax.hessian's jacfwd(jacrev(...)) still composes.

sph_legendre_p, sph_harm_y : The degree and order are static Python ints and may not be arrays, unlike SciPy's, which broadcast them. The returned shape is therefore \(\theta\)'s own rather than a broadcast against \(n\) and \(m\). This is not only a simplification: it is what makes the values right. jax.scipy.special.sph_harm_y takes them as arrays and pairs n[i] with theta[i] positionally, so a length-1 degree against a batch of angles is correct at index 0 and wrong everywhere else, by up to 1.18 absolute for \(n \le 3\). (Why.)

sph_legendre_p, sph_harm_y : Derivatives are finite at both poles, where jax.scipy.special.sph_harm_y's are nan for every \(n \ge 1\). Upstream differentiates \(\sqrt{1 - \cos^2\theta}\), whose slope is infinite at \(\theta = 0, \pi\); the surviving \(0 \cdot \infty\) is nan, even though \(\partial_\theta Y_n^m\) is perfectly finite there. Here \(\sin^m\theta\) is built as an integer power of \(\sin\theta\), so the square root is never formed and lax.integer_pow supplies an exact derivative at zero. The \(n = m = 0\) case, where the function is constant and the true derivative is zero, was nan upstream too and has since been fixed; \(n \ge 1\) has not.

sph_legendre_p : Outside SciPy's documented \(\theta \in [0, \pi]\) the two differ by a sign for odd \(m\). SciPy carries \((1-u^2)^{m/2}\), which is \(\lvert\sin\theta\rvert^m\) and therefore even about \(\theta = 0\); this carries \(\sin^m\theta\), the analytic continuation, which is odd for odd \(m\). On \([0, \pi]\) they agree to the tolerance above. The continuation is the reason the one-sided pole derivative is a real number rather than a cusp — a central difference through \(\theta = 0\) against SciPy reports zero for every odd \(m\) regardless of the truth.

sph_legendre_p : Finite and accurate past \(m = 200\), where the unnormalized \(p_m^m = (-1)^m(2m-1)!!\) overflows float64 near \(m = 90\). The normalization is folded into the recurrence seed, computed in log space, so that quantity is never materialized; built the obvious way — recurrence first, \(N_{lm}\) afterwards — everything from \(m = 90\) up is $\infty \cdot 0 = $ nan, and moderate \(m\) loses about two digits to cancellation. Worst measured error on the diagonal is \(1.5\times10^{-13}\) absolute at \(n = m = 200\).

sph_harm_y_cart, sph_harm_y_cart_all : The direction is assumed already normalized and is not normalized internally. That is deliberate twice over: it avoids repeating a caller's own work, and it leaves the \(r = 0\) policy with the caller. A zero vector is well defined here — \(Y_0^0 = N_{00}\), and zero for every \(m \ge 1\) — and remains twice differentiable, which normalizing internally would turn into nan. A non-unit direction is not rejected; it simply evaluates the formula, which is the harmonic scaled by a power of the norm.

sph_harm_y_cart, sph_harm_y_cart_all : These exist because the Cartesian gradient on the z-axis is exactly 0.0 for every \(m \ge 1\) term of any harmonic evaluated through \(\theta\) and \(\phi\) — neither angle has a directional derivative there — against a non-zero true limit. That is a property of the coordinates, not of an implementation, so SciPy and JAX share it and no upstream fix would remove it. Verified against a finite difference taken along \(x\) from the pole, agreeing to \(10^{-5}\) with the step size used.

sph_harm_y_cart_all_terms : Same values, same indexing and the same layout as sph_harm_y_cart_all, returned as a nested tuple of arrays rather than one stacked array. Reach for it whenever you are going to reduce over the table — summing \(\sum_{lm} c_{lm} Y_l^m\), say. Indexing a stacked table stops XLA folding each term into the reduction as it is produced, so the whole table is materialized: measured on a multipole expansion at \(n = 12\) over a million directions, 17.7 s through the stacked form against 10 ms through these terms, for identical values. It is also the one public function here not wrapped in jax.jit, deliberately — a jitted function returning a pytree materializes each leaf at the call boundary, which is the fusion it exists to preserve. jit around it instead.

sph_harm_y_cart_all : Shape and index layout follow scipy.special.sph_harm_y_all: the result is \((n+1, 2m+1, \ldots)\) and negative orders live at the end of the second axis, so Y[l, -k] reaches them by ordinary negative indexing. Entries with \(\lvert j \rvert > i\) are zero, since no such harmonic exists.

eval_gegenbauers : alpha and x broadcast against each other, and the orders are stacked on a new leading axis — so scalars give \((n+1,)\) and arrays give \((n+1, \ldots)\). It was scalar-only before, guarded by a ValueError, because the orders were glued together with hstack, which concatenates rather than stacks once the operands have an axis of their own.

comb : Returns 0 for k > N, k < 0 and N < 0, matching scipy.special.comb. comb(inf, 0) is 1, comb(inf, k >= 1) is inf and comb(inf, inf) is nan, all as in SciPy. A subnormal k gets the k = 0 answer of 1, where the limit is inf — at every width except float16, whose subnormals are ordinary normal float32 numbers by the time the computation runs, so float16 alone returns the mathematically correct inf. The rest: XLA flushes it, and no bit test can tell a flushed subnormal from a zero — one that reads the bits is right eagerly and collapses under jit, which is worse than a uniform floor. (is_negative survives the same treatment only because a flushed negative keeps its sign bit as -0.0.) The inexact variant only; there is no exact=True path. Why not jax.scipy.special.comb? It computes the same log-gamma difference and inherits its collapse: measured against mpmath, it is \(3.7\times10^{-6}\) out at \(N = 10^{10}\), wrong by a factor of \(10^{23}\) at \(10^{16}\), returns 1.0 at \(10^{18}\) and nan at DBL_MAX. It is also 70x out in bfloat16 and admits a negative subnormal k. Below \(N = 1000\) the two agree bit for bit, since that is the same formula. Two formulas are stitched at \(N = 1000\) in float64, and at \(N = 10\) in float32 and narrower — the log-gamma cancellation grows in units of eps, so it bites \(10^9\) times sooner there, while the Beta form's own error is dtype-independent. float16 and bfloat16 are computed in float32 and rounded back. The two are: a log-gamma difference below the join, which is the more accurate there but cancels catastrophically above, and a Beta-function form above, which does not cancel. The step in value across the join is a few times \(10^{-13}\). From about \(N = 1/\varepsilon\) (\(4.5\times10^{15}\) in float64, \(8.4\times10^{6}\) in float32) the Beta form is evaluated from its asymptotic form instead — the switch is \(\mathrm{small}(\mathrm{small}+1) < 2\,\mathrm{big}\,\varepsilon\), i.e. wherever the dropped correction and its slope are both below an eps. The \(+1\) is load-bearing: the correction is \(\mathrm{small}(\mathrm{small}-1)/(2\,\mathrm{big})\), which vanishes identically at \(k = 0\) and \(k = N\) (where \(\mathrm{small}\) is 1) while its derivative, \(1/(2(N+1))\), does not — so testing the value alone fired this branch at both ends of the \(k\) range for every \(N\) past the join and left jax.grad(comb) low by exactly that, \(6.7\times10^{-5}\) relative at \(N = 1001\) and \(1.4\times10^{-2}\) in float32 at \(N = 11\). Widening the branch is what makes the derivative right at large \(N\), because jax.scipy.special.betaln divides its smaller argument by its larger one and XLA flushes that quotient to zero once it is subnormal — which silently cost a factor of \(e^{-2}\), 86%, at comb(N, 1). In float32, comb(FLT_MAX, 1) is inf rather than FLT_MAX: the value is exp of a logarithm near 88.7, which rounds up out of range. comb holds to \(2.2\times10^{-13}\) out to DBL_MAX, worst near \(N \approx 7\times10^{151}\). That is a floor, not a tuning choice: the value is the exp of a logarithm near 699, and that round trip alone costs \(\sim 8\times10^{-14}\) whatever the formula does.

comb : comb(N, 0) and comb(N, N) are 1 to \(6\times10^{-15}\) rather than exactly 1 for \(1000 < N \lesssim 3\times10^{6}\) in float64 — 28 ulps, well inside the stated rtol. This is a deliberate trade. The asymptotic branch returns \(\exp(\log(N+1) - \log(N+1))\), which is exactly 1; routing those two points through betaln instead is what makes the derivative right there, and the derivative was previously wrong by \(6.7\times10^{-5}\). Exactness in the value at two points is worth less than a correct slope along the whole line. In float32 the two branches agree and nothing changed — but see the gammaln note below, which costs float32 far more than this.

comb : Without x64, comb(N, 0) and comb(N, N) fall up to 8 ulps short of the exact 1 they should be, and comb(5, 2) is 9.999995. The cause is jax.scipy.special.gammaln returning \(2^{-21}\) instead of 0 at a float32 argument of exactly 1 — an upstream defect at an exact zero of the function, filed as #29, and \(\exp(-\log\Gamma(1))\) is exactly what comb(N, 0) reduces to. The shortfall is not one number: measured across \(N\) from 0 to FLT_MAX it is 0, 4, 6 or 8 ulps depending on \(N\), as the two branch cross-overs and gammaln's own rounding move beneath it, and it is 0 from about \(N = 8.9\times10^{6}\) where the asymptotic form takes over. float16 and bfloat16 round it away at every \(N\); float64 is exact throughout.

comb : The second derivative with respect to \(N\) has the wrong sign from \(N \approx 1.3\times10^{154}\) in float64 and \(N \approx 1.8\times10^{19}\) in float32. Both surviving terms are of size \(1/N^2\), which is subnormal at those magnitudes, and XLA flushes subnormals: jax.grad(jax.grad(log))(1e154) is -0.0 where the true \(-10^{-308}\) is a perfectly representable denormal. The dominant negative term vanishes and the positive one is left standing, so the magnitude is right and the sign is inverted. No regrouping recovers it — the quantity itself is below the platform's floor — and only a rule that never goes through log space could, which comb cannot have: the closed form \(\binom{N}{k}(\psi(N-k+1) - \psi(k+1))\) cancels catastrophically, measuring \(1.8\times10^{-5}\) at \(N = 10^{10}\) against differentiating the implementation's \(10^{-15}\). The first derivative is unaffected everywhere. Where the value has already overflowed to inf, the first derivative is nan in reverse mode and inf in forward mode; there is no right answer past the dtype's range, and the two modes are simply not obliged to agree about which infinity they lost.

eval_gegenbauer : Batched evaluation is not bit-identical to a scalar call — XLA re-associates the lax.scan recurrence once there is a leading axis, so an array or vmap argument can differ by an ulp in both value and gradient, and the absolute figure above roughly doubles. A subnormal alpha is read as zero, so the limit at \(x = \pm\infty\) comes back 0 rather than \(\pm\infty\) — uniformly at every order, and identically in both entry points. A bit-level sign test recovers it eagerly and not otherwise: both entry points are unconditionally jax.jit, and inside that fused kernel the operand has already been flushed, so the test reached the same 0 by a longer route while claiming to have fixed the case. It has been withdrawn; this is the same platform floor the subnormal x below sits on.

spence : The second derivative is inf across the subnormal band, where it is finite and representable for most of it (\(4.5\times10^{307}\) at \(z = 2.2\times10^{-308}\)). Separating the pole from the subnormals needs a bit-level test, and one does not survive XLA's fusion — dtype.exactly_zero is correct in isolation and collapses when a select is its only consumer, which is how the first derivative came to be \(-\infty\) across the same band under jit while eager was right. The first derivative is now exact there, in both modes and for complex input as well as real — log_no_flush takes a complex argument apart by components, which keep their bits, rather than bitcasting the pair. It reconstructs whenever either component is subnormal, not only when both are: with a real part of exactly tiny beside a subnormal imaginary part the two are within a factor of two, and requiring both cost the logarithm its entire imaginary part, returning an argument of 0 where \(\pi/4\) was right. The second is a documented ceiling.

eval_gegenbauer : A subnormal x is flushed by XLA in the multiply itself, so 2*alpha*x loses it, and a subnormal alpha is read as zero, so the limit at \(x = \pm\infty\) comes back 0 at every order. It is recoverable in principle: regrouping the product around the mantissa bits gets eval_gegenbauer(1, 1e300, 5e-324) exactly right. It is not worth what it costs. The regrouping needs a lax.optimization_barrier to survive XLA's reassociation, and that barrier blocks the fusion the recurrence depends on — 1.58x slower over 128 points — while the where over bitcast integers returns a strongly typed result, so ordinary calls stop being weakly typed and promote differently downstream. Both costs fall on every caller; the loss falls only where \(\alpha\) is around \(10^{285}\), since below that the lost quantity is itself of denormal size. So the result is the value at x = 0. Same XLA floor as gamma below \(x \approx -170.6\). Note SciPy returns nan for several of these where spexial returns the correct limit.

gamma : Delegates the value to jax.scipy.special.gamma, so it cannot drift from upstream; spexial supplies only the derivative. Accepts real and complex input, the latter from jax 0.10.2 (below that JAX branches on floor(x) and raises). Returns inf at x = 0 and nan at the negative integers, matching JAX and scipy.special.gamma from 1.18 — SciPy is not a stable reference at the poles, returning inf everywhere up to 1.14. There is no near-pole blow-up: measured error stays at \(10^{-16}\)–\(8\times10^{-14}\) right up to \(10^{-8}\) from a pole. Below \(x \approx -170.6\) the true value is subnormal and XLA on CPU flushes it to zero, so spexial returns 0 where SciPy returns a denormal. In the other direction a subnormal argument is supported, which upstream does not manage: \(\Gamma(x) \to 1/x\) there, still representable for the factor of about two between tiny and \(1/\mathrm{max}\) — in float32 the reachable band \(2.9\times10^{-39}\) to \(1.2\times10^{-38}\) — where jax.scipy.special.gamma returns inf. The lower edge is exactly \(1/\mathrm{max}\), \(2.9387\times10^{-39}\) in float32. The substitution is made across the whole subnormal band, at every width. An earlier revision gated it on upstream returning a non-finite value, on the strength of a float16 comparison that had been read backwards: over the 767 float16 subnormals where jax.scipy.special.gamma is finite, this branch is the closer of the two to mpmath at 690 of them and worse at 5, worst case \(4.9\times10^{-4}\) against upstream's \(4.2\times10^{-3}\). The logarithm behind it is evaluated at the default float width rather than the argument's, because a bfloat16 subnormal is a float32 subnormal too — the two share an exponent range — and \(\log\lvert x\rvert \approx -87\) has no room in a dtype whose spacing there is \(0.5\). The first derivative is accurate to \(1.6\times10^{-13}\) everywhere tested, but the second and higher are not usable on the negative axis: \(\Gamma''\) routes through jax.scipy.special.digamma's own derivative, and JAX's trigamma diverges from the truth from about \(x = -7.5\) (at \(x = -10.5\) it is wrong by \(10^{11}\) relative, and at \(-20.5\) it has the wrong sign). This is upstream — jax.grad(jax.grad(jax.scipy.special.gamma)) returns the identical wrong number — but spexial inherits it. gamma(+inf) is inf, matching SciPy; gamma(-inf) is nan — Gamma has a pole at every negative integer, so the limit does not exist, and SciPy's -inf is not something to copy.

k0, k1, k2, k0e, k1e, k2e : Worst relative error is \(2.0 \times 10^{-7}\) for k0, \(1.8 \times 10^{-7}\) for k1 and \(1.3 \times 10^{-7}\) for k2, at \(z = 8.998353\). The profile there is spiky, not smooth — half a thousandth away, at \(z = 8.9932\), the error is only \(7.9\times10^{-8}\) — so a coarse grid understates it. Away from that point it falls in stages: ~\(8\times10^{-9}\) out to \(z = 15\), ~\(2\times10^{-13}\) to \(z = 30\), ~\(10^{-15}\) beyond. (Why the peak is where it is.)

Upper limit. k0, k1 and k2 underflow to 0 above \(z \approx 705.3\), where the true value is smaller than the smallest normal double. That number is float64's, and the ceiling is a property of the dtype, not of the function: it is \(z \approx 85.3\) in float32, \(85.2\) in bfloat16 (which shares float32's exponent range) and \(z \approx 16.2\) in float16. A float16 caller reading 705 would be four decades out. SciPy is not uniform here: scipy.special.kn underflows at \(z \approx 698\), while k0 and k1 return denormals out to \(z \approx 745\). Use k0e, k1e or k2e above \(z \approx 705\): \(e^z K_n(z)\) decays only as \(1/\sqrt z\) and stays accurate to ~\(10^{-15}\) up to DBL_MAX (\(1.798\times10^{308}\)), where scipy.special.kve returns nan.

In float32, the cross-over moves to \(z = 4.65\), the worst error rises to \(7.1\times10^{-3}\), and the unscaled ceiling drops to \(z \approx 85.3\) — the point where \(K_0\) falls below float32's smallest normal, \(1.18\times10^{-38}\). float16 and bfloat16 are computed in float32 and rounded back to the caller's dtype. The scaled forms have no ceiling at any width.

Derivatives of the unscaled family are exact to second order wherever the values are, the \(z = 699\)–\(705\) tail included, and grad(grad) at the pole is inf for all six. Three caveats:

• A subnormal argument is supported. XLA on CPU flushes one to zero, which made every k0 below tiny return inf; the logarithm is taken from the mantissa bits instead, so k0 matches SciPy across the band and at \(5\times10^{-324}\) returns 744.556 where SciPy overflows. k1 follows \(1/z\) there and is exact wherever \(1/z\) is representable — a band about a factor of two wide, since \(\mathrm{tiny} \times \mathrm{max} \approx 2\) in any IEEE format, which in float32 is \(2.9\times10^{-39}\) to \(1.2\times10^{-38}\). Past it, and for k2 throughout, the true value overflows and inf is the right answer.
• For k1 and k2, orders three and above lose accuracy as a ramp rather than a step: \(3\times10^{-9}\) from \(z \approx 686\) through 690, then a jump: \(3\times10^{-6}\) by 695 and \(7\times10^{-4}\) from 700 on. Orders 1 and 2 stay at \(2\times10^{-16}\), and k0 is unaffected at every order — it holds \(1.6\times10^{-16}\) through the same band.
• The scaled forms lose relative accuracy to cancellation, since \((e^zK_0)' = G_0 - G_1\) subtracts two nearly equal numbers: the first derivative passes \(10^{-6}\) relative at \(z \approx 3\times10^{10}\) and is meaningless by \(10^{16}\); the second costs about two decades near the cross-over (\(3.6\times10^{-5}\) at \(z = 8.9\)). Absolute error stays at machine precision, and the values remain exact to \(1.6\times10^{-16}\) at \(10^{300}\) — it is only the ratio that degrades.
• The third derivative of k0e/k1e/k2e is unreliable for any small \(z\), not merely at the pole: at \(z = 10^{-200}\) and \(10^{-100}\) one mode returns \(0\) or nan where the unscaled functions correctly give \(-\infty\), and a plausible finite \(0\) is worse than a disagreement. The two modes disagree in both directions — at \(z = 10^{-300}\) reverse mode gives k0e the correct \(-\infty\) while forward mode gives nan, and for k1e/k2e the reverse holds — so neither mode is the reliable one there. Second order at the pole is guaranteed; third is not, in either mode.

polylog : At the pole \(z = 1\) the value is inf for \(n = 1\) and finite above it, but the second derivative is nan from \(n = 2\) — the identity \(\mathrm{Li}_n' = \mathrm{Li}_{n-1}/z\) becomes \(\infty - \infty\) there, where the limit is \(+\infty\). The neighbourhood is correct (\(\mathrm{Li}_2''(1 - 10^{-6}) = 999987.2\)). Broadcasts over z (how). Raises ValueError for an order below 1, and for complex z — two of the three branches take the real part of a complex intermediate, which is exact for real z and would silently discard a genuine imaginary part. A non-integer order, including a whole-number float such as polylog(2.0, z), is a TypeError from the runtime type checker rather than a ValueError. For \(\lvert z \rvert \ge 2\) the order is capped at 60 by the Bernoulli table the inversion formula needs; past that the result is nan. Smaller \(\lvert z \rvert\) has no such ceiling: \(\Gamma(n+1)\) overflowing sends the harmonic term to zero, which is its correct limit, so polylog(500, 1.5) is exact. polylog(1, 1) is the pole and returns inf.

zeta : Bernoulli numbers are computed from exact fractions.Fraction arithmetic, not jax.scipy.special.bernoulli. At and above \(n = 54\) the result is the constant 1.0, which is not an approximation: \(\zeta(n) - 1 \approx 2^{-n}\) falls below half an eps of 1 once \(n > 53\), so every double-precision value from there up is 1.0. Taking the constant also steps around jax.scipy.special.zeta, which returns nan above \(n \approx 10^{15}\); spexial is correct at every magnitude including inf, matching SciPy. jax.grad(zeta) is meaningful except at two sets of points: the tabulated integers \(0 \ge n \ge -60\), where the value is a table lookup that carries no information about how \(\zeta\) varies between entries, and the negative even integers at any magnitude, which are a constant 0. Both report a finite number that is not \(\zeta'\). Everywhere else — non-integers, and the odd integers past the table — the functional equation is differentiated and the gradient is genuine to about \(10^{-14}\).

spence : Accepts real and complex argument; jax.scipy.special.spence is real-only and raises on complex. Its derivative, \(\log z/(1-z)\), is supplied analytically — evaluated as the limit \(-1\) at \(z = 1\), where the closed form is \(0/0\), and \(-\infty\) at \(z = 0\) (on the complex path the derivative itself is a clean \(-\infty + 0i\), but the JVP then multiplies it by the tangent, and a complex multiply forms \(\infty \cdot 0\) in its cross term: the imaginary part comes back nan eagerly and 0 under jit, since whether XLA folds the multiply by 1 is a property of the graph. Either way it is an infinity at a pole) — which matters beyond speed: JAX's own spence differentiates to nan across roughly \(1 < z < 2\), where spexial is exact. Do not compare against SciPy's complex spence near \(z = 3 \pm \sqrt3\). It returns 0.01125 at \(3 - \sqrt3\) where the true value is \(-0.25186\); spexial returns the true value. Use mpmath.polylog(2, 1 - z) as the reference at those two points. float16 and bfloat16 are computed in float32 and rounded back, as the Bessel functions are; jax.scipy.special.spence rejects both outright. At \(z = 0\) the first derivative is \(-\infty\) and the second \(+\infty\), both the true limits, but the third and higher orders return 0 there rather than diverging — \(z = 0\) is a genuine pole, where every order diverges and a substituted constant is the only value available, so fixing one order simply moves the problem to the next. That 0 covers the whole flushed band, not the single point: the guard behind it is the flushed comparison, which is True for every subnormal. Above the band the third derivative is \(-1/z^2\), correct to \(10^{-7}\) at \(z = 10^{-4}\) and returning \(-\infty\) once that has overflowed, from about \(z = 1.5\times10^{-154}\) down — one binade wider than the true overflow at \(7.5\times10^{-155}\), and -inf rather than the nan the fused \(1/(z(1-z))\) used to give. Neither quotient in the second-derivative rule is fused, and both halves matter: differentiating \(a/b\) forms \(a\,\mathrm{d}b/b^2\), so a squared denominator becomes a fourth power one order up. Splitting only the first of the two left \(\log z/(1-z)^2\) raising \((1-z)\) to the fourth, which overflows from \(z = 8.2\times10^{76}\) — dropping the dominant term and giving the third derivative the wrong sign in every mode at once, from \(z = 3.1\times10^{9}\) in float32. SciPy's real path is unaffected and agrees everywhere. There is no upper limit on the real argument. Above \(z = 1/\mathrm{tiny}\) (\(4.5\times10^{307}\) in float64, \(8.5\times10^{37}\) in float32) the value comes from the inversion formula \(-\pi^2/6 - \log(z)^2/2\) rather than from jax.scipy.special.spence, whose opening \(1/z\) is subnormal there and is flushed into the wrong branch, giving nan across the top two binades of every width. The dropped term is of order \(1/z\), i.e. \(10^{-313}\) of the value at the threshold, so the substitute agrees with SciPy to the last bit at \(1/\mathrm{tiny}\), \(10^{308}\) and DBL_MAX. spence(inf) stays nan, as in SciPy, and so does its derivative — \(\log z/(1-z)\) tends to \(0^-\) there, but a nan value with a 0 derivative would be the less coherent pair, and SciPy's nan at infinity is a convention rather than a limit.

zeta coverage against SciPy¤

Every real argument is covered, and SciPy agrees throughout — including where the answer is ±inf. \(\lvert\zeta(n)\rvert\) passes DBL_MAX at \(n \approx -260.2\), so from there down only the trivial zeros and their immediate neighbours are finite; that is the true value overflowing a double, not a limitation here. Accuracy is worst just off a negative even integer, where the \(\sin(\pi n/2)\) of the functional equation is near a zero of its own: \(1.8\times10^{-13}\) at \(n = -99.99\). SciPy is \(2\times10^{-4}\) there, so this is the better of the two.

Reference implementations¤

Outside the supported domain¤

Traced JAX code cannot raise. A function called outside its domain returns nan or inf, and that value propagates through jit, vmap and grad. Validate inputs before the call, or test the output with jnp.isnan. For why the library behaves this way, see About domain edges.