Elementary functions of balls

If x is inexact, prec is first lowered to a few bits beyond the accuracy its radius allows (relative for |x| < 1, absolute otherwise), and the radius is added to both outputs at the end (as |\sin'|, |\cos'| \le 1). A ball inside (-1, 1) whose radius exceeds a quarter of its midpoint gives [0 \pm 2^e] and 1 \pm 2^{2e-1} for |x| < 2^e; a wider ball elsewhere gives [0 \pm 1].

For a midpoint with z leading zero bits (|m| < 2^{-z}): if 2z \ge \mathit{prec} + 3 the results are m and 1, and if 4z \ge \mathit{prec} + 3 they are m - m^3/6 and 1 - m^2/2, with the Taylor remainders in the radius. Otherwise, for |m| < 1, the fixed-point kernel is called on |m| at \mathit{prec} + z bits plus guard bits (or a series for the reduced argument when z is large). For |m| \ge 1 the argument is reduced modulo \pi/2, with the quotient chosen so that the remainder is nonnegative, folded to [0, \pi/4], and the outputs are permuted and negated according to the quotient modulo 4. The reduction adds at most 3 ulps to the kernel's bound and is not repeated at a higher precision near zeros.

After the reduction |m| = a\pi/2 + sv, v \in [0, \pi/4], the result is s \tan v for even a and -s/\tan v for odd a (one Newton reciprocal). Near a zero or a pole v is small, so the kernel runs with the leading zero bits of v as extra precision, and the reduction is repeated at that precision when the first pass had fewer limbs. An exact m is never a multiple of \pi/2, and by the irrationality measure of \pi a b-bit m lies at least about 2^{-7.2 b} away from one; the repetition is capped at 8b + 1024 extra bits (status 0 beyond). For an inexact x no repetition is made once v falls below the radius.

The radius r of x enters through a pole test: with T enclosing \tan m and U = r(1 + |T|) \ge r/|\cos m|, U < 1 implies that the ball contains no pole, and then |\tan y - \tan m| \le r(1 + T^2)/(1 - U)^2.

The midpoint is dyadic, and its units bit and fraction bits give |m| = a/2 + st with a \in \{0, 1, 2, 3\}, s = \pm 1 and t \in [0, 1/4] exactly. Then v = \pi t enters the formulas of the reduction modulo \pi/2; t = 0 gives 0 or \pm 1 (and a pole of the tangent for odd a, status 0), and t = 1/4 the exact \tan(\pi/4) = 1. Otherwise t has z leading zero bits and the kernel runs at z extra bits. The radius r of x is a radius \pi r of the argument; \pi r \ge 1/4 gives [0 \pm 1] for the sine and cosine.

The fraction is reduced to lowest terms and folded by the symmetries to \pi t/q with 0 \le 2t \le q. The nonzero values of \sin(\pi t/q) and \cos(\pi t/q) are at least 1/q, so _mp_real_sin_cos_pi_ui_div_ui at n + \lceil(\operatorname{bits}(q) + 16)/\mathtt{FLINT\_BITS}\rceil limbs gives the relative accuracy. The tangent uses the closed forms \sqrt 3, 1/\sqrt 3, \sqrt 2 \mp 1, 2 \mp \sqrt 3 for q = 3, 6, 8, 12 and \sqrt{5 \mp 2\sqrt 5}, \sqrt{1 \mp 2/\sqrt 5} for q = 5, 10 (from 12 limbs), the Chebyshev sine and cosine and one division for the denominators where they win, and the tangent kernel otherwise.

The radius \rho of x enters as the factor [1 \pm (e^\rho - 1)] (exp), as \rho/(m - \rho) (log, by the mean value theorem), and as \rho, or \rho/(|m| - \rho)^2 when the ball stays beyond 1 in absolute value (atan). An exp ball too wide for two correct bits gives [0 \pm e^{m + \rho}].

