Fredrik Johansson

COMPUTING ELEMENTARY FUNCTIONS USING MULTI-PRIME ARGUMENT REDUCTION

We describe an algorithm for arbitrary-precision computation of the elementary functions (exp, log, sin, atan, etc.) which, after a cheap precomputation, gives roughly a factor-two speedup over previous state-of-the-art algorithms at precision from a few thousand bits up to millions of bits. Following an idea of Schönhage, we perform argument reduction using Diophantine combinations of logarithms of primes; our contribution is to use a large set of primes instead of a single pair, aided by a fast algorithm to solve the associated integer relation problem. We also list new, optimized Machin-like formulas for the necessary logarithm and arctangent precomputations. 1 M(B) is the complexity of B-bit multiplication. We can take M(B) = O(B log B) [HvdH21]. 2 The logarithmic and trigonometric functions have analogous algorithms; alternatively, they can be computed from the exponential via connection formulas and root-finding for inverses. For a more comprehensive overview of techniques for elementary function evaluation, see Smith [Smi89], Muller [Mul16], Brent and Zimmermann [BZ11], and Arndt [Arn10]. 1 3

In binary arithmetic, we only need to evaluate 3 c since multiplying by a power of two is free. This optimization is not a vital ingredient of the algorithm, however.

4 Unfortunately, the only published records of Schönhage's algorithm are two seminar talk abstracts which are light on details. The abstracts do mention the possibility of combining three primes instead of a single pair "for an improved design", but there is no hint of a practical algorithm working with arbitrarily large n, r and B, which will be presented here.

Introduction

There are two families of competitive algorithms for arbitrary-precision computation of elementary functions: the first uses Taylor series together with argument reduction and needs O(M(B) log 2 B) time for B-bit precision [START_REF] Brent | The complexity of multiple-precision arithmetic[END_REF], while the second is based on the arithmetic-geometric mean (AGM) iteration for elliptic integrals and achieves complexity O(M(B) log B) [START_REF] Brent | Fast multiple-precision evaluation of elementary functions[END_REF]. 1 Due to constant-factor overheads, optimized implementations of Taylor series tend to perform better than the AGM for practical sizes of B, possibly even for B in the billions.

The degree of argument reduction is a crucial tuning parameter in Taylor series methods. For example, the standard algorithm for the exponential function 2 amounts to choosing a reduction parameter r ≥ 0 and evaluating exp(x) = (exp(t)) 2 r , t = x/2 r .

(1)

If |x| < 1, this costs r squarings plus the summation of N ≤ B/r terms of the series for exp(t), or better N/2 terms for s = sinh(t) (with exp(t) = s + √ s 2 + 1). In moderate precision (up to around B = 10 4 ) the series evaluation costs O( √ N ) multiplications; for quasilinear complexity as B → ∞, we use the "bit-burst algorithm": we write exp(t) = exp(t 1 ) • exp(t 2 ) • • • where t j extracts 2 j bits in the binary expansion of t and evaluate each exp(t j ) series using binary splitting.

Asymptotically, r should grow at most logarithmically with B, or the O(rM(B)) time spent on squarings will dominate. In practice, the best r will be of order 10 to 100 (varying with B) and these r squarings may account for a large fraction of the work to evaluate the function. This prompts the question: can we reduce the argument to size 2 -r without the cost of r squarings?

The only known solution relies on precomputation. For example, we need only a single multiplication for r-bit reduction if we have a precomputed table of exp(j/2 r ), 0 ≤ j < 2 r , or m multiplications with an m-partite table of m2 r/m entries. Tables of this kind are useful up to a few thousand bits [START_REF] Johansson | Efficient implementation of elementary functions in the medium-precision range[END_REF], but they are rarely used at higher precision since they yield diminishing returns as the space and precomputation time increases linearly with B and exponentially with r. Most commonly, arbitrary-precision software will only cache higher-precision values of the constants π and log(2) computed at runtime, used for an initial reduction to ensure |x| < 1.

Schönhage's method. In 2006, Schönhage [START_REF] Schönhage | Fast algorithms for computing exp, ln, sin[END_REF][START_REF] Schönhage | Fast medium precision computation of elementary functions[END_REF] proposed a method to compute elementary functions using "diophantine combinations of incommensurable logarithms" which avoids the problem with large tables. The idea is as follows: given a real number x, we determine integers c, d such that

x ≈ c log(2) + d log(3) (2)
within some tolerance 2 -r (it is a standard result in Diophantine approximation that such c, d exist for any r). We can then use the argument reduction formula exp(x) = exp(t) 2 c 3 d , t = x -c log(2) -d log(3).

(3)

There is an analogous formula for complex x and for trigonometric functions using Gaussian primes. The advantage of Schönhage's method is that we only need to precompute or cache the two constants log(2) and log(3) to high precision while the rational power product 2 c 3 d can be computed on the fly using binary exponentiation. If 3 c < 2 B , this step costs O(M(B)). 3 Schönhage seems to have considered this method useful only for B in the range from around 50 to 3000 bits (in his words, "medium precision"). The problem is that the coefficients c, d in (3) grow exponentially with the desired amount of reduction. Indeed, solutions with |t| < 2 -r will generally have c, d = O(2 r/2 ). It is also not obvious how to compute the coefficients c and d for a given x; we can use a lookup table for small r, but this retains the exponential scaling problem.

Our contribution. In this work, we describe a version of Schönhage's algorithm in which we perform reduction using a basis of n primes, where n is arbitrary and in practice may be 10 or more. The coefficients (power-product exponents) will then only have magnitude around O(2 r/n ), allowing much greater reduction than with a single pair of primes. 4 Section 2 presents an algorithm for quickly finding an approximating linear combination of several logarithms, which is a prerequisite for making the method practical. Section 3 describes the main algorithm for elementary functions in more detail. Section 4 discusses use of Machin-like formulas for fast precomputation of logarithms or arctangents, where we tabulate new optimized multi-evaluation formulas for special sets of values.

