#4856 Bigfloat atan/asin/acos/atanh/asinh/sinh/tanh lose several ulp: the big-float-* evaluators use no guard digits

[David Scherfgen]

Spotted by Claude:

Summary

Applying an inverse-trigonometric or hyperbolic function to a bigfloat is several ulp less accurate than obtaining the same mathematical quantity through the ordinary bfloat path. The big-float-* functions that back these operations evaluate their fp… primitive at the working precision with no guard digits, whereas the bfloat(f(x)) path goes through the guarded *fp… evaluators (which add 8 guard digits). atan is the clearest demonstration.

Reproducer

block([fpprec:32],
  a: atan(bfloat(sqrt(2))),   /* big-float-atan  — UNGUARDED */
  b: bfloat(atan(sqrt(2))),   /* *fpatan         — GUARDED   */
  [a, b, a - b]);

a and b are the same value: atan(√2) at 32-digit precision, both evaluated from the identical input bfloat(sqrt(2)). Therefore a - b must be 0b0.

Expected: a - b = 0b0.

Actual: a - b is a nonzero bigfloat of magnitude ≈ 6b-33 — about four units in the last place at fpprec:32 (observed -6.16b-33 on an SBCL build; the exact trailing digits are build-dependent, but the several-ulp discrepancy is robust).

b is the correct value. atan(√2) is the "magic angle":

atan(sqrt(2)) = 0.95531661812450927816385710251575775424341...

Rounded to 32 significant digits that is 9.5531661812450927816385710251576b-1, which is exactly what b (the guarded path) returns. a (the unguarded big-float-atan) is ~4 ulp off.

Root cause

atan(<bigfloat>) is dispatched through the *big-float-op* table (src/trigi.lisp) to big-float-atan, whose real-argument branch is just:

(bcons (fpatan (cdr x)))

It calls fpatan at the current fpprec with no extra precision. fpatan(√2) internally forms %pi/2 - atan(1/√2) and sums a ~70-term series; with no guard digits the accumulated rounding is ~5 ulp, and the result is additionally hypersensitive (a 1-ulp change in the underlying adder shifts it by several ulp).

By contrast, bfloat(atan(sqrt(2))) keeps atan(sqrt(2)) symbolic and evaluates it via %atan's floatprog property → atanbigfloat → *fpatan, which does it correctly:

(defun *fpatan (a y)
  (fpend (let ((fpprec (+ 8. fpprec)))     ; 8 guard digits
           (if (null y)
               (if ($bfloatp a) (fpatan (cdr ($bfloat a)))
                   (list '(%atan) a))
               (fpatan2 ...)))))

— it raises fpprec by 8, evaluates, and rounds back once with fpend. big-float-log does the same (it binds fpprec (+ fpprec extra) with extra = 8). The big-float-* inverse-trig/hyperbolic evaluators never received that treatment.

Scope

The identical "no guard digits" pattern is present in every one of these src/float.lisp evaluators — each calls its fp… primitive directly at the working precision:

big-float-atan   -> fpatan
big-float-asin   -> fpasin
big-float-acos   -> fpacos
big-float-atanh  -> fpatanh
big-float-asinh  -> fpasinh
big-float-sinh   -> fpsinh
big-float-tanh   -> fptanh

so asin, acos, atanh, asinh, sinh, and tanh of bigfloats are affected as well. (big-float-sqrt → fproot has only 2 guard bits, which is why bfloat(sqrt(2)) itself is ~1 ulp high — arguably worth bumping too.)

Suggested fix

Give each real-argument branch guard digits the way *fpatan and big-float-log already do: raise fpprec by ~8, evaluate the primitive, and round the (real or complex) result back to fpprec. For atan the branch can simply reuse the existing guarded evaluator, (*fpatan x nil).

Implementer gotcha: fpatanh, fpasin-core, complex-atanh, and complex-sqrt read the constants *bfhalf* / *bfmhalf* directly (e.g. (fptimes* (cdr *bfhalf*) …)). These are precision-pinned globals set only by fpprec1, so naively bumping fpprec in a let makes (cdr *bfhalf*) evaluate to 1/2 · 2^-8 — e.g. atanh(bfloat(1/2)) then comes out exactly 256× too small. Those reads should go through a precision-aware accessor analogous to fpone (e.g. an fphalf/fpmhalf that returns the cached constant when its precision matches and recomputes 2^(fpprec-1) otherwise).