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.