Our implementation results presented in section 5 show that the new version of Schönhage's algorithm scales remarkably well: we can quickly reduce the argument to magnitude 2 -r where we may have r ≥ 100 at moderately high precision (a few thousand bits) and perhaps r ≥ 500 at millions of bits. When n is chosen optimally, the new algorithm runs roughly twice as fast as the best previous elementary function implementations (both Taylor and AGM-based) for bit precisions B from a few thousand up to millions. The storage requirements (nB bits) and precomputation time (on par with one or a few extra function evaluations) are modest enough that the method is ideal as a default algorithm in arbitrary-precision software over a large range of precisions.

Historical note. With the exception of Schönhage's work, we are not aware of any previous investigations into algorithms of this kind for arbitrary-precision computation of elementary functions of real and complex arguments. However, the underlying idea of exploiting differences between logarithms of prime numbers in a computational setting goes back at least to Briggs' 1624 Arithmetica logarithmica [START_REF] Briggs | Arithmetica logarithmica[END_REF][START_REF] Roegel | A reconstruction of the tables of Briggs[END_REF]. Briggs used a version of this trick when extending tables of logarithms of integers. We revisit this topic in section 4.

Integer relations

We consider the following inhomogeneous integer relation problem: given real numbers x and α 1 , . . . , α n and a tolerance 2 -r , find a vector (c 1 , . . . , c n ) ∈ Z n with small coefficients such that

x ≈ c 1 α 1 + . . . c n α n (4)
with error at most 2 -r . We assume that the equation c 1 α 1 + . . . + c n α n = 0 has no solution over the integers. In the special case where P = {p 1 , . . . p n } is a set of prime numbers and α i = log(p i ), solving (4) will find a P -smooth rational approximation

exp(x) ≈ p c1 1 • • • p cn n ∈ Q (5)
with small numerator and denominator. Integer relation problems can be solved using lattice reduction algorithms like LLL [START_REF] Arjen K Lenstra | Factoring polynomials with rational coefficients[END_REF][START_REF] Cohen | A Course in Computational Algebraic Number Theory[END_REF]. However, directly solving

c 0 x + c 1 α 1 + . . . + c n α n ≈ 0 (6)
will generally introduce a denominator c 0 = 1, requiring a c 0 -th root extraction on the right-hand side of (5). In any case, running LLL each time we want to evaluate an elementary function will be too slow. Algorithm 1 solves these issues by precomputing solutions to the homogeneous equation c 1 α 1 + . . . c n α n ≈ 0 and using these relations to solve the inhomogeneous version (4) through iterated reduction.

Analysis of Algorithm 1. We assume heuristically that each step in the precomputation phase (1) succeeds to find a relation (7) with ε i within a small factor of ±C -i and with coefficients (d i,1 , . . . , d i,n ) of magnitude O(C i/n ). We will simply observe that this always seems to be the case in practice; a rigorous justification would require further analysis.

It can happen that picking the first integer relation computed by LLL yields the same relation consecutively (ε i = ε i+1 ). In that case, we can just pick a different relation (while keeping the ε i sorted) or skip the duplicate relation. However, a Algorithm 1 Approximate x ∈ R to within 2 -r by a linear combination

x ≈ c 1 α 1 + . . . + c n α n , c i ∈ Z
given α i ∈ R which are linearly independent over Q. Alternatively, find a good approximation subject to some size constraint f (c 1 , . . . , c n ) ≤ M .

(1) Precomputation (independent of x): choose a real convergence factor C > 1.

For i = 1, 2, . . ., LLL-reduce

     1 0 . . . 0 ⌊C i α 1 + 1 2 ⌋) 0 1 . . . 0 ⌊C i α 2 + 1 2 ⌋) . . . . . . . . . . . . . . . 0 0 . . . 1 ⌊C i α n + 1 2 ⌋)      .
This yields an approximate integer relation

ε i = d i,1 α 1 + . . . d i,n α n , ε i = O(C -i ). (7)
(In fact, it yields n such relations; we can choose any one of them.) We store tables of the coefficients d i,j and floating-point approximations of the errors ε i . We stop after the first i where

|ε i | < 2 -r .
(2) Reduction (given x).

• Let (c 1 , . . . , c n ) = (0, . . . , 0).

• For i = 1, 2, . . ., compute m i = ⌊x/ε i + 1/2⌋ and update:

(c 1 , . . . , c n ) ← (c 1 + m i d i,1 , . . . , c n + m i d i,n ),
x ← x -m i ε i . Stop and return the relation (c 1 , . . . , c n ) when |x| < 2 -r or when the next update will give f (c 1 , . . . , c n ) > M .

decrease by much more than a factor C between successive step should be avoided as it will result in larger output coefficients.

Phase (1) terminates when i = N ≈ r log(2)/ log(C). The multiplier m i computed in each step of the phase (2) reduction has magnitude around C. The coefficients (c 1 , . . . , c n ) at the end of phase (2) will therefore have magnitude around

N i=1 C i/n+1 = C 1/n+1 C 1/n -1 2 r/n -1 ≈ Cn log(C) 2 r/n (8)
or perhaps a bit smaller than this since on average there can be some cancellation. The prefactor C/ log(C) is minimized when C = e, or in other words it is theoretically optimal to force ε i = Θ(exp(-i)). However, this prefactor does not vary strongly with C, and a choice like C = 2 (one bit per step) or C = 10 (one decimal per step) may be convenient.

Step i of phase (1) requires LLL-reducing a matrix with β-bit entries where β = O(i). The standard complexity bound for LLL is O(n 5+ε β 2+ε ), so phase (2) costs O(n 5+ε r 3+ε ). 5 In our application, the tables generated in phase (1) are small (a few kilobytes) and do not need to be generated at runtime, so it suffices to note that the computations are feasible for ranges of n and r of interest; for empirical results, see section 5.

Phase (2) requires O(nr) arithmetic operations with r-bit numbers, for a running time of O(nr 2+ε ). It is convenient to treat x and α i as fixed-point numbers with r-bit fractional part. As an optimization, we can work with a machine-precision (53-bit) floating-point approximations of x ′ and the errors ε i . We periodically recompute

x ′ = x -(c 1 α 1 + . . . + c n α n ) accurately from the full-precision values only when this approximation runs out of precision, essentially every 53/ log 2 (C) steps. The resulting algorithm has very low overhead. We will not consider asymptotic complexity improvements since r will be moderate (a small multiple of the word size) in our application.

Computation of elementary functions

Given x ∈ R and a set of prime numbers P = {p 1 , p 2 , . . . , p n }, the algorithm described in the previous section allows us to find integers c 1 , . . . , c n such that

t = x -(c 1 log(p 1 ) + . . . + c n log(p n )) (9)
is small, after which we can evaluate the real exponential function as

exp(x) = exp(t) p c1 1 • • • p cn n . (10)
Algorithm 2 describes the procedure in some more detail.

Algorithm 2 Computation of exp(x) for x ∈ R to B-bit precision using argument reduction by precomputed logarithms of primes.

(1) Precomputation (independent of x): select a set of prime numbers P = {p 1 , p 2 , . . . , p n } with p 1 = 2. Compute log(p 1 ), . . . , log(p n ) to B-bit precision.

(2) Using Algorithm 1, find an integer relation x ≈ c 1 log(p 1 ) + . . . + c n log(p n ), attempting to make the error as small as possible subject to c 1 , . . . , c n P ≤ B. This step can use low precision (about r bits where 2 -r is the target reduction, in practice no more than a few machine words). (3) Compute the power product v/w = p c2 2 • • • p cn n as an exact fraction, using binary splitting to recursively split the set of primes in half and using binary exponentiation to compute the individual powers. (4) Calculate t = x -(c 1 log(p 1 ) + . . . + c n log(p n )) using the precomputed logarithms. (5) Compute u = exp(t) using Taylor series: depending on B, either use rectangular splitting for the sinh series or use the bit-burst decomposition t = t 1 +t 2 +. . . with binary splitting (see e.g. [START_REF] Brent | Modern Computer Arithmetic[END_REF] for details). (6) Return 2 c1 uv/w.

Remarks. The bottleneck in the argument reduction is the cost of evaluating the power product p c1 1 • • • p cn n ∈ Q. How large coefficients (exponents) should we allow? A reasonable heuristic, implemented in Algorithm 2, is to choose coefficients such that the weighted norm

ν = c 1 , . . . , c n P = n i=1 pi =2 |c i | log 2 (p i ) (11)
is smaller than B: this ensures that the rational power product p c1 1 • • • p cn n has numerator and denominator bounded by B bits. We discount the prime 2 in the norm with the assumption that we factor out powers of two when performing binary arithmetic. If |x| > 1, we should use log(2) alone for the first reduction in Algorithm 1 so that the corresponding exponentiation is free.

We note that when computing the power product, there is no need to compute GCDs since the numerator and denominator are coprime by construction.

There is not much to say about numerical issues; essentially, we need about log 2 ( i |c i | log(p i )) guard bits to compensate for cancellation in the subtraction, which in practice always will be less than one extra machine word. If |x| ≫ 1, we need an additional log 2 |x| guard bits for the accurate removal of log(2).

3.1. Numerical example. We illustrate computing exp(x) to 10000 digits (or B = 33220 bits) where x = √ 2 -1, using n = 13 primes. The following Pari/GP output effectively shows the precomputations of phase (1) of Algorithm 1 with convergence rate C = 10. Since 2 -100 ≈ 7.9 • 10 -31 , reducing by 32 relations with C = 10 is equivalent to r = 100 squarings in (1).

We prepend the relation [1, 0, . . .] for an initial reduction by log(2), and we can eliminate the duplicate entries.

The We compute the reduced argument t = x -log(2 c1 v/w) ≈ -1.57 • 10 -32 by subtracting a linear combination of precomputed logarithms. Now taking 148 terms of the Taylor series for sinh(t) yields an error smaller than 10 -10000 . Evaluating this Taylor series using rectangular splitting costs roughly 2 √ 148 ≈ 24 full 10000-digit multiplications, and this makes up the bulk of the time in the exp(x) evaluation.

For comparison, computing exp(x) using (1) without precomputation, it is optimal to perform r ≈ 20 squarings after which we need 555 terms of the sinh series, for a cost of r+2 √ 555 ≈ 67 multiplications.6 Alternatively, computing log(x) using the AGM requires 25 iterations, where each iteration a n+1 , b n+1 = (a n + b n )/2, √ a n b n costs at least as much as two multiplications.

Counting arithmetic operations alone, we can thus expect Algorithm 2 to be at least twice as fast as either method in this example. As we will see in section 5, this back-of-the-envelope estimate is quite accurate.

Trigonometric functions.

We can compute the real trigonometric functions via the exponential function of a pure imaginary argument, using Gaussian primes a + bi ∈ Z[i] for reduction. Enumerated in order of norm a 2 + b 2 , the nonreal Gaussian primes are 1 + i, 2 + i, 3 + 2i, 4 + i, 5 + 2i, 6 + i, 5 + 4i, 7 + 2i, 6 + 5i, . . .

where we have discarded entries that are equivalent under conjugation, negation or transposition of real and imaginary parts (we choose here, arbitrarily, the representatives in the first quadrant and with a ≥ b).

The role of the logarithms log(p) is now assumed by the irreducible angles

α = 1 i [log(a + bi) -log(a -bi)] = 2 atan b a (13)
which define rotations by e iα = (a + bi)/(a -bi) on the unit circle. We have the argument reduction formula