• exp: |m| = q \log 2 + t, t \in [0, \log 2), with q from two limbs of 1/\log 2 and one multiply-subtract by \log 2; for m < 0 the quotient is rounded up and the remainder complemented, so that \exp(m) = 2^{-q}\exp(t) directly. Small |m| use 1 + m, 1 + m + m^2/2, ... or a series for the reduced argument.
• log: with m = 2^E u, u \in [1/2, 1), \log m = (E - 1)\log 2 + \operatorname{log1p}(2u - 1) up to 600 limbs and E \log 2 - (-\log u) by _mp_real_neglog_newton above; near 1, 2\operatorname{atanh}(s/(2 \pm s)) for s = |m - 1|, or the first terms of \log(1 + D).
• atan: \operatorname{atan}(m) = \pi/4 + \operatorname{atan}(s/(2+s)) for 1 < m < 2, s = m - 1, and \pi/2 - \operatorname{atan}(1/m) for m \ge 2; \operatorname{atan}(1) = \pi/4.

At few limbs the per-size kernels are called directly on the argument's limbs in place.

The domains are: 1 + y > 0 for log1p, |y| < 1 for atanh, and for asin, acos and acosh |y| \le 1 resp. y \ge 1 for an exact argument and strictly inside for an inexact one. mp_real_expm1_bits throws like mp_real_exp_bits for y \ge 2^{\mathtt{FLINT\_BITS} - 5}, and mp_real_sinh_cosh_bits for |y| \ge 2^{\mathtt{FLINT\_BITS} - 5}; expm1 of a large negative argument gives -1 with a tiny radius.

The functions are compositions of mp_real_exp_bits, mp_real_log_bits and mp_real_atan_bits by formulas without cancellation: \operatorname{log1p}(y) = \log(1 + y) with 1 + y formed exactly (or to \mathit{prec} + z bits for |y| < 2^{-z}); \operatorname{expm1} from \exp at \mathit{prec} + z bits; \sinh and \cosh from \operatorname{expm1}(\pm m); \tanh(m) = \operatorname{expm1}(2m)/(\operatorname{expm1}(2m) + 2); \operatorname{asin}(m) = \operatorname{atan}(m/\sqrt{(1 - m)(1 + m)}) and \operatorname{acos}(m) = 2\operatorname{atan}(\sqrt{(1 - m)/(1 + m)}) with 1 \pm m exact; \operatorname{asinh}(m) = \operatorname{log1p}(|m| + m^2/(1 + \sqrt{1 + m^2})); \operatorname{acosh}(m) = \operatorname{log1p}(d + \sqrt{d(m + 1)}) with d = m - 1 exact; and \operatorname{atanh}(m) = \operatorname{log1p}(2|m|/(1 - |m|))/2. Tiny and huge arguments use the leading terms with a bounded remainder (m, \log(2|m|), \pm 1), and \operatorname{asin}(\pm 1), \operatorname{acos}(\pm 1) and \operatorname{acosh}(1) are exact. The radius \rho enters through a bound for the variation over the ball, for example \rho \cosh(|m| + \rho) for \sinh and \rho/(1 - |m| - \rho) for atanh.

• mp_real_exp_notab_log2 reduces by \log 2 (computing and caching \log 2) and calls _mp_real_exp_notab; it is the fastest when \log 2 is already cached.
• mp_real_exp_notab_squaring needs no \log 2: \exp(|m|) = \exp(|m| 2^{-s})^{2^s} with |m| 2^{-s} < 1/2, the reduced value from _mp_real_exp_notab and s ball squarings (with s + 16 guard bits), and one reciprocal for m < 0. It is faster than the reduction when \log 2 must be computed.
• mp_real_exp_agm computes y_0 = \exp(m) to a fraction of the precision by the squarings and corrects it by one Newton–Taylor step on the AGM logarithm of y_0 (Elementary functions in mp_real), in O(M(n) \log n). It calls the squarings below about 30 limbs and overtakes them at about 15000 limbs.