cos(x) + i sin(x) = exp(ix) = exp(i(x -cα)) (a + bi) c (a -bi) c , c ∈ Z (14)
which can be iterated over a combination of Gaussian primes. Algorithm 3 computes cos(x) and sin(x) together using this method.

Algorithm 3 Computation of cos(x) + i sin(x) = exp(ix) for x ∈ R to B-bit precision using argument reduction by precomputed irreducible angles.

(1) Precomputation (independent of x): select a set of Gaussian prime numbers

Q = {a 1 + b 1 i, . . . , a n + b n i} from (12) with a 1 + b 1 = 1 + i. Compute 2 atan(b 1 /a 1 ), . . . , 2 atan(b n /a n ) to B-bit precision.
(2) Using Algorithm 1, find an integer relation x ≈ c 1 2 atan(b 1 /a 1 ) + . . . + c n 2 atan(b n /a n ), attempting to make the error as small as possible subject to c 1 , . . . , c n Q ≤ B. This step can use low precision (about r bits where 2 -r is the target reduction, in practice no more than a few machine words).

(3) Compute the power product

v w = (a 2 + b 2 i) c2 • • • (a n + b n i) cn (a 2 -b 2 i) c2 • • • (a n -b n i) cn ∈ Q(i) (15)
as an exact fraction, using binary splitting to recursively split the set of primes in half and using binary exponentiation to compute the individual powers. (4) Calculate t = x -(c 1 2 atan(b 1 /a 1 ) + . . . + c n 2 atan(b 1 /a 1 )) using the precomputed arctangents. (5) Compute u = exp(it) using Taylor series (depending on B, either using rectangular splitting for the sin series or using the bit-burst decomposition t = t 1 + t 2 + . . . with binary splitting). (6) Return i c1 uv/w.

Remarks. Here, a suitable norm is

ν = c 1 , . . . , c n Q = n j=1 pj =1+i |c j | log 2 (a 2 j + b 2 j ). ( 16
)
The special prime 2 in the argument reduction for the real exponential is here replaced by the Gaussian prime 1 + i, for which

(1 + i) c (1 -i) c = i c (17)
can be evaluated in constant time; the angle reduction corresponds to removal of multiples of π/2. We only need to compute the factors in the numerator of the right-hand side of (15) since the remaining product can be obtained via complex conjugation. As in the real case, all factors are coprime so we can multiply numerators and denominators using arithmetic in Z[i] without the need for GCDs.

We can save a marginal amount of work (essentially in the last division) if we want either the sine of the cosine alone, or if we want tan(x).

3.3. Inverse functions. The formulas above can be transposed to compute the inverse functions. For example,

log(x) = log x p c1 1 • • • p cn n + (c 1 log(p 1 ) + . . . + c n log(p n )). ( 18
)
For the complex logarithm or arctangent, we need to be careful about selecting the correct branches.

As an alternative, we recall the standard method of implementing the inverse functions using Newton iteration, starting from an low-precision approximation obtained with any other algorithm. The constant-factor overhead of Newton iteration can be reduced with an m-th order method derived from the addition formula for the exponential function [Arn10, section 32.1]. If y = log(x) + ε, then log(x) = y + log(1 + δ), δ = x exp(-y) -1.

(

) 19
We first compute y ≈ log(x) at precision B/m (calling the same algorithm recursively until we hit the basecase range) so that the unknown error ε is O(2 -B/m ). Then, we evaluate (19) at precision B using the Taylor series for log(1+δ) truncated to order O(δ m ). This gives us log(x) with error O(2 -B ).

The inverse trigonometric functions can be computed analogously via the arctangent: if y = atan(x) + ε, then

atan(x) = y + atan(δ), δ = x -t 1 + tx = cx -s c + sx , t = tan(y) = s c = sin(y) cos(y) . (20)
With a suitably chosen m (between 5 and 15, say) and rectangular splitting for the short Taylor series evaluation, the inverse functions are perhaps 10%-30% more expensive than the forward functions with this method.

Precomputation of logarithms and arctangents

The precomputation of logarithms and arctangents of small integer or rational arguments is best done using binary splitting evaluation of trigonometric and hyperbolic arctangent series

atan 1 x = ∞ k=0 (-1) k (2k + 1) 1 x 2k+1 , atanh 1 x = ∞ k=0 1 (2k + 1) 1 x 2k+1 . ( 21
)
We want the arguments x in (21) to be integers, and ideally large integers so that the series converge rapidly. It is not a good idea to use the primes p or Gaussian integer tangents b/a directly as input since convergence will be slow; it is better to recycle values and evaluate differences of arguments (Briggs' method). For example, if we have already computed log(2), we can compute logarithms of successive primes using [START_REF] Gourdon | The logarithmic constant: Log 2[END_REF]

log(p) = log(2) + 1 2 log p -1 2 + log p + 1 2 + atanh 1 2p 2 -1 . ( 22
)
Methods to reduce arctangents to sums of more rapidly convergent arctangent series have been studied by Gauss, Lehmer, Todd and others [Leh38, Tod49, Wet96]. The prototype is Machin's formula Given the first n primes, we will thus look for a set of integers X = {x 1 , x 2 , . . . , x n }, as large as possible, such that there is an integer relation

   log(p 1 ) . . . log(p n )    = M    2 atanh(1/x 1 ) . . . 2 atanh(1/x n )    , M ∈ Q n×n ( 24
)
or similarly (with different X and M ) for Gaussian primes

   atan(b 1 /a 1 ) . . . atan(b n /a n )    = M    atan(1/x 1 ) . . . atan(1/x n )    , M ∈ Q n×n . ( 25
)
For example, the primes P = {2, 3} admit the simultaneous Machin-like formulas log(2) = 4 atanh(1/7) + 2 atanh(1/17), log(3) = 6 atanh(1/7) + 4 atanh(1/17), i.e.

X = {7, 17}, M = 2 1 3 2 .
The following method to find relations goes back to Gauss who used it to search for generalizations of Machin's formula. Arndt [Arn10, section 32.4] also discusses the application of simultaneous computation of logarithms of several primes.

The search space for candidate sets X in ( 24) and ( 25) is a priori infinite, but it can be narrowed down as follows. Let P = {p 1 , . . . , p n }. Since 2 atanh(1/x) = log(x + 1) -log(x -1) = log

x + 1

x -1 , we try to write each p ∈ P as a power-product of P -smooth rational numbers of the form (x + 1)/(x -1). We will thus look for solutions X of (25) of the form

X ⊆ Y, Y = {x : x 2 -1 is P -smooth}, (26)
i.e. such that both x + 1 and x -1 are P -smooth. Similarly, we look for solutions of (25) of the form

X ⊆ Z, Z = {x : x 2 + 1 is Q-smooth} ( 27
)
where Q is the set of norms {a 2 1 + a 2 1 , . . . , a 2 n + b 2 n }. It is a nontrivial fact that the sets Y and Z are finite for each fixed set of primes P or Q. For the 25 first primes p < 100, the set Y has 16223 elements which have been tabulated by Luca and Najman [START_REF] Luca | On the largest prime factor of x 2 -1[END_REF][START_REF] Luca | Errata to "On the largest prime factor of x 2 -1[END_REF]; the largest element 7 is x = 19182937474703818751 with

x -1 = 2 • 5 5 • 11 • 19 • 23 2 • 29 • 59 4 • 79, x + 1 = 2 22 • 3 • 17 3 • 37 • 41 • 43 • 67 • 71.
For the first 22 Gaussian primes, having norms a 2 + b 2 < 100, the set Z has 811 elements which have been tabulated by Najman [START_REF] Najman | Smooth values of some quadratic polynomials[END_REF]; the largest element is x = 69971515635443 with

x 2 + 1 = 2 • 5 5 • 17 • 37 • 41 2 • 53 2 • 89 • 97 3 • 137 2 • 173.
7 Knowing this upper bound, the Luca-Najman table can be reproduced with a brute force enumeration of 97-smooth numbers x -1 with x ≤ 19182937474703818751, during which one saves the values x for which trial division shows that x + 1 is 97-smooth. This computation takes two hours on a 2022-era laptop. Reproducing the table Z takes one minute.

Given a candidate superset Y = {y 1 , . . . , y r } or Z = {z 1 , . . . , z s }, we can find a formula X with large entries using linear algebra:

• Let X = {}, and let R be an initially empty (0 × n) matrix.

• For x = y r , y r-1 , . . . or x = z r , z s-1 , . . . in order of decreasing magnitude, let E = (e 1 , . . . , e n ) be the vector of exponents in the factorization of the rational number

(x + 1)/(x -1) = p e1 1 • • • p en n ,
respectively,

x 2 + 1 = (a 2 1 + b 2 1 ) e1 , . . . , (a 2 n + b 2 n ) en .
• If E is linearly independent of the rows of R, add x to X and adjoin the row E to the top of R; otherwise continue with the next candidate x. • When R has n linearly independent rows, we have found a complete basis X and the relation matrix is given by M = R -1 .

Tables 1 and2 give the Machin-like formulas found with this method using the exhaustive Luca-Najman tables for Y and Z. We list only the set X since the matrix M is easy to recover with linear algebra (in fact, we can recover it using LLL without performing any factorization). The corresponding Lehmer measure µ(X) = x∈X 1/ log 10 (|x|) gives an estimate of efficiency (lower is better). 4.2. Remarks about the tables. We conjecture that the formulas in Tables 1 and 2 are the best possible (in the Lehmer sense) n-term formulas for the respective sets of n primes or Gaussian primes.

Apart from the first few entries which are well known, we are not aware of a previous tabulation of this kind. There is an extensive literature about Machin-like formulas for computing π alone, but little about computing several arctangents simultaneously. There are some preexisting tables for logarithms, but they are not optimal. Arndt [START_REF] Arndt | Matters Computational: ideas, algorithms, source code[END_REF] gives a slightly less efficient formula for the 13 primes up to 41 with µ(X) = 1.48450, which appears to have been chosen subject to the constraint max(X) < 2 32 . Gourdon and Sebah [START_REF] Gourdon | The logarithmic constant: Log 2[END_REF] give a much less efficient formula for the first 25 primes derived from (22), with µ(X) > 7.45186.

The claim that the formulas in Tables 1 and2 are optimal comes with several caveats. We can achieve lower Lehmer measures if we add more arctangents. Indeed, the formula for P = {2, 3, 5, 7} has a lower Lehmer measure than the formulas for {2}, {2, 3} and {2, 3, 5}, so we may just as well compute four logarithms if we want the first one or three. A more efficient formula for log(2) alone is the threeterm X = {26, 4801, 8749} with µ(X) = 1.23205 which however cannot be used to compute log(3), log(5) or log(7) (the set X 2 -1 is 7-smooth but does not yield a relation for either 3, 5 or 7). The 1-term formula for atan(1) = π/4 has infinite Lehmer measure while Machin's formula (23), which follows from the 13-smooth factorizations 5 2 + 1 = 2 • 13 and 239 2 + 1 = 2 • 13 4 , achieves µ(X) = 1.85112.

In practice µ(X) is not necessarily an accurate measure of efficiency: it overestimates the benefits of increasing x, essentially because the running time in binary splitting tends to be dominated by the top-level multiplications which are independent of the number of leaf nodes. It is therefore likely an advantage to keep the number of arctangents close to n.

A curiosity is that in the logarithm relations, we have det(R) = ±1 and therefore M ∈ Z n×n for the first 21 sets of primes P , but for P containing the primes up to 79, 83, 89 and 97 respectively the determinants are -2, -6, -4 and -4.

Implementation results

The algorithms have been implemented in Arb [START_REF] Johansson | Arb: Efficient arbitrary-precision midpoint-radius interval arithmetic[END_REF] version 2.23. The following results were obtained with Arb 2.23 linked against GMP 6.2.1 [GMP21], MPFR 4.1.0 [FHL + 07], and FLINT 2.9 [START_REF] Hart | Fast Library for Number Theory: An Introduction[END_REF], running on an AMD Ryzen 7 PRO 5850U (Zen3).

5.1. Default implementations with fixed n. Previously, all elementary functions in Arb used Taylor series with precomputed lookup tables up to B = 4608 bits. The tables are m-partite giving r-bit reduction with r ≤ 14 and m ≤ 2, requiring 236 KB of fixed storage [START_REF] Johansson | Efficient implementation of elementary functions in the medium-precision range[END_REF]. At higher precision, the previous implementations used argument reduction based on repeated argument-halving (requiring squaring or square roots) together with rectangular splitting or bit-burst evaluation of Taylor series, with the exception of log which wrapped the AGM-based logarithm in MPFR. To the author's knowledge, these were the fastest arbitrary-precision implementations of elementary functions available in public software libraries prior to this work.

In Arb 2.23, all the elementary functions were rewritten to use the new algorithm with the fixed number n = 13 of primes, starting from a precision between B = 2240 bits (for exp) and B = 3400 bits (for atan) up to B = 4000000 bits (just over one million digits). The Newton iterations ( 19) and ( 20) are used to reduce log and atan to the exponential and trigonometric functions. The B-bit precomputations of logarithms and arctangents are done at runtime using the n = 13 Machin-like formulas of Table 1 andTable 2.

We compare timings for the old and new implementations in Table 3.

Remarks. The average speedup is around a factor two (1.3× to 2.4×) over a large range of precisions. The typical slowdown for a first function call is also roughly a factor two, i.e. the precomputation takes about as long as a single extra function call. 8 This is clearly a worthwhile tradeoff for most applications; e.g. for a numerical integration b a f (x)dx where the integrand f will be evaluated many times, we do observe a factor-two speedup in the relevant precision ranges.

The relatively large speedup for atan is explained by the fact that the traditional argument reduction method involves repeated square roots which are a significant constant factor more expensive than the squarings for exp.

The relatively small speedup for sin and cos is explained by the fact that traditional argument reduction method only requires real squarings (via the half-angle formula for cos), while the new method uses complex arithmetic.

Previously, the AGM-based logarithm was neck and neck with the Taylor series for exp at any precision (these algorithms were therefore roughly interchangeable if one were to use Newton iteration to compute one function from the other). With the new algorithm, Taylor series have a clear lead.

The default parameter n = 13 was chosen to optimize performance around a few thousand digits, this range being more important for typical applications than Since the tables are small (a few KiB) and independent of B, they can be precomputed once and for all, so the timings (here essentially just exercising FLINT's LLL implementation) are not really relevant. Indeed, in the previously discussed default implementation of elementary functions, the n = 13 tables are stored as trivial technical restriction which we do not bother to lift since there would be a pay-off only for multi-million B. static arrays written down in the source code. However, the timings are reasonable enough that tables could be generated at runtime in applications that will perform a large number of function evaluations.

5.3. Function evaluation with variable n and B. Table 5 shows timings for the computation of the exponential function and trigonometric functions for different combinations of precision B and number of primes n. The n = 0 reference timings correspond to the old algorithm without precomputation, in which repeated squaring will be used instead. At lower precisions, using 10-20 primes seems to be optimal. It is interesting to note that roughly a factor-two speedup can be achieved across a huge range of precision when n increases with B. It seems likely that n = 128 or more primes could be useful at extreme precision, though the precomputation will increase proportionally.

We used the Machin-like formulas from Table 1 and Table 2 only up to n = 25 or n = 22; for n = 32 and n = 64 we fall back on less optimized formulas, which results in a noticeably slower precomputation.

Extensions and further work

We conclude with some ideas for future research.

6.1. Complexity analysis and fine-tuning. It would be interesting to perform a more precise complexity analysis. Under some assumptions about the underlying arithmetic, it should be possible to obtain a theoretical prediction for the optimal number of primes n as a function of the bit precision B, with an estimate of the possible speedup when B → ∞.

There are a large number of free parameters in the new algorithm (the number of primes n, the choice of primes, the precise setup of the precomputed relations ε i , the allowed size of the power product, choices in the subsequent Taylor series evaluation...). Timings can fluctuate depending on small adjustments to these parameters and with different values of the argument x. It is plausible that a consistent speedup can be obtained by simply tuning all the parameters more carefully. 6.2. Complex arguments. All elementary functions of complex arguments can be decomposed into real exponentials, logarithms, sines and arctangents after separating real and imaginary parts. An interesting alternative would be to compute exp(z) or log(z) directly over C, reducing z with respect to complex lattices generated by pairs of Gaussian primes. We do not know whether this presents any advantages over separating the components. 6.3. p-adic numbers. The same methods should work in the p-adic setting. For the p-adic exponential and logarithm, we can choose a basis of n prime numbers p i = p and use LLL to precompute relations n i=1 c i log(p i ) = O(p i ) for i = 1, 2, . . . , r. We can then use these relations to reduce the argument to order O(p r ) before evaluating the function using Taylor series or the p-adic bit-burst method [START_REF] Caruso | Fast evaluation of some p-adic transcendental functions[END_REF]. We have not attempted to analyze or implement this algorithm.

6.4. More Machin-like formulas. It would be useful to have larger tables of optimized Machin-like formulas for multi-evaluation of logarithms and arctangents. In practice, formulas need not be optimal as long as they are "good enough"; for example, one could restrict the search space to 64-bit arctangent arguments x. Nevertheless, a large-scale computation of theoretically optimal tables would be an interesting challenge of its own.

  Simultaneous Machin-like formulas. If we have the option of computing the set of values log(p 1 ), . . . , log(p n ) or atan(b 1 /a 1 ), . . . , atan(b n /a n ) in any order (not necessarily one by one), then we can try to look for optimized simultaneous Machin-like formulas[START_REF] Arndt | Matters Computational: ideas, algorithms, source code[END_REF].

  19 463 • 37 634 2 274 • 3 414 • 5 187 • 7 314 • 11 211 • 17 392 • 23 36 • 29 369 • 31 231 where the numerator and denominator have 7679 and 7678 bits, comfortably smaller than B.

	or	phase (2) reduction in Algorithm 1 with x = [-274, -414, -187, -314, -211, 651, -392, 463, -36, -369, -231, 634, 0] √ 2 -1 now yields the relation
		exp(x) ≈	2 c1 v w	=	13 651 •

Table 1 .

 1 n-term Machin formulas {atanh(1/x) : x ∈ X} for simultaneous computation of log(p) for the first n primes p ∈ P .

	n	P	X

Table 2 .

 2 n-term Machin formulas {atan(1/x) : x ∈ X} for simultaneous computation of the irreducible angles atan(b/a) for the first n nonreal Gaussian primes a + bi, having norms a 2 + b 2 ∈ Q.

	n	Q	X

Table 3 .

 3 Time to compute elementary functions to D decimal digits (B ≈ 3.32D) with Arb 2.23. Old is the time in seconds with the new algorithm disabled. New is the time in seconds with the new algorithm enabled, using the fixed default number n = 13 of primes. First is the time for a first function call, and Repeat is the time for repeated calls (with logarithms and other data already cached). We show average timings for 100 uniformly random input x ∈ (0, 2).

			exp(x)	log(x)	(cos(x), sin(x))	atan(x)
	D		First	Repeat	First	Repeat	First	Repeat	First	Repeat
	1000	Old	2.92e-05 2.91e-05 0.000145 3.69e-05 3.49e-05 3.49e-05 3.52e-05 3.52e-05
		New	0.000182 2.04e-05 0.000188 2.58e-05 0.00019 2.84e-05 3.52e-05 3.52e-05
		Speedup	0.16×	1.43×	0.77×	1.43×	0.18×	1.23×	1.00×	1.00×
	2000	Old	0.000103 0.000101 0.000367 0.000110 0.000217 9.92e-05 0.000423 0.000217
		New	0.000480 4.9e-05 0.000500 6.07e-05 0.000542 7.92e-05 0.000564 9.83e-05
		Speedup	0.22×	2.06×	0.73×	1.81×	0.40×	1.25×	0.75×	2.21×
	4000	Old	0.000355 0.000353 0.00103 0.000348 0.000511 0.000341 0.000915 0.000660
		New	0.00107 0.000149 0.00111 0.000187 0.00119 0.000211 0.00124 0.000269
		Speedup	0.33×	2.37×	0.93×	1.86×	0.43×	1.62×	0.74×	2.45×
	10000	Old	0.00185 0.00168 0.00439 0.00166	0.0022	0.00177 0.00323 0.00272
		New	0.00384 0.000826 0.00418 0.000977 0.00417 0.000935 0.00461 0.00122
		Speedup	0.48×	2.03×	1.05×	1.70×	0.53×	1.89×	0.70×	2.23×
	100000	Old	0.0541	0.0536	0.143	0.0632	0.0880	0.0818	0.0957	0.0896
		New	0.107	0.0354	0.114	0.0377	0.129	0.0509	0.140	0.0586
		Speedup	0.51×	1.52×	1.25×	1.68×	0.68×	1.61×	0.68×	1.53×
	1000000	Old	1.10	1.09	2.84	1.36	1.66	1.61	2.02	1.97
		New	2.18	0.864	2.31	0.982	2.83	1.25	3.02	1.58
		Speedup	0.51×	1.26×	1.23×	1.39×	0.59×	1.29×	0.67×	1.25×

Table 4 .

 4 Static precomputation of reduction tables: phase (1) of Algorithm 1.

	α1, . . . , αn	n	Smallest εi	Max r	Data	Time
	Logarithms	2	ε7 = +1.82 • 10 -5	15	0.2 KiB 0.0000514 s
		4	ε11 = -1.46 • 10 -14	45	0.3 KiB	0.000228 s
		8	ε33 = +7.66 • 10 -33	106	1.1 KiB	0.00249 s
		16	ε67 = +5.18 • 10 -71	233	3.2 KiB	0.0447 s
		32 ε144 = -1.51 • 10 -141	467	11 KiB	1.24 s
		64 ε268 = -4.42 • 10 -266	881	38 KiB	34.2 s
	Arctangents	2	ε7 = -4.75 • 10 -5	14	0.2 KiB 0.0000472 s
		4	ε14 = -2.95 • 10 -15	48	0.4 KiB	0.000248 s
		8	ε33 = +6.43 • 10 -33	106	1.1 KiB	0.00256 s
		16	ε64 = +1.77 • 10 -71	235	3.0 KiB	0.0448 s
		32 ε143 = +1.70 • 10 -140	464	11 KiB	1.22 s
		64 ε270 = +1.42 • 10 -267	886	38 KiB	34.6 s

Table 5 .

 5 Computation of the exponential function and the trigonometric functions. The argument is taken to be x = √ 2 -1. Precomp is the time (in seconds) to precompute n logarithms or arctangents for use at B-bit precision. The cached logarithms or arctangents take up Data space. Time is the time to evaluate the function once this data has been precomputed. The argument is reduced to size 2 -r .

					exp(x)		cos(x) + i sin(x) = exp(ix)
	B	n	Data	Precomp	r	Time	Precomp	r	Time
	3333	0				2.89e-05			3.56e-05
		2	0.8 KiB	5.33e-05	11	2.88e-05	6.34e-05	11	3.49e-05
		4	1.6 KiB	5.42e-05	15	2.71e-05	7.35e-05	22	2.74e-05
		8	3.3 KiB	7.61e-05	32	2.06e-05	9.65e-05	33	2.72e-05
		16	6.5 KiB	0.000131	73	1.78e-05	0.000136	37	2.89e-05
		32	13.0 KiB	0.000268	60	1.97e-05	0.000411	38	2.92e-05
		64	26.0 KiB	0.000605	60	2.2e-05	0.00104	38	3.15e-05
	10000	0				0.000202			0.000207
		2	2.4 KiB	0.000238	11	0.000183 0.000281	13	0.000209
		4	4.9 KiB	0.000240	27	0.000137 0.000333	30	0.000159
		8	9.8 KiB	0.000335	52	0.000106 0.000412	41	0.000144
		16	19.5 KiB	0.000579	83	8.48e-05	0.000633	61	0.000114
		32	39.1 KiB	0.00123	86	8.75e-05	0.00187	47	0.000129
		64	78.1 KiB	0.00270	72	9.71e-05	0.00468	47	0.000131
	33333	0				0.00166			0.00178
		2	8.1 KiB	0.00135	18	0.00135	0.0016	13	0.00167
		4	16.3 KiB	0.00136	44	0.00107	0.00186	30	0.00133
		8	32.6 KiB	0.00199	56	0.000938 0.00239	65	0.00110
		16	65.1 KiB	0.00330	89	0.000748 0.00371	90	0.000932
		32 130.2 KiB 0.00683	139 0.000637 0.0103	138 0.000841
		64 260.4 KiB 0.0152	168 0.000614 0.0256	63	0.00103
	100000	0				0.00895			0.0125
		2	24.4 KiB	0.00679	18	0.00747	0.00786	17	0.0119
		4	48.8 KiB	0.0068	44	0.00638	0.00922	40	0.00987
		8	97.7 KiB	0.00977	71	0.00565	0.0119	65	0.00754
		16 195.3 KiB 0.0164	106 0.00534	0.0179	90	0.00625
		32 390.6 KiB 0.0337	161 0.00445	0.0491	138 0.00523
		64 781.2 KiB 0.0755	240 0.00383	0.125	126 0.00612
	1000000	0				0.221			0.337
		2	244.1 KiB 0.159	18	0.195	0.187	17	0.322
		4	488.3 KiB 0.159	47	0.175	0.219	40	0.295
		8	976.6 KiB 0.228	99	0.154	0.271	96	0.273
		16	1.9 MiB	0.37	142 0.140	0.419	118 0.260
		32	3.8 MiB	0.77	161 0.136	1.14	171 0.255
		64	7.6 MiB	1.72	454 0.120	2.91	391 0.178
	10000000	0				4.36			6.50
		2	2.4 MiB	3.02	18	3.89	3.56	17	6.18
		4	4.8 MiB	3.01	47	3.53	4.1	40	5.75
		8	9.5 MiB	4.14	110 3.18	5.03	109 5.24
		16	19.1 MiB	6.57	222 2.90	7.49	203 5.13
		32	38.1 MiB	13.8	338 2.61	20.6	348 4.64
		64	76.3 MiB	31.3	551 2.39	53.4	592 4.50

The factor r 3+ε can be improved to r 2+ε using a quasilinear version of LLL[START_REF] Novocin | An LLL-reduction algorithm with quasi-linear time complexity[END_REF].

This estimate is not completely accurate because a squaring is somewhat cheaper than a multiplication (theoretically requiring 2/3 as much work). The same remark also concerns series evaluation, where some operations are squarings. We also mention that computing exp(x/2 r ) with the bit-burst algorithm might be faster than using the sinh series at this level of precision, though probably not by much; we use the sinh series here for the purposes of illustration since the analysis is simpler.

millions of digits. As shown below, it is possible to achieve larger speedup at very high precision by choosing a larger n. 4 shows sample results for the precomputation phase of Algorithm 1 to generate tables of approximate relations over n logarithms or arctangents.

Precomputation of reduction tables. Table

Here we choose the convergence factor C = 10 (each approximate relation ε i adds one decimal) and we terminate before the first relation with a coefficient |d i,j | ≥ 2 15 . This bound was chosen for convenience of storing table entries in 16-bit integers; it is also a reasonable cutoff since larger exponents will pay off only for multi-million B (as we will see below). We test the method up to n = 64, where the smallest tabulated ε i corresponds to an argument reduction of more than r = 800 bits. 9 9 Part of the implementation uses machine-precision floating-point numbers with a limited exponent range, making |ε i | < 2 -1024 ≈ 10 -300 inaccessible. Like the 16-bit limit, this is again a

Acknowledgements

The author was present at RISC in 2011 where Arnold Schönhage gave one of the talks [START_REF] Schönhage | Fast medium precision computation of elementary functions[END_REF] presenting his original "medium-precision" version of the algorithm using a pair of primes. Ironically, the author has no memory of the event beyond the published talk abstract; the inspiration for the present work came much later, with Machin-like formulas for logarithms as the starting point, and the details herein were developed independently. Nevertheless, Schönhage certainly deserves credit for the core idea. We have tried unsuccessfully to contact Schönhage (who is now retired) for notes about his version of the algorithm.

The author learned about the process to find Machin-like formulas thanks to MathOverflow comments by Douglas Zare and the user "Meij" [START_REF] Zare | MathOverflow answer: Machin-like formulas for logarithms[END_REF] explaining the method and pointing to the relevant chapter in Arndt's book.

Algorithm 1 was inspired by a comment by Simon Puchert in 2018 proposing an iterative argument reduction using smooth fractions of the form (m + 1)/m. We have substantially improved this algorithm by using LLL to look for arbitrary smooth fractions close to 1 instead of restricting to a set of fractions of special form, and by working with the logarithmic forms during reduction.

The author was supported by ANR grant ANR-20-CE48-0014-02 NuSCAP.