Fredrik Johansson
email: fredrik.johansson@gmail.com

Marc Mezzarobba

FAST AND RIGOROUS ARBITRARY-PRECISION COMPUTATION OF GAUSS-LEGENDRE QUADRATURE NODES AND WEIGHTS

Keywords: Legendre polynomials, Gauss-Legendre quadrature, arbitrary-precision arithmetic, interval arithmetic AMS subject classifications. 65Y99, 65G99, 33C45

We describe a strategy for rigorous arbitrary-precision evaluation of Legendre polynomials on the unit interval and its application in the generation of Gauss-Legendre quadrature rules. Our focus is on making the evaluation practical for a wide range of realistic parameters, corresponding to the requirements of numerical integration to an accuracy of about 100 to 100 000 bits. Our algorithm combines the summation by rectangular splitting of several types of expansions in terms of hypergeometric series with a fixed-point implementation of Bonnet's three-term recurrence relation. We then compute rigorous enclosures of the Gauss-Legendre nodes and weights using the interval Newton method. We provide rigorous error bounds for all steps of the algorithm. The approach is validated by an implementation in the Arb library, which achieves order-of-magnitude speedups over previous code for computing Gauss-Legendre rules suitable for precisions in the thousands of bits.

1. Introduction. The Legendre polynomials P n (x) are the sequence of orthogonal polynomials with respect to the unit weight on the interval (-1, 1), normalized so that P n (1) = 1. Like other classical orthogonal polynomials, Legendre polynomials satisfy a three-term recurrence, in this case the relation [START_REF] Antonov | An estimate of the remainder in the expansion of the generating function for the Legendre polynomials (generalization and improvement of Bernstein's inequality)[END_REF] (n + 1)P n+1 (x) -(2n + 1)xP n (x) + nP n-1 (x) = 0, n ≥ 1, also known as Bonnet's formula, and a second order differential equation, here

(1 -x 2 )P n (x) -2xP n (x) + n(n + 1)P n (x) = 0.

The definition implies that P n has n roots all located in (-1, 1). Perhaps the most important application of Legendre polynomials is the Gauss-Legendre quadrature rule (3)

1 -1 f (x)dx ≈ n-1 i=0 w i f (x i ), w i = 2 (1 -x 2 i )[P n (x i )] 2 ,
where the nodes x i are the roots of P n . The quantity w i is called the weight associated with the node x i . For some applications in computer algebra, number theory, mathematical physics, and experimental mathematics, it is necessary to compute integrals to an accuracy of hundreds of digits, and occasionally even tens of thousands of digits [START_REF] Bailey | High-precision numerical integration: Progress and challenges[END_REF][START_REF] Bailey | Integrals of the Ising class[END_REF][START_REF] Broadhurst | Feynman integrals, L-series and Kloosterman moments[END_REF][START_REF] Molin | Numerical Integration and L-Functions Computations[END_REF]. The Gauss-Legendre formula (3) achieves an accuracy of p bits using n = O(p) evaluation points if f is analytic on a neighborhood of (-1, 1), and if the neighborhood is large (that is, if the path of integration is well isolated from any singularities of f ), then the constants hidden in the O notation are close to the best achievable by any quadrature rule [START_REF] Kowalski | Is Gauss quadrature optimal for analytic functions?[END_REF]. This quality is related to the fact that (3) maximizes the order of accuracy among n-point quadrature rules for integrating polynomials, being exact when f is any polynomial of degree up to 2n -1; as a result, the accuracy is also excellent for analytic integrands that are well approximated by polynomials. 1In general, the error in (3) can be bounded in terms of sup x∈(-1,1) |f (2n) (x)|, or, if f is analytic on an elliptical domain D with foci at ±1, in terms of sup z∈D |f (z)| and the semi-axes of the ellipse. Even when the conditions are not ideal for using [START_REF] Bailey | Integrals of the Ising class[END_REF] directly, rapid convergence is often possible by combining [START_REF] Bailey | Integrals of the Ising class[END_REF] with adaptive subdivision of the integration path [START_REF] Petras | Self-validating integration and approximation of piecewise analytic functions[END_REF]. We give some elements of comparison between Gauss-Legendre quadrature and alternative methods, such as Clenshaw-Curtis quadrature, in section 8.

The Gauss-Legendre scheme has the drawback that the quadrature nodes and weights are somewhat inconvenient to compute. Indeed, P n becomes highly oscillatory for large n and hence presents difficulties for naive root-finding and polynomial evaluation methods. The classical Golub-Welsch algorithm avoids accuracy problems by formulating the task of computing the nodes as finding the eigenvalues of a tridiagonal matrix [START_REF] Golub | Calculation of Gauss quadrature rules[END_REF], but this approach is still too slow to be practical for large n.

In the last decade, several authors have contributed to the development of asymptotic methods that permit computing any individual node and weight for arbitrarily large n in O(1) time, culminating in the 2014 work by Bogaert [START_REF] Bogaert | O(1) computation of Legendre polynomials and Gauss-Legendre nodes and weights for parallel computing[END_REF][START_REF] Hale | Fast and accurate computation of Gauss-Legendre and Gauss-Jacobi quadrature nodes and weights[END_REF][START_REF] Bogaert | Iteration-free computation of Gauss-Legendre quadrature nodes and weights[END_REF]. For a review of this progress, see Townsend [START_REF] Townsend | The race for high order Gauss-Legendre quadrature[END_REF]. Of course, the "O(1)" bound assumes that a fixed level of precision is used. In the prevailing literature this generally means 53-bit IEEE 754 floating-point arithmetic. In addition, the available O(1) implementations rely in part on heuristic error estimates without rigorously proved bounds.

The literature on arbitrary precision or rigorous evaluation is comparatively limited. Petras [START_REF] Petras | On the computation of the Gauss-Legendre quadrature formula with a given precision[END_REF] gave explicit bounds for the error |x (i) k -x k | when the roots x k of the Legendre polynomial P n are approximated using Newton iteration [START_REF] Bailey | ARPREC: An arbitrary precision computation package[END_REF] x

(i+1) k = x (i) k - P n (x (i) k ) P n (x (i) k )
provided that the initial values x (0) k are computed by a certain asymptotic formula. However, Petras did not address the numerical evaluation of P n (x). Fousse [START_REF] Fousse | Accurate multiple-precision Gauss-Legendre quadrature[END_REF] discussed the rigorous implementation of Gauss-Legendre quadrature using generic polynomial root isolation methods together with interval Newton iteration for root refinement, but did not study fast methods for large n. Code for high-precision Gauss-Legendre quadrature rules can also be found in packages such as Pari/GP [START_REF]PARI/GP version 2.9[END_REF] and ARPREC [START_REF] Bailey | ARPREC: An arbitrary precision computation package[END_REF], but without error analysis and without special techniques for large n.

In the present article, we are interested in computing Gauss-Legendre nodes and weights to precisions p significantly larger than machine precision typically in the hundreds to thousands of bits , with rigorous error bounds. As mentioned earlier, certain applications require enclosing the values of integrals to accuracies of this order, and it is often reasonable to use quadrature rules of degree n that grows roughly linearly with p for this purpose. For example, Johansson and Blagouchine [START_REF] Johansson | Computing Stieltjes constants using complex integration[END_REF] study the computation of Stieltjes constants to precisions of hundreds of digits using complex integration, building among other things on the work described in the present paper.

If we assume that the precision p varies, basic arithmetic operations are no longer constant-time. It is well-known that addition, multiplication and division of p-bit numbers (of bounded exponent) can be performed in O(p) operations [START_REF] Brent | Modern Computer Arithmetic[END_REF], where the notation O(•) means that we neglect logarithmic factors. It is then clear that any node and weight can be computed to p-bit accuracy in O(np) time, by performing O(log p) Newton iterations (4) from an appropriate initial value. As a consequence, the full set of nodes and weights for the degree-n quadrature rule can be computed in O(n 2 p) time. For numerical integration of analytic functions where we typically have p ≈ n, a better (indeed, optimal) estimate than the classical O(n 3 ) bound is possible.

Theorem 1. If p ∼ α n for some fixed α, then the Gauss-Legendre nodes and weights of degree n can be computed to p-bit accuracy in O(n 2 ) (equivalently, O(p 2 )) bit operations.

Proof sketch. Using the formulae in [START_REF] Petras | On the computation of the Gauss-Legendre quadrature formula with a given precision[END_REF], we can compute good initial values for Newton iteration in O(n) bit operations. The Newton iteration can be performed for all roots simultaneously using fast multipoint evaluation, which costs O(np) bit operations. Fast multipoint evaluation is numerically unstable and generically loses O(n) bits of accuracy, but we can compensate for this loss by using O(n) guard bits [START_REF] Kobel | Fast approximate polynomial multipoint evaluation and applications[END_REF]. Since p ∼ αn by assumption, this does not change the complexity bound.

Completing the details of the proof is a technical exercise. Despite being elegant in theory, the algorithm behind Theorem 1 has a high overhead in practice, in large part due to the need to work with greatly increased precision. Working with expanded polynomials and processing all roots simultaneously also results in high memory usage and makes parallelization difficult. As discussed in subsection 5.2 below, an approach based on the "bit-burst" evaluation method for hypergeometric series leads to a similar complexity bound and may prove more practical for extremely large p, but likely not for p ≤ 10 6 . We can achieve a slightly worse but still subcubic complexity of O(n 5/2 ) by employing fast multipoint evaluation in a completely different way to compute P n values in isolation [START_REF] Johansson | Evaluating parametric holonomic sequences using rectangular splitting[END_REF], but unfortunately that algorithm also has high overhead.

In this work, we develop rigorous and more practical alternatives to the asymptotically fast algorithm outlined above. Our main contribution is to give a complete evaluation strategy for Legendre polynomials on [-1, 1] in ball arithmetic [START_REF] Van Der Hoeven | Ball arithmetic[END_REF][START_REF] Johansson | Arb: efficient arbitrary-precision midpoint-radius interval arithmetic[END_REF]. Computing the Gauss-Legendre nodes, then, is a relatively simple application of the results in [START_REF] Petras | On the computation of the Gauss-Legendre quadrature formula with a given precision[END_REF] together with the interval Newton method [START_REF] Moore | Methods and applications of interval analysis[END_REF]. For generating Gauss-Legendre quadrature rules with n ∼ αp, our algorithm has an asymptotic bit complexity of O(n 3 ) like classical methods, but much lower overhead. For parameters p, n ≤ 10 5 which are most relevant to applications, the observed running time is effectively subcubic. Our algorithm outperforms that of Theorem 1 for practically any combination of n and p in that range. Furthermore, if p = O(1), the complexity reduces to O(n) as in the machine-precision implementations by Bogaert [START_REF] Bogaert | Iteration-free computation of Gauss-Legendre quadrature nodes and weights[END_REF] and others.

The remainder of this paper is organized as follows. Section 2 gives an overview of our algorithm for evaluating Legendre polynomials. This is a hybrid algorithm that switches between different methods, detailed in the following sections, depending on the values of n, p, and x. In section 3, we prove practical error bounds for the three-term recurrence [START_REF] Antonov | An estimate of the remainder in the expansion of the generating function for the Legendre polynomials (generalization and improvement of Bernstein's inequality)[END_REF], which can be efficiently implemented in fixed-point arithmetic. This method is ideal for n and p up to a few hundred. For larger n or p, we use a fast method for evaluation of hypergeometric series. Section 4 discusses the hypergeometric series expansions that are preferable for different inputs and precision (including a well-known asymptotic expansion for large n), and section 5 their efficient evaluation. In section 6, we propose a strategy to select the best formula for any combination of n, p, x. Section 7 presents benchmark results that compare the performance of our algorithm to some previous implementations as well as the asymptotically fast algorithm in Theorem 1. Finally, section 8 reviews the viability of Gauss-Legendre quadrature compared to other methods for extremely high precision integration.

Our code for evaluating Legendre polynomials and computing Gauss-Legendre nodes and weights is freely available as part of the Arb library [19]2 .

2. General strategy. We work in the framework of midpoint-radius interval arithmetic, also called ball arithmetic [START_REF] Van Der Hoeven | Ball arithmetic[END_REF][START_REF] Johansson | Arb: efficient arbitrary-precision midpoint-radius interval arithmetic[END_REF]. In general, given an integer n and a ball x = [m ± r] = [m -r, m + r], we want to evaluate P n (x) at x, yielding an enclosure y = [m ± r ] such that P n (ξ) ∈ y holds for all ξ ∈ x.

We restrict our attention to real x ∈ [-1, 1], which is the most interesting part of the domain for applications. Since P n (-x) = (-1) n P n (x), we can further restrict to 0 ≤ x ≤ 1. Bogaert [START_REF] Bogaert | Iteration-free computation of Gauss-Legendre quadrature nodes and weights[END_REF] suggests working with P n (cos(θ)) instead of P n (x) directly to improve numerical stability for x close to ±1. This is not necessary in arbitraryprecision arithmetic since a slight precision increase (of the order of O(log n) bits, since the distance between two successive roots of P n close to ±1 is about 1/n 2 ) works as well.

We note that, for rigorous evaluation of P n (z) with complex z as well as Legendre functions of non-integer order n, generic methods for the hypergeometric 2 F 1 function are applicable if n is not extremely large; see [START_REF] Johansson | Computing hypergeometric functions rigorously[END_REF]. Real |x| > 1 can also be handled easily using naive methods.

In view of the use of Newton's method to compute the roots, we also need to evaluate the derivative P n (x), typically at the same time as P n (x) itself. A simple option is to deduce P n (x) from P n (x) and P n-1 (x) using ( 5)

(x 2 -1)P n (x) = n xP n (x) -P n-1 (x) .
When x is close to ±1, though, this formula involves a cancellation of about |log 2 (1-x)| bits in the subtraction, followed by a division by x 2 -1. Therefore, a direct evaluation of P n (x) may be preferable to reduce the working precision. We use either of these strategies depending on the values of n, p, and x.

Our evaluation algorithms rely on ball arithmetic internally to propagate the error bounds up to the final result. Therefore, to ensure that the enclosure output by our evaluation algorithm contains the image of the input, it is enough to have bounds on the truncation errors associated to each of the approximate expressions of P n that we use. The corresponding bounds are stated in equations ( 19), [START_REF] Olver | NIST Handbook of Mathematical Functions[END_REF], and (29).

To limit overestimation and computational overhead, we deviate from the direct use of ball arithmetic on two occasions. First, Algorithm 1 does not keep track of round-off errors internally: instead, we prove an a priori bound on the accumulated error (Corollary 6) and add it to the radius of the output ball after calling that algorithm. Second, since some methods would produce unsatisfactorily large enclosures when executed on input balls x = [m ± r] of radius r > 0, we evaluate P n (m) (with higher internal precision if necessary) and P n (x) and use a first-order bound

max ξ∈x |P n (ξ) -P n (m)| ≤ r max ξ∈x |P n (ξ)|
to separately bound the propagated error. Similarly, we use a bound for P n to compute a reasonably tight enclosure for P n ([m ± r]). Suitable bounds are given in Proposition 3. Denote by [z n ]f (z) the coefficient of index n in a power series f (z), and write f (z) z f (z) for two power series f, f if f has nonnegative coefficients and

|[z n ]f (z)| ≤ [z n ] f (z) for all n.
Lemma 2. If f , g, f , ĝ are such that f (z) z f (z) and g(z) z ĝ(z), then z 0 f z z 0 f and f (z)g(z) z f (z)ĝ(z). Proposition 3. The following bounds hold for -1 ≤ x ≤ 1:

|P n (x)| ≤ min 2 3/2 √ π √ n (1 -x 2 ) 3/4 , n(n + 1) 2 , ( 6
)
|P n (x)| ≤ min 2 5/2 √ π n 3/2 (1 -x 2 ) 5/4 , (n -1)n(n + 1)(n + 2) 8 . (7)
Proof. It is classical that Legendre polynomials are given by the generating series ( 8)

F (x, z) = ∞ n=0 P n (x)z n = 1 √ 1 -2xz + z 2 .
Differentiation with respect to x yields

∞ n=0 P n (x)z n = zF (x, z) 1 -2xz + z 2 , ∞ n=0 P n (x)z n = 3z 2 F (x, z) (1 -2xz + z 2 ) 2 .
Set θ = arccos x, so that the roots of 1 -2xz + z 2 are e ±iθ . Then, in the notation of Lemma 2, we have the bound

1 1 -2xz + z 2 = 1 2i sin θ 1 z -e iθ - 1 z -e -iθ = ∞ n=0 sin (n + 1)θ sin(θ) z n z sin(θ) -1 1 -z .
In addition, Bernstein's inequality for the Legendre polynomials [37] (see also [START_REF] Chow | A Bernstein-type inequality for the Jacobi polynomial[END_REF]) combined with the logarithmic convexity of the Gamma function yields

|P n (x)| ≤ √ 2 √ π 1 √ sin θ 1 n + 1/2 ≤ √ 2 √ π 1 √ sin θ Γ(n + 1/2) Γ(n + 1) ,
and hence

F (x, z) z √ 2 √ sin θ ∞ n=0 1 √ π Γ(n + 1/2) Γ(n + 1) z n = √ 2 √ sin θ 1 √ 1 -z .
By Lemma 2, these bounds combine into

dF dx z √ 2 sin(θ) 3/2 1 (1 -z) 3/2 , d 2 F dx 2 z √ 2 sin(θ) 5/2 3z 2 (1 -z) 5/2 . Since [z n ](1 -z) -k-1/2 = Γ(n + 1/2)/(Γ(k + 1/2)Γ(n -k + 1)
) and using again the logarithmic convexity of Γ, we conclude that

|P n (x)| ≤ √ 2 sin(θ) 3/2 2 √ π Γ(n + 1/2) Γ(n) ≤ 2 3/2 √ π √ n sin(θ) 3/2 , |P n (x)| ≤ 2 5/2 √ π n 3/2 (1 -x 2 ) 5/4 .
The result follows since all derivatives of Legendre polynomials reach their maximum at x = 1 (or by using the bounds (z-e ±iθ ) -1 , F (x, z) z (1-z) -1 and Lemma 2). Remark 4. By the same reasoning, the inequality

|P (k) n (x)| ≤ 2 k+1/2 √ π n k-1/2 (1 -x 2 ) (2n+1)/4
actually holds for all k. Unfortunately, it seems to overestimate the envelope of |P Based on these reductions, we assume from now on that x is a floating-point number with 0 ≤ x ≤ 1. Our main algorithm for evaluating P n at x combines the following methods:

• the iterative computation of P n (x) via the three-term recurrence (1),

• the asymptotic expansion (17

) of P n (x) as n → ∞,
• the usual expanded expression ( 23), ( 24) of P n in the monomial basis,

• the analogous terminating expansion ( 27) at 1. All three expansions can be written as hypergeometric series, i.e., sums of the form

k c k ξ k where c k /c k-1 is a rational function of k.
The constraints and heuristics used to select between these methods are described in detail below. Roughly speaking, the three-term recurrence is used for small index n and precision p, when x is not too close to 1; the asymptotic series when n is large enough, again with x not too close to 1; the expansion at 0 for large p unless x is close to 1; and finally the expansion at 1 in the remaining cases when x is close to 1.

For an evaluation at p-bit precision, we choose parameters such as truncation orders and internal working precision to target an absolute error of 2 -p-p for some p = O(log n), corresponding to a relative error of about 2 -p measured with respect to monotone envelopes for P n (x) and P n (x) as in [START_REF] Bogaert | O(1) computation of Legendre polynomials and Gauss-Legendre nodes and weights for parallel computing[END_REF]. The relative error of a computed ball for P n (x) where x is near a zero x k can be arbitrarily large, but the relative error of P n (x) near x k will be small, which is sufficient for the Newton iteration method. Since the output consists of a ball, we also have the option of catching a result with large relative error and repeating the evaluation with a higher precision as needed.

3. Basecase recurrence. For small n, a straightforward way to compute P n (x) is to apply the three-term recurrence (1), starting from P 0 (x) = 1 and P 1 (x) = x. Computing P n (x) by this method takes about (M(t) + O(t)) n bit operations, where t is the working precision and M(t) denotes the cost of t-bit multiplication. It is thus attractive for small n and t, especially when both P n (x) and P n (x) are needed, since we can get P n-1 (x) at no additional cost.

Fix x ∈ [-1, 1], and let p n = P n (x). Bonnet's formula (1) gives ( 9)

p n+1 = 1 n + 1 (2n + 1)xp n -np n-1 , n ≥ 0.
In a direct implementation of this recurrence in ball arithmetic, the width of the enclosures would roughly double at every iteration, requiring to increase the internal working precision by O(n) bits. We avoid this issue by performing an a priori round-off error analysis (to be presented now) of the evaluation that yields a less pessimistic bound on the accumulated error. Additionally, the static error bound allows us to implement the recurrence in fixed-point arithmetic, avoiding the overhead of floatingpoint and interval operations. Suppose x = x 2 -t with x ∈ Z is a given fixed-point number. Let u denote the integer truncation of a real number u (note that this is not the same thing as rounding to the nearest integer; however, any rounding function would do). The integer sequence (p n ) defined by

(10) p0 = 2 t , p1 = x, pn+1 = 1 n + 1 (2n + 1) xp n 2 -t -np n-1
is easy to compute using only integer arithmetic, and pn 2 -t is an approximation of p n . Algorithm 1 provides a complete C implementation using GMP [START_REF] Granlund | the GMP development team, GNU MP: The GNU Multiple Precision Arithmetic Library[END_REF]. As a small optimization, we delay the division by n + 1 until we have accumulated a denominator of the size of a machine word.

Algorithm 1 Evaluation of Legendre polynomials in GMP fixed-point arithmetic

Input: An integer x and t ≥ 0 such that |2 -t x| ≤ 1, and n ≥ 1 mpz_set(q, x); q0 = x 6:

Output: p, q such that |2 -t p -Pn-1(2 -t x)|, |2 -t q -Pn(2 -t x)| ≤ (0.75 (n + 1)(n + 2) + 1) 2 -
for (k = 1; k < n; k++) { 7:

mpz_mul(tmp, q, x); mpz_tdiv_q_2exp(tmp, tmp, t);

x q k-1 2 -t 8:
mpz_mul_si(p, p, -k*k) mpz_tdiv_q_ui(q, q, den); } else den = denlo;

d k = (k + 1) d k-1 16: } 17:
mpz_tdiv_q_ui(p, p, den/n); mpz_tdiv_q_ui(q, q, den); To bound the difference |p n 2 -t -p n |, we analyze the effect on the result of a small perturbation in each iteration of [START_REF] Broadhurst | Feynman integrals, L-series and Kloosterman moments[END_REF]. The bound is based on a classical linearity argument (compare, e.g., [START_REF] Wimp | Computation with Recurrence Relations[END_REF]) combined with majorant series techniques.

Proposition 5. Suppose that a real sequence (p n ) n≥-1 satisfies p0 = 1 and

(11) pn+1 = 1 n + 1 (2n + 1)xp n -np n-1 + ε n , n ≥ 0.
for arbitrary real numbers ε n with |ε n | ≤ ε for all n. Then we have

|p n -P n (x)| ≤ (n + 1)(n + 2) 4 ε
for all n ≥ 0.

Proof. Let δ n = pn -p n and η n = (n + 1)ε n . Subtracting ( 9) from [START_REF] Chudnovsky | Computer algebra in the service of mathematical physics and number theory[END_REF] gives

(12) (n + 1)δ n+1 = (2n + 1)xδ n -nδ n-1 + η n ,
with δ 0 = 0. Consider the formal generating series δ(z) = n≥0 δ n z n and η(z) = n≥0 η n z n . Noting that ( 12) holds for all n ∈ Z if the sequences (δ n ) and (η n ) are extended by 0 for n < 0 and using the relations

∞ n=-∞ f n-1 z n = z ∞ n=-∞ f n z n , ∞ n=-∞ nf n z n = z d dz ∞ n=-∞ f n z n ,
we see that [START_REF] Fousse | Accurate multiple-precision Gauss-Legendre quadrature[END_REF] translates into

(1 -2xz + z 2 )z d dz δ(z) = z(x -z)δ(z) + zη(z).
The solution of this differential equation with δ(0) = 0 reads, cf. [START_REF] Brent | Modern Computer Arithmetic[END_REF],

δ(z) = p(z) z 0 η(w) p(w) dw, p(z) = ∞ n=0 p n z n = F (x, z) = 1 √ 1 -2xz + z 2 .
This is an exact expression of the "global" error δ in terms of the "local" errors ε n .

Since

|p n | = |P n (x)| ≤ 1 and |η n | ≤ (n + 1)ε, it follows by Lemma 2 that |δ n | = [z n ] p(z) z 0 η(w) p(w) dw ≤ [z n ] 1 1 -z z 0 ε (1 -w) 2 1 1 -w dw
and therefore

|δ n | ≤ [z n ] 1 2 ε (1 -z) 3 = (n + 1)(n + 2) 4 ε.
Corollary 6. Suppose that x = x 2 -t for some t ≥ 0 and x ∈ Z. The sequence (p n ) n≥0 defined by (10) satisfies [START_REF] Golub | Calculation of Gauss quadrature rules[END_REF] |p n 2 -t -p n | ≤ 0.75 (n + 1)(n + 1) 2 -t , n ≥ 0.

Furthermore, the quantities p, q returned by Algorithm 1 are such that

(14) |p -2 t P n-1 (x)|, |q -2 t P n (x)| ≤ 0.75(n + 1)(n + 2) + 1.
Proof. We can write

pn+1 = 1 n + 1 (2n + 1)(x pn 2 -t + α n ) -n pn-1 + β n
for some α n , β n of absolute value at most one, and hence

pn+1 = 1 n + 1 (2n + 1) x pn 2 -t -n pn-1 + ε n , ε n = 2n + 1 n + 1 α n + β n .
where |ε n | ≤ 3. Proposition 5 applied to pn = pn 2 -t then provides the bound [START_REF] Golub | Calculation of Gauss quadrature rules[END_REF].

Turning to Algorithm 1, let p 0 , q 0 , d 0 denote the values of the variables p, q, den before the loop, and p k , q k , d k their values at the end of iteration k. Consider the sequence pk = 2 -t q k-1 /d k-1 , k ≥ 1, extended by p0 = 1 and an arbitrary (finite) p-1 . For all k ≥ 1, depending whether the conditional branch is taken, we have one of the systems of equations

p k = q k-1 , q k = (2k + 1) xq k-1 2 -t -k 2 p k-1 , d k = (k + 1) d k-1 (15) p k = q k-1 d k-1 , q k = (2k + 1) xq k-1 2 -t -k 2 p k-1 d k-1 , d k = k + 1. (16)
In both cases, we can write

p k d k = q k-1 (k + 1)d k-1 + α k k + 1 , q k d k = (2k + 1)(xq k-1 + β k ) -k 2 p k-1 (k + 1)d k-1 + γ k k + 1 with |α k |, |β k |, |γ k | ≤ 1. The first equation implies 2 -t k p k-1 /d k-1 = pk-1 + α k-1 2 -t for k ≥ 2.
Since the latter equality also holds for k = 1 with α 0 = 0, we can substitute it in the second equation, yielding

pk+1 = (2k + 1)xp k -k pk-1 k + 1 + 2 -t -kα k-1 k + 1 + (2k + 1)β k (k + 1)d k-1 + γ k k + 1 .
This relation holds for k ≥ 1, and we extend it to k = 0 by setting β 0 = γ 0 = 0. Thus, pk also satisfies [START_REF] Chudnovsky | Computer algebra in the service of mathematical physics and number theory[END_REF] with |ε n | ≤ 3 • 2 -t , and Proposition 5 applies. The final values of q and p are respectively 2 t pn and np n-1 /d n-1 = 2 t pn-1 + α n-1 , whence the bound [START_REF] Granlund | the GMP development team, GNU MP: The GNU Multiple Precision Arithmetic Library[END_REF].

We do not use asymptotically faster evaluation techniques for large n in combination with this recurrence, since the series expansions to be presented next perform very well in this case.

4. Series expansions. For large n or p, we employ series expansions of P n (x) with respect to either n or x rather than the algorithm from the previous section. The coefficients of the series are also computed by recurrence, but fewer than n terms will typically be required. Let us now review the various series expansions that we are using (an asymptotic expansion as n → ∞, series expansions at x = 0 and x = 1), before discussing their efficient evaluation in the next section. P n (cos(θ)) = 2 π sin(θ)

1/2 K-1 k=0 C n,k cos(α n,k (θ)) sin k (θ) + ξ n,K (θ)
(for arbitrary K ≥ 1), where

C n,k = [Γ(k + 1/2)] 2 Γ(n + 1) π2 k Γ(n + k + 3/2)Γ(k + 1) , α n,k (θ) = (n + k + 1/2)θ -(k + 1/2)π/2, (18)
and the error term satisfies

(19) |ξ n,K (θ)| < 2 2 π sin(θ) 1/2 C n,K sin K (θ)
.

The coefficients C n,k are a hypergeometric sequence with

(20) C n,k C n,k-1 = (2k -1) 2 4k(2n + 2k + 1) , C n,0 = 1 √ π 4 n (n + 1 2 ) 2n n .
To evaluate the error bound, we can use the following inequality deduced from ( 18):

(21) C n,k ≤ 1 πn 1/2 k!n! 2 k (n + k)! ≤ 1 πn 1/2 k! (2n) k , n, k ≥ 1.
The condition |ξ n,K (θ)| ≤ 2 -p is satisfiable as soon as sin(θ) ≥ (p + 3) ln(2)/(2n), as we can see by choosing K = 2n sin(θ) and using the inequality k! ≤ k k+1/2 e 1-k . Since the Gauss-Legendre nodes are distributed linearly in θ, the asymptotic expansion gives a candidate algorithm for all but about (log(2)/π) p/n of the nodes as p/n → 0. It is in fact a convergent series when 2 sin(θ) > 1, which allows evaluating P n (x) to unbounded accuracy for fixed n when 1 6 π < θ < 5 6 π. The particular form (17) must be used for this purpose; there is a slightly different version of the expansion which is asymptotic to P n (x) (for fixed K when n → ∞) but paradoxically converges to 2P n (x) (for fixed n when K → ∞); see [START_REF] Olver | Asymptotics and Special Functions[END_REF]Section 10.3] and [START_REF] Olver | NIST Handbook of Mathematical Functions[END_REF]Section 18.15(iii)].

We can restate [START_REF] Johansson | Efficient implementation of elementary functions in the medium-precision range[END_REF] as a hypergeometric series by working with complex numbers. Letting ω = 1 -(x/y)i, with x = cos(θ) and y = sin(θ) as usual, we have

(1 -i)(x + iy) n+1/2 ω k = √ 2 e -i π 4 e i(n+1/2)θ e ik(θ-π/2) y -k = √ 2 exp(i α n,k (θ)) sin k (θ) ,
and hence

(22) P n (x) = 1 √ πy Re (1 -i)(x + yi) n+1/2 K-1 k=0 C n,k ω k + ξ n,K (θ).
This eliminates the explicit trigonometric functions and permits using Algorithm 2 below to evaluate the series. The evaluation of ( 22) in ball arithmetic is numerically stable, and we therefore only need to add a few guard bits to the working precision.

Expansion at zero.

If n = 2d is even, the expansion of P n (x) in the monomial basis reads

P 2d (x) = (-1) d d k=0 (-1) k 2 n n d -k n + 2k n x 2k = (-1) d 2 2d 2d d d k=0 A -1 (d, k)(-x 2 ) k , (23)
and if n = 2d + 1 is odd, we have

P 2d+1 (x) = (-1) d x d k=0 (-1) k 2 n n d -k n + 2k + 1 n x 2k = (-1) d (d + 1) 2 2d+1 2d + 2 d + 1 x d k=0 A +1 (d, k)(-x 2 ) k (24)
where the hypergeometric sequences A ±1 can be defined by A ±1 (d, 0) = 1 and (25)

A σ (d, k) A σ (d, k -1) = (d -k + 1)(2d + 2k + σ) k(2k + σ) , σ ∈ {-1, +1}.
At very high precision, we evaluate the full polynomials, where [START_REF] Molin | Numerical Integration and L-Functions Computations[END_REF] and ( 24) have the advantage compared to other expansions of requiring only n/2 terms due to the odd-even form. At lower precision p, the high order terms will be smaller than 2 -p when |x| is small, and we can truncate the series accordingly and add a bound for the omitted terms to the radius of the computed ball. When the series are truncated after the k = K -1 term (for any K < d + 1), comparison with a geometric series shows that the error is bounded by the first omitted term times a simple factor. Proposition 7. For σ ∈ {-1, +1}, the error when truncating the bottom sum in [START_REF] Molin | Numerical Integration and L-Functions Computations[END_REF] or (24) (with prefactors removed) after the k = K -1 term satisfies

(26) d k=K A σ (d, k)(-x 2 ) k ≤ A σ (d, K)|x| 2K 1 -α , α = |x| 2 (d -K + 1)(2d + 2K + σ) K(2K + σ) provided that α < 1.
For bounding A σ (d, K) in this expression, and for selecting an appropriate truncation point K, we use the binomial closed forms [START_REF] Molin | Numerical Integration and L-Functions Computations[END_REF], [START_REF] Moore | Methods and applications of interval analysis[END_REF] together with the remarks in subsection 5.3.

The alternating series ( 23) and ( 24) may suffer from significant cancellation, which requires use of increased precision. We can estimate the magnitude by noting that no cancellation occurs if x is an imaginary number. Solving the majorizing recurrence

f n = 2|z|f n-1 + f n-2 with f 0 = 1, f 1 = |z| shows that |P n (z)| ≤ |P n (i|z|)| ≤ |z| + 1 + |z| 2 n .
Therefore, the possible cancellation assuming that

|P n (x)| ≈ 1 is about p A = n log 2 |x| + 1 + |x| 2 bits (which is at most n log 2 (1 + √ 2) ≈ 1.
27 n), so using ball arithmetic with about p + p A bits of working precision for the series evaluation gives p-bit accuracy.

c n,k c n,k-1 = (n -k + 1)(n + k) k 2 .
As in the previous section, we can truncate [START_REF] Petras | On the computation of the Gauss-Legendre quadrature formula with a given precision[END_REF] and bound the error by comparison with a geometric series.

Proposition 8. The error when truncating [START_REF] Petras | On the computation of the Gauss-Legendre quadrature formula with a given precision[END_REF] after the k = K -1 term satisfies

(29) n k=K c n,k u k ≤ c n,K |u| K 1 -α , α = |u| (n -K)(n + K + 1) (K + 1) 2 provided that α < 1.
For u ≥ 0, [START_REF] Petras | On the computation of the Gauss-Legendre quadrature formula with a given precision[END_REF] does not suffer from cancellation. For u < 0, we can estimate the amount of cancellation from the magnitude of P n (x ) where |u| = (x -1)/2. For not too large x ≥ 1, a very good approximation is

P n (x ) ≤ 2 ∞ k=0 n 2k k! 2 |u| k = 2 I 0 (2n |u|) ≤ 2 e 2n √ |u| .
We therefore need about 2n max(0, -u)/ ln(2) bits of increased precision. We can compute P n from P n and P n-1 , but since this involves a division by 1 -x 2 , it is better to evaluate P n directly when x is close to 1. We have

P n (x) = n-1 k=0 c n,k u k where c n,k = (k + 1)c n,k+1 /2 satisfies (30) c n,0 = n(n + 1) 2 , c n,k c n,k-1 = (n -k)(n + k + 1) k(k + 1) . Since c n,k ≤ nc n,k+1 , the analog (31)
n k=K c n,k u k ≤ n n K + 1 n + K + 1 K + 1 |u| K 1 1 -α
of Proposition 8 holds with u and α as above.

5. Fast evaluation of series expansions. All these series expansions are amenable to fast evaluation techniques specific to multiple-precision arithmetic. Using such fast summation algorithms is critical for achieving good performance at high precision. We now discuss the algorithm that we use for evaluating the series of the previous section.

Rectangular splitting.

We use rectangular splitting [START_REF] Smith | Efficient multiple-precision evaluation of elementary functions[END_REF][START_REF] Johansson | Evaluating parametric holonomic sequences using rectangular splitting[END_REF] to evaluate hypergeometric series with rational parameters where the argument x is a high-precision number. This reduces evaluating a K-term series to O(K) cheap scalar operations (additions and multiplications or divisions by small integer coefficients) and about 2 √ K expensive nonscalar operations (general multiplications), whereas direct evaluation of the hypergeometric recurrence uses O(K) expensive operations.

Algorithm 2 presents our version of rectangular splitting for the present application. We implement the various series expansions by defining the functions p(k), q(k) used in steps 6 and 17 according to formulae (20), ( 25), ( 28), [START_REF] Takahasi | Double exponential formulas for numerical integration[END_REF].

This algorithm is a generalization of the method for evaluating Taylor series of elementary functions given in [START_REF] Johansson | Efficient implementation of elementary functions in the medium-precision range[END_REF], which combines rectangular splitting with partially unrolling the recurrence to reduce the number of scalar divisions (which in practice are more costly than scalar multiplications). The terms are computed in the reverse direction to allow using Horner's rule for the outer multiplications.

Our code uses ball arithmetic for x and s so that no error analysis is needed, and we use a bignum type for c (so no overflow can occur regardless of u). For low precision a faster implementation would be possible using fixed-point arithmetic with tight control of the word-level operations as was done for elementary functions in [START_REF] Johansson | Efficient implementation of elementary functions in the medium-precision range[END_REF].

The algorithm contains two tuning parameters 3 . The splitting parameter m controls the number m of multiplications for powers versus the number K/m of multiplications for Horner's rule. The choice m ≈ √ K is optimal, but when evaluating two series for the same x (in our case, to compute both P n (x) and P n (x)), the table of powers can be reused, and then m ≈ √ 2K minimizes the total cost. The unrolling parameter u controls the number of coefficients to collect on a single denominator, reducing the number of divisions to N/u. Ideally, u should be chosen so that b j=a p(j) and b j=a q(j) fit in 1 or 2 machine words. The example value u = 4 is a reasonable default, but as an optimization, one might vary u for each iteration of the main loop to ensure that c always fits in a specific number of words.

Algorithm 2 Evaluation of hypergeometric series using rectangular splitting

Input: An arbitrary x, recurrence data p, q ∈ Z[k], integer K ≥ 0, offset Ω ∈ {0, 1}. Output: s = K-1 k=Ω x k k j=Ω p(j)/q(j)

1: m ← √ K ; precompute [1, x, x 2 , . . . , x m ]
Tuning parameter: any m ≥ 1 can be used

2: s ← 0; k ← K -1 3: while k ≥ Ω do 4: u ← min(4, k + 1 -Ω)
Tuning parameter: any 1 ≤ u ≤ k + 1 -Ω can be used c ← (c/p(k))q(k) Exact small integer division

18: k ← k -1 19:
end while 20:

s ← s/c 21: end while 22: return s

The redundant parameter Ω is a small convenience in the pseudocode. Setting Ω = 1 and adding the constant term separately avoids having to make a special case to prevent division by zero when q(0) = 0.

Due to the scalar operations, rectangular splitting ultimately requires O(K) arithmetic operations with O(p) bit complexity each just like straightforward evaluation of the recurrence, so it is not a genuine asymptotic improvement, but it is an improvement in practice and can give more than a factor 100 speedup at very high precision. It is possible to genuinely reduce the complexity of evaluating a hypergeometric sequence to O( √ K log(K)) arithmetic operations using a baby-step giant-step method that employs fast multipoint evaluation, but in practice rectangular splitting performs better until both K and p exceed 10 6 (see [START_REF] Johansson | Evaluating parametric holonomic sequences using rectangular splitting[END_REF]).

A note on the bit-burst method.

Another technique for fast evaluation of hypergeometric series, called binary splitting, would be useful when p is large and the argument x is a rational number with small numerator and denominator, but this case is not relevant for our application. Binary splitting also forms the basis of the bit-burst method [START_REF] Chudnovsky | Computer algebra in the service of mathematical physics and number theory[END_REF]Section 4], which permits evaluating any fixed hypergeometric series at any fixed point without the restriction to simple rational numbers of plain binary splitting to absolute precision p in only O(p) bit operations. Yet, we do not use this method either in our implementation, due to its large overhead.

Indeed, computing P n (x) by the bit-burst method method requires O(log p) analytic continuation steps, each of which entails two evaluations of general solutions of the Legendre differential equation [START_REF] Bailey | High-precision numerical integration: Progress and challenges[END_REF]. These solutions are defined by unit initial values at some intermediate point x 0 , and can be represented as a power series of radius of convergence 1 -|x| whose coefficients obey recurrences of order two. This is to be compared with a single series, given by a first-order recurrence and typically converging faster, for the expansions considered in subsection 4.1 to subsection 4.3. Thus, the bit-burst method is unlikely to be competitive in the range of precision we are interested in, especially when the asymptotic series can be used.

Nevertheless, it can be shown that the solutions with unit initial values at x 0 of (2) occuring at intermediate analytic continuation steps have Taylor coefficients c k at x 0 bounded by (n 2 /(1 -|x 0 |)) O(k) uniformly in n and x 0 . As a consequence, the asymptotic cost of computing any individual root of P n (x) by the bit-burst method and Newton iteration is O(p) log(n) O (1) . For computing all the roots, this approach matches the O(np) estimate of Theorem 1, while allowing for parallelization and requiring less memory. It may hence provide an alternative to multipoint evaluation worth investigating for precisions in the millions of bits.

Binomial coefficients.

The prefactors of both the series expansion at x = 0 and the asymptotic series contain the central binomial coefficient 2n n . We need to compute this factor efficiently for any n and precision p. Since 2n n ≈ 4 n , it is best to use an exact algorithm when n < Cp for some small constant C > 1/2. We use the binomial function provided by GMP for n < 6p + 200 and otherwise use an asymptotic series for 2n n with error bounds given in [7, Corollaries 1 and 2]. We also need to quickly estimate the magnitude of binomial coefficients for error bounds of series truncations. We have the binary entropy estimate

log 2 n k ≤ nG(k/n), G(x) = -x log 2 (x) -(1 -x) log 2 (1 -x)
and the equivalent form

(32) n k ≤ n n -k n-k n k k = n n k k (n -k) n-k .
The function G(x) can be evaluated cheaply with a precomputed lookup table. A coarse estimate is sufficient, since overestimating log 2 n k by a few percent only adds a few percent to the running time.

6. Algorithm selection. We first use a set of cutoffs found experimentally to decide whether to use the basecase recurrence or one of the series expansions. The recurrence is mainly faster for some combinations of p < 1 000, n < 400 when computing (P n (x), P n (x)) simultaneously and for some combinations of p < 500, n < 100 when computing P n (x) alone (in all cases subject to some boundaries ε < x < 1 -ε); the actual optimal regions are complicated due to differences in overhead between fixed-point integer arithmetic and ball arithmetic for the respective algorithm implementations. For the actual cutoffs used, we refer to the source code.

To select between the series expansion at x = 0, the expansion at x = 1, and the asymptotic series, the following heuristic is used. For each algorithm A, we estimate the evaluation cost as C A = K A (p + p A ) where K A is the number of terms required by algorithm A (K A = ∞ if A is the asymptotic series and it does not converge to the required accuracy), p is the precision goal, and p A is the extra precision required by algorithm A due to internal cancellation. For the asymptotic series, we multiply the cost by an extra factor 2 as a penalty for using complex numbers. In the end, we select the algorithm with the lowest C A .

We select K A and estimate p A heuristically using machine precision floating-point computations, working with logarithmic magnitudes to avoid underflow and overflow. 1. Time to evaluate (Pn(x), P n (x)) as the argument x = cos(θ) varies (x = 1 at θ = 0 and x = 0 at θ = π/2), here with n = 10 000, for precision p somewhat smaller than n (top plot) and somewhat larger (bottom plot). The variable θ is used for the horizontal scale in this picture to follow the distribution of the roots (which are clustered near x = 1) linearly.

For example, when A is the asymptotic series, we search for a small K A such that log(K A !/(2n sin θ) K A ) ≤ -(p + p A ) for some small p A , in accordance with [START_REF] Kobel | Fast approximate polynomial multipoint evaluation and applications[END_REF]. During the actual evaluation of the series expansions, K A and p A are then given; we compute rigorous upper bounds for the truncation error via ( 19), ( 21), ( 26), ( 29), ( 31), (32) using floating-point arithmetic with directed rounding, while additional rounding errors are tracked by the ball arithmetic.

The assumption that the running time is a bilinear function of K A and p + p A is not completely realistic, but this cost estimate nonetheless captures the correct asymptotics when

x → 0, x → 1, n → ∞, p → ∞
separately, and hopefully will not be too inaccurate in the transition regions. This is verified empirically. Figure 1 illustrates how the time to evaluate (P n (x), P n (x)) varies with x when the automatic algorithm selection is used. Here, we have timed the case n = 10 000 for two different p. For large n and p n (top plot), a sharp peak appears at the transition between the series expansion at x = 1 and the asymptotic expansion which is used for most x. This peak tends to become taller but narrower for larger n. We could presumably get rid of the peak by implementing another algorithm specifically for the transition region, but the area under the peak is so small compared to the median baseline that computing all the roots would not be sped up much. For p somewhat larger than n (bottom plot), we observe a smooth transition between the series at x = 1 near the left of the picture and the series at x = 0 used over the most of the range.

Benchmarks.

As already mentioned, our code for computing values and roots of Legendre polynomials is part of the Arb library, available from http://arblib.org/. The benchmark results given in this section were obtained using prerelease builds of version 2.13 of Arb. The implementation of the algorithms of this article is located in the files arb_hypgeom/legendre_p_ui_* of the Arb source tree. In subsection 7.1, we also use some ad hoc code available from this paper's public git repository4 when comparing the main algorithm with variants absent from the Arb implementation. The source code for the experiments themselves can be found in the same git repository. Except where otherwise noted, we ran the programs under 64-bit Linux on a laptop with a 1.90 GHz Intel Core i5-4300U CPU using a single core.

7.1. Polynomial evaluation. Figure 2 compares the performance of different methods for evaluating (P n (x), P n (x)) on a set of n/2 points distributed like the positive roots of P n (x) to simulate one stage of Newton iteration at p-bit precision. The time for the three-term recurrence is set to 1, i.e., we divide the other timings by this measurement. The following methods are timed:

• Our hybrid method (from here on called the "main algorithm") with automatic selection between the three-term recurrence and different series expansions. • The main algorithm without the three-term recurrence as the basecase, i.e., using series expansions even for very small n. • Fast multipoint evaluation of the expanded polynomials P n (x) and P n (x). As in subsection 4.2, we expand P n ( √ x) for even n and P n ( √ x)/ √ x for odd n and evaluate at x 2 since this halves the amount of work. The polynomial coefficients are generated using the hypergeometric recurrence and the fast multipoint evaluation is done using _arb_poly_evaluate_vec_fast_precomp (where the "precomp" suffix indicates that the same product tree is used for both P n and P n ). The fast multipoint evaluation is done with 2.9n guard bits, which was found experimentally to be sufficient for full accuracy. The crossover point between the three-term recurrence and series expansions usually occurs around n ≈ 10 2 -10 3 (it can be as low as n ≈ 10 if p is much larger). For modest n, the three-term recurrence is much faster than the hypergeometric series (typically by a factor 3-4) due to working with negligible extra precision and thanks to the low overhead of fixed-point arithmetic. This low overhead is very useful for typical evaluation of Legendre polynomials and generation of quadrature nodes for one or a few machine words of precision. The crossover point could be lowered slightly if we used a similarly optimized fixed-point implementation for the hypergeometric series. When p is fixed (top left in Figure 2), the main algorithm is a factor O(n) faster than the three-term recurrence since the asymptotic expansion converges to sufficient accuracy after O(1) terms for all sufficiently large n. With the constant precision p = 64, the main algorithm is 3.0 times faster for n = 10 3 and 30 times faster for n = 10 4 . Conversely, fast multipoint evaluation with constant p is a factor O(n) slower than our algorithm due to the higher internal precision.

When p ∝ n (the three remaining plots in Figure 2), the main algorithm appears to show the same O(n) speedup over the three-term recurrence after the crossover point, at least initially. This speedup should level off asymptotically, but in practice this only occurs for n larger than 10 4 where we have already gained a factor 10 or more. The leveling off is visible in the bottom right figure (p = 10n).

Fast multipoint evaluation gives a true asymptotic O(n) speedup, but since it has much higher overhead, it only starts to give an improvement over the main algorithm from n ≈ 10 4 and for p larger than n. When p = n/10, it appears that fast multipoint evaluation will only break even for n much larger than 10 5 . We conclude that fast multipoint evaluation would be worthwhile only for the last few Newton iterations when computing quadrature nodes for exceptionally high precision. Since independent evaluations are more convenient and easy to parallelize, the fast multipoint evaluation method currently seems to have limited practical value for this application. 7.2. Quadrature nodes. The function arb_hypgeom_legendre_p_ui_root(x, w, n, k, p) sets the output variable x to a ball containing the root of P n with index k (we use the indexing 0 ≤ k < n, with k = 0 giving the root closest to 1), computed to p-bit precision. It also sets w to the corresponding quadrature weight. We use the formulae in [27, Theorem 1(c)] to compute an initial enclosure with roughly machine precision, followed by refinements with the interval Newton method at doubling precision steps for very high precision. We deduce the quadrature weights thanks to the classical expression as functions of the nodes recalled in equation [START_REF] Bailey | Integrals of the Ising class[END_REF].

Table 1 shows timings for computing degree-n Gauss-Legendre rules to p-bit n is needed if the path of integration is close to a singularity or if the integrand is highly oscillatory. In such cases, bisecting the interval a few times to reduce the necessary n is often a better tradeoff. On the other hand, since the time to generate nodes with our code only grows linearly with n beyond n ≈ p, increasing the degree further is viable, and potentially useful if the integrand is expensive to evaluate.

In the present work, we refrain from a more detailed discussion of adaptive integration strategies and the computation of error bounds for the integral itself. However, we mention that Arb contains an implementation of a version of the Petras algorithm [START_REF] Petras | Self-validating integration and approximation of piecewise analytic functions[END_REF] for rigorous integration. This code uses both adaptive path subdivision and Gauss-Legendre quadrature with an adaptive choice of n up to n ≈ 0.5p by default, with degree increments n ≈ 2 k/2 and automatic caching of the nodes for fast repeated integrations. Node generation takes at most a few seconds for a first integration at 1 000-digit precision and a few milliseconds for 100-digit precision.

8. Gauss-Legendre versus Clenshaw-Curtis and the double exponential method. The Clenshaw-Curtis and double exponential (tanh-sinh) quadrature schemes have received much attention as alternatives to Gauss-Legendre quadrature for numerical integration with very high precision [START_REF] Takahasi | Double exponential formulas for numerical integration[END_REF][START_REF] Bailey | High-precision numerical integration: Progress and challenges[END_REF][START_REF] Trefethen | Is Gauss quadrature better than Clenshaw-Curtis?[END_REF]. Both schemes typically require a constant factor more evaluation points than Gauss-Legendre rules for equivalent accuracy, but the nodes and weights are easier to compute. Gauss-Legendre quadrature is therefore the best choice when the integrand is expensive to evaluate or when nodes can be precomputed for several integrations. It is of some interest to compare the relative costs empirically.

Here we assume an analytic integrand with singularities well isolated from the finite path of integration so that Gauss-Legendre quadrature is a good choice to begin with. As observed in [START_REF] Trefethen | Is Gauss quadrature better than Clenshaw-Curtis?[END_REF], Clenshaw-Curtis often converges with identical rate to Gauss-Legendre for less well-behaved integrands, and the double exponential method is far superior to either Clenshaw-Curtis or Gauss-Legendre for analytic integrands with endpoint singularities.

Clenshaw-Curtis quadrature uses the Chebyshev nodes cos(πk/n), 0 ≤ k ≤ n, and the corresponding weights can be expressed by a discrete cosine transform which takes O(n log n) arithmetic operations to compute by an FFT. As a rule of thumb, 2npoint Clenshaw-Curtis quadrature gives the same accuracy as n-point Gauss-Legendre quadrature (for instance, the 384-point Clenshaw-Curtis rule gives errors of 10 -229 , 10 -294 and 10 -154 for the three integrals in Table 3). As a point of comparison with Tables 1 and3, Arb computes a length-2 048 FFT with 1 000-digit precision in 0.09 seconds and a length-32 768 FFT with 10 000-digit precision in 36 seconds. The precomputation for Clenshaw-Curtis quadrature is therefore roughly a factor 20 cheaper than for Gauss-Legendre quadrature with our algorithm, while subsequent integration with cached weights is twice as expensive for Clenshaw-Curtis.

Double exponential quadrature uses the change of variables x = tanh( 1 2 π sinh t) to convert an integral on (-1, 1) to the interval (-∞, +∞) in such a way that the trapezoidal rule ∞ -∞ f (t)dt ≈ h n k=-n f (hk) converges exponentially fast. One generally chooses the discretization parameter as h = 2 -j so that both the evaluation points and weights can be recycled for successive levels j = 1, 2, 3 . . ., and n is chosen so that the tail of the infinite series is smaller than 2 -p . The 2n + 1 nodes and weights can be computed with n + O(1) exponential function evaluations and O(n) arithmetic operations. Double exponential quadrature with Cn evaluation points typically achieves the same accuracy as n-point Gauss-Legendre quadrature, where C is slightly larger than for Clenshaw-Curtis, e.g. C ≈ 5; see Table 4. The time to compute nodes and weights is comparable to Clenshaw-Curtis quadrature (around 0.2 seconds for 1 000-digit precision and two minutes for 10 000-digit precision).

In summary, for integration with precision in the neighborhood of 10 3 to 10 4 digits, computing n nodes with our algorithm is about an order of magnitude more expensive than performing n elementary function (e.g. exp or log) evaluations or computing an FFT of length n. This makes Gauss-Legendre quadrature competitive for computing more than m integrals (assuming that nodes and weights are cached), for a single integral requiring splitting into m subintervals, or for a single integral when the integrand costs more than m elementary function evaluations, where m ≈ 10 1 .

The picture becomes more complicated when accounting for the method used to estimate errors in an adaptive integration algorithm. One drawback of the Gauss-Legendre scheme is that the nodes are not nested, so that an adaptive strategy that repeatedly doubles the quadrature degree requires twice as many function evaluations as the degree of the final level. This drawback disappears if the error is estimated by extrapolation or if an error bound is computed a priori as in the Petras algorithm [START_REF] Petras | Self-validating integration and approximation of piecewise analytic functions[END_REF]. 9. Conclusion. In [START_REF] Bailey | High-precision numerical integration: Progress and challenges[END_REF], it was claimed that "There is no known scheme for generating Gaussian abscissa-weight pairs that avoids [the] quadratic dependence on n. High-precision abscissas and weights, once computed, may be stored for future use. But for truly extreme-precision calculations -i.e., several thousand digits or morethe cost of computing them even once becomes prohibitive".

In this quote, "quadratic dependence" refers to the number of arithmetic operations. We may remark that using asymptotic expansions in the evaluation of Legendre polynomials avoids the quadratic dependence on n for fixed precision p, and Theorem 1 avoids the implied cubic time dependence on n when n = O(p). In fact, Theorem 1 implies that Gauss-Legendre, Clenshaw-Curtis and double exponential quadrature have the same quasi-optimal asymptotic bit complexity, up to logarithmic factors.

Our experiments show that the algorithm in Theorem 1 hardly is worthwhile. However, the hybrid method described in sections 2 to 5 does achieve a significant speedup for practical p and n which allows us to compute Gauss-Legendre quadrature rules for 1 000-digit integration in 1-2 seconds and for 10 000-digit integration in 10-20 minutes on a single core. This is not prohibitively expensive compared to the repeated evaluation of typical integrands, especially if several integrations are needed. Parallelization is also trivial since all roots are computed independently.

A natural extension of this work would be to consider Gaussian quadrature rules for different weight functions. The techniques should transfer to other classical orthogonal polynomials (Jacobi, Hermite, Laguerre, etc.) which likewise have hypergeometric expansions and satisfy three-term recurrence relations. The main obstacle might be to obtain large-n asymptotic expansions with suitable error bounds.

  about 2 k in the region where it is smaller than P (k) n (1).

14 :

 14 den = k+1; d k = k + 1 15:

4. 1 .

 1 Asymptotic series. For fixed |x| < 1 or equivalently x = cos(θ) with 0 < θ < π, an asymptotic expansion for P n (x) as n → ∞ can be given as [6, Eq. 3.4][START_REF] Johansson | Efficient implementation of elementary functions in the medium-precision range[END_REF]

4. 3 .

 3 Expansion at one. Expanding at x = 1 yields (27) P n (x) = n k=0 c n,k u k , c n,k = n k n + k kwhere u = (x-1)/2. The coefficients c n,k are hypergeometric with initial value c n,0 = 1 and term ratio[START_REF] Petras | Self-validating integration and approximation of piecewise analytic functions[END_REF]

5 :while k ≥ a do 8 :r ← k mod m 9 :s ← s + c • x r Using precomputed power of x 13 : end if 14 :

 5891314 (a, b) ← (k -u + 1, k) if k = b then 10: s ← c • (s + x r )Using precomputed power of x if r = 0 and k = 0 then 15:s ← s • x m Using precomputed power of x

Fig.

  Fig.1. Time to evaluate (Pn(x), P n (x)) as the argument x = cos(θ) varies (x = 1 at θ = 0 and x = 0 at θ = π/2), here with n = 10 000, for precision p somewhat smaller than n (top plot) and somewhat larger (bottom plot). The variable θ is used for the horizontal scale in this picture to follow the distribution of the roots (which are clustered near x = 1) linearly.

Fig. 2 .

 2 Fig.2. Performance comparison of various methods to evaluate (Pn(x), P n (x)) to p-bit precision for a set of n/2 points 0 < x < 1 distributed like the roots of Pn. The y axis (relative time) shows the time divided by the time using the three-term recurrence in fixed-point arithmetic.

Table 1

 1 Time in seconds to compute the degree-n Gauss-Legendre rules with p-bit precision using our code. (For large n and p, the time was estimated by computing a subset of the nodes and weights.)

	n \ p	64	256	1 024	3 333	33 333
	20 0.000133 0.000229 0.000510 0.00121 0.0198
	50 0.000450 0.000870 0.00212 0.00520 0.0710
	100 0.00138	0.00310	0.00720	0.0163	0.191
	200 0.00550	0.0111	0.0267	0.0550	0.589
	500	0.0236	0.0610	0.164	0.325	2.61
	1 000	0.0530	0.145	0.584	1.238	9.21
	2 000	0.0860	0.298	1.12	4.20	32.6
	5 000	0.191	0.665	2.67	14.3	181
	10 000	0.350	1.26	4.93	26.6	674
	100 000	3.60	12.2	41.3	212	13 637
	1 000 000	58.0	146	411	1 850	103 960

Table 3

 3 Left columns: time in seconds to generate 1 000-digit quadrature rules for the degrees n used by ARPREC. Right columns: for three different integrals, the error | 1 -1 f (x)dx -n-1 k=0 w k f (x k )| of the degree-n quadrature rule, and the time to evaluate this degree-n approximation of the integral at 1 000-digit precision in Arb given the nodes and weights (x k , w k ).

	n ARPREC Our code	1 -1 log(2+x)dx	1 -1 Ai(10x)dx	1 -1 Γ(1+ix)dx
				Error	Time	Error Time Error Time
	12	0.00520	0.000592	10 -14		10 -1		10 -8
	24	0.0189	0.00171	10 -28		10 -9		10 -17
	48	0.0629	0.00507	10 -56		10 -34		10 -36
	96	0.251	0.0163	10 -111		10 -105		10 -73
	192	0.974	0.0532	10 -222		10 -284 0.075 10 -146
	384	3.83	0.195	10 -441 0.023 10 -721 0.15 10 -293	1.3
	768	15.2	0.763	10 -881 0.045	< ε	0.29 10 -588	2.5
	1 536	60.9	2.82	< ε	0.091			< ε	5.0
	3 072	241	9.55
	6 144	1 013	18.3

Table 4

 4 Left columns: step sizes h, number of evaluation points, and time to compute nodes for double exponential quadrature with Arb at 1 000-digit precision. Right columns: error and evaluation time given precomputed nodes.

	h	2n + 1 Time	1 -1 log(2+x)dx	1 -1 Ai(10x)dx	1 -1 Γ(1+ix)dx
				Error	Time	Error Time Error Time
	2 -7	1 989	0.07 10 -407	0.12	10 -423 0.93 10 -314	6.3
	2 -8	3 977	0.14 10 -814	0.25	10 -909 1.75 10 -630	13.0
	2 -9	7 955	0.27	< ε	0.55	< ε	3.49	< ε	25.1

However, statements about the near-optimality of Gauss-Legendre quadrature must not be over-interpreted. Indeed, the rate of convergence of Gauss-Legendre quadrature is not optimal asymptotically when n → ∞ for analytic f on a fixed neighborhood, being improvable by a small factor[START_REF] Trefethen | Six myths of polynomial interpolation and quadrature[END_REF].

http://arblib.org/

Let us stress that the choice of these parameters only affects the performance of the algorithm, not its correctness. The same remark holds every other time we resort to heuristics in this article.

https://github.com/fredrik-johansson/legendrepaper/

Acknowledgements. We are indebted to Nick Trefethen and two anonymous referees for their useful comments. We also thank Bruno Salvy for pointing out Petras' paper [27]. Marc Mezzarobba was supported in part by ANR grant ANR-14-CE25-0018-01 (FastRelax).

Table 2

Time in seconds for Pari/GP to compute the degree-n Gauss-Legendre quadrature rules with p-bit precision. The numbers in parentheses show the speedup of our code compared to Pari/GP. n \ p 64 256 1 024 3 333 33 333 20 0.00059 (×4.4) 0.00070 (×3.1) 0.0015 (×2.9) 0.0035 (×2.9) 0.078 (×3.9) 50 0.0043 (×9.6) 0.0050 (×5.8) 0.010 (×4.9) 0.022 (×4. precision by calling this function repeatedly with 0 ≤ k < n/2. Table 2 compares our code to the intnumgaussinit function in Pari/GP which uses a generic polynomial root isolation strategy followed by Newton iteration for high precision refinement. The improvement is most dramatic for small p and large n where we benefit from using asymptotic expansions, but we also obtain a consistent speedup for large p.

For low precision and large n, our implementation is about three orders of magnitude slower than the machine precision code by Bogaert [START_REF] Bogaert | Iteration-free computation of Gauss-Legendre quadrature nodes and weights[END_REF] which is reported to compute the nodes and weights for n = 10 6 in 0.02 seconds on four cores. This difference is reasonable since we use arbitrary-precision arithmetic, compute rigorous error bounds, and evaluate the Legendre polynomials explicitly whereas Bogaert uses a more sophisticated asymptotic development for both the nodes and the weights.

We also note that we can compute 53-bit floating-point values with provably correct rounding in about the same time as the 64-bit values, using Ziv's strategy of increasing the precision. For a ball with relative radius just larger than 2 -64 , there is less than a 1% probability that the correct 53-bit rounding cannot be determined, in which case that particular node can be recomputed with a few more bits.

Fousse [START_REF] Fousse | Accurate multiple-precision Gauss-Legendre quadrature[END_REF] reports a few timings for smaller n and high precision obtained on a 2.40 GHz AMD Opteron 250 CPU. For example, n = 80, p = 500 takes 0.14 seconds (our implementation takes 0.029 seconds) and n = 556, p = 5 000 takes 17 seconds (our implementation takes 0.53 seconds). Of course, these timings are not directly comparable since different CPUs were used.

The mathinit program included with version 2.2.19 of D. H. Bailey's ARPREC library generates Gauss-Legendre quadrature nodes using Newton iteration together with the three-term recurrence for evaluating Legendre polynomials [START_REF] Bailey | ARPREC: An arbitrary precision computation package[END_REF][START_REF] Bailey | High-precision numerical integration: Progress and challenges[END_REF]. With default parameters, this program computes the rules of degree n = 3 • 2 i+1 for 1 ≤ i ≤ 10 at 3 408 bits of precision, intended as a precomputation for performing degree-adaptive numerical integrations with up to 1 000 decimal digit accuracy. This takes about 1 300 seconds in total (our implementation takes 32 seconds). A breakdown for each degree level is shown in Table 3.

Table 3 also shows the approximation error and the evaluation time (not counting the computation of the nodes and weights) for the degree-n approximations of three different integrals, illustrating the relative costs and realistic requirements for n. As motivation for the third integral, we might think of a segment of a Mellin-Barnes integral. The log, Airy and gamma function implementations in Arb are used.

The last few degree levels (with n roughly larger than the number of decimal digits) used by ARPREC tend to be dispensable for well-behaved integrands. A larger


O A González-Estrada

J F Cáceres-Ramírez

W Pinto-Hernández

Virtual tool for training and evaluation of problem-based engineering subjects

In Educational Technology, one of the challenges in the use of virtual tools for engineering education is the assessment of problem-based subjects. In this work, a software generator of exercises is developed, to be contained in a virtual platform, which seeks to improve the learning process and the evaluation of the subjects based on problems such as, e.g., Mechanics of Materials. The student can use it for training, with content corresponding to the subject of the course, monitoring their performance and obtaining effective feedback. Additionally, the teacher has the possibility to periodically conduct formative assessments that allow him to recognize the concepts in which the students are failing. This tool can provide an objective instrument to normalize grading across groups directed by different lecturers. The software was written in the programming languages Java, JavaScript and HTML for the presentation of the contents, and can be used from any electronic device. This tool can be used to enhance classroom learning in face to face, blended or online learning.

Introduction

Since its inception, the World Wide Web has undergone several transformations thanks to technological advances and the needs that different types of users demand. Web 2.0 is the conception of the network as a platform, it gives people an active role, and encourage participation and usability [START_REF] Dinucci | Fragmented Future[END_REF]. These features have promoted the boom of Educational Technology aimed to enhance the teacherstudent experience, improve knowledge management [START_REF] Lache | Las tecnologías de información y comunicación como prácticas de referencia en la gestión de conocimiento : una revisión sistemática de la literatura[END_REF], and facilitate the learning process [START_REF] Moreno-Chaustre | Modelo de evaluación para valorar el cambio en las prácticas docentes con TIC[END_REF]. Within Educational Technology, electronic learning (e-learning) consists in the use of new multimedia and Internet technologies to improve the quality of learning by facilitating access to resources and services, as well as to the collaboration and remote exchange. In e-learning, the formative process takes place through a kind of virtual environment in which the teacher-student interaction takes place as well as the activities of the students with the learning materials. The e-learning serves as a means for face-to-face, online or blended teaching.

Educational Technology has many advantages, it can be used to simplify and reduce the cost for labs used to demonstrate concepts [START_REF] Machado-Mercado | Una herramienta computacional didáctica para el análisis cinemático de mecanismos planos de cuatro barras[END_REF], massify access to knowledge with initiatives as massive online open courses (MOOC), improve feedback for effective formative assessment [START_REF] Gikandi | Online formative assessment in higher education: A review of the literature[END_REF], or even help vulnerable populations by facilitating access to remote places affected by security issues [START_REF] González-Zabala | Propuesta de un sistema de información para apoyar la operación de modelos de educación en situación de emergencia : caso círculos de aprendizaje[END_REF]. Challenges for developing countries include improving Information and Communication Technologies (ICT) infrastructure, facilitating access and use of the technological resources, and encouraging the adoption of new training methods [START_REF] Douglas | Designing an ELearning Portal for Developing Countries : An Action Design Approach[END_REF].

We can find numerous efforts that have been reported in the literature to support the teachinglearning process of various subjects. In [START_REF] Torres | Software para la Enseñanza Didáctica de Mecánica de Materiales[END_REF] a program created in Macromedia Flash v8.0 focuses on the understanding of the laboratory tests that accompany the theory of the course Mechanics of Materials. It includes tests of traction, compression, shear, Poisson's effect, torsion, extensometry, bending, columns and photo-elasticity. The student has different types of evaluations that not only allow him to evaluate his learning but also prepare him for the tests that are normally carried out in theory courses.

In [START_REF] Guerra | Software educativo para el trabajo con matrices[END_REF], a didactic tool is developed with the objective of supporting the teaching-learning process of Linear Algebra in engineering degrees. When entering the problem data in the software, the operation is processed and the procedure is visualized, showing step by step how the problem is solved. Moreover, it contains theoretical material to review the concepts that the student needs to develop the exercises. The proposed system was implemented in the C ++ Builder.

SimulPhysics [START_REF] Aristizábal | Uso de las NTIC para apoyar la enseñanza de la física básica para ingenieros: experiencia en la Universidad Nacional de Colombia sede Medellín[END_REF] is a multiplatform software, written in Java, that contains a set of interactive simulations on different topics of physics that supports the assimilation of the concepts developed in the theoretical part of the course. On the platform, we can find simulations on different topics that are seen in the subject of Basic Physics, e.g., mechanics, optics, waves, etc.

Celosías paramétricas [START_REF] Mansilla | Celosías Paramétricas[END_REF] is an Android application, targeted to Education and Learning, that automatically generates trusses by combining parameters related to the geometry, material, supports and loads, and shows the influence of each of them on its structural behavior: reaction forces, stresses, and deformations. It allows the optimization of the structural design by determining the values that produce the lowest stresses and deformations for the highest loads, while minimizing the weight.

MDSolids [START_REF] Philpot | MDSolids[END_REF] is an educational software package dedicated to the introduction of the Mechanics of Materials course. The software has an intuitive interface that guides the user to understand important facts that are found in different types of problems. It helps to visualize the nature of internal stresses and strains.

Regarding assessment, Hernández Pina [START_REF] Pina | La evaluación de los alumnos en el contexto de la evaluación de la calidad de las universidades[END_REF] investigates the learning process in the university context, and concludes by stating that the way in which the teacher raises the evaluation affects learning approaches (superficial or deep) and the quality of such learning. Quantitative evaluation strategies lead to superficial learning approaches, while formative and qualitative strategies can produce deep learning and high-performance approaches [START_REF] Sadler | Formative assessment and the design of instructional systems[END_REF]. In [START_REF] Restrepo | Generador automático de tareas como apoyo a los procesos de evaluación[END_REF], a project carried out by the School of Engineering of the EAFIT University implements an interactive software in a virtual platform. With this tool the students can carry out tests, clarifying doubts and preparing them for the final test. It incorporates a rating system and a problem generator algorithm for the course of Statics.

Following the development of the new Information and Communication Technologies (ICTs), we propose a software tool for a virtual educational environment which could be used to improve the assessment of problem-based subjects and provide a standardized marking system that could be used across classes and campuses. The software generates problems on the different topics of the subject Mechanics of Materials, giving the student an interactive tool that allows him to learn by solving problems that help him to strengthen his analytical skills and increase his cognitive abilities. Moreover, the teacher could perform controls and follow performance indicators, which can be used as a reference of the level of study and understanding of their students on the content seen in class, and the competences defined in the curriculum.

The paper is organized as follows: in the next section, we exposed the methodology, including the knowledge base and performance indicators, which are used to define the constraints for the design of the tool. In section 3, we present details related to the implementation of the tool in the Learning Management System. Conclusions are given in the final section.

Methods

The methodology is focused on the subject Mechanics of Materials, although it could be easily extended to other problem-based subjects in engineering curricula. The project is based on the framework of Educational Technology, which refers to the systematic way of conceiving, applying and evaluating the set of teaching and learning processes, considering both technical and human resources and the interactions between them, as a way to obtain a more effective education. In this regard, we define a strategy to tackle different issues detected during the learning and assessment process in the schools of engineering at Universidad Industrial de Santander. First, when evaluating problem-based subjects, the teacher should pay attention to the process followed by the student to reach the answer. Normally, the problem solution is divided into steps, and a weighting is given to each of the partial steps. However, if the student fails early in the process, the teacher should evaluate the remaining steps considering the wrong values obtained at the beginning, and their impact in the final answer. This is particularly difficult for large groups, which finally yields nonhomogeneous, subjective marks. This difference is accentuated if we consider groups led by different teachers. Moreover, regarding teachers, there is a notable difference in the assessment performed by tenured and non-tenured professors, mostly due to the correlation between the grades and the teacher evaluation, which affects contract extensions for non-tenured professors. There are even instances where the content given by the teacher is incomplete and we end up with groups of students with different level of competences reached. This non-standardized marking system gives rise to problems for the mobility of students across campuses, and helps to increase retention rates.

We start by defining the knowledge base that encompasses the area of Mechanics of Materials, putting special attention on the competences required at the end of the course. Next, we work on the performance indicators for the established competences, which help to define the constraints for the design of the tool.

Subject description and knowledge by competences

The subject Mechanics of Materials [START_REF] Beer | Mechanics of Materials[END_REF] studies the mechanics of deformable solids by simplified linear elastic models, relating in the design process the parameters of geometry, forces, and material. The mechanical engineer works with static and moving mechanical structures and systems. The design and calculation of these systems involve considering all the loads that act on each of its elements and the intensity of internal force (stresses) generated by these loads. In mechanical design, the stresses can be related to the mechanical strength of the material to assess the possibility of failure for a given element. In general, a machine element can fail because the stresses exceed the strength of the material, or because excessive deformation is occurring that prevents its normal operation [START_REF] Ayestarán | Resistencia de vigas esbeltas de acero inoxidable bajo cargas concentradas mediante elementos finitos[END_REF]. The calculation of these stresses and deformations involves considering the internal forces and the geometric features of the element.

At the end of the course the student is expected to:

• Recognize the type of stress produced by a force or moment applied to a machine or structural element, according to the way it acts on the element. • Calculate stresses and deformations that occur at any point of an element subjected to a specific type of load. • Select the critical point where the highest stress occurs in any element subjected to several types of load simultaneously. • Define the plane stress state at the critical point of any element subjected to different types of external load. • Analyze the possible failure of an element subjected to different types of loads using different criteria according to the type of material, and propose changes to the geometry accordingly. The curriculum is divided into four fundamental thematic blocks. The first block addresses the introduction to the concepts of stress and strain, and the state of axial loading. The second block deals with elements subjected to pure torsion. The third block describes pure bending, typical of beams. The fourth block refers to multiaxial stress, combined load, pressure vessels and design failure criteria.

Performance indicators

The performance of the student is evaluated with quantifiable indicators, related to the competences, which provide a global assessment of his understanding of the subject. These are:

• Identifies the type of load produced by a force or a moment applied on a machine or structural element and evaluates the corresponding stress and deformation. • Recognizes the differences between the different types of stresses and deformations that appear in elastic problems. • Identifies the critical point where the highest stresses occur in any mechanical element subjected to several types of load simultaneously. • Calculates the principal stresses, maximum shear stress and equivalent normal stress at the critical point of any element subjected to different types of external load, according to a specified failure theory.

Design of the tool

In this section, the thematic and evaluative content of the course is identified, as well as the parameters that can be varied to generate multiple versions of the same exercises without modifying their knowledge structure. To carry out this process the following steps were followed: Definition of the learning objective: to be able to define the exercise that is going to be proposed, we must know what we want to achieve with it, i.e., identify the cognitive competences that the student must acquire during the course. Here, we refer back to the competences listed in section 2.1.

Definition of the exercises: a review of the curriculum of the subject Mechanics of Materials must be done, and the main topics of the subject are identified: axial loading, bending, torsion and combined load. After considering the thematic content, the bibliography presented by the teachers is reviewed in order to look for the most representative exercises. For each topic, a set of exercises is proposed in order to cover the main knowledge components of the cognitive competences.

Solution and generalization of the exercise: the selected exercises must be solved to create a base algorithm that follows the steps taken to obtain the answer, in order to generate several versions of the same problem by modifying certain parameters but preserving the knowledge component. We define and bound the parameters that can vary to enforce randomness in the data of the problem. In Table 1 we can see an example of the design of an exercise that was included in the software database.

The problem template in Table 1 considers the following points: (i) the topic and objective of the exercise describing the knowledge component, (ii) the problem statement where the questions to be answered are proposed, (iii) the list of parameters acting as variables of the problem, i.e., the problem data randomly generated by the software that must have a previously defined range to avoid errors in the answers, (iv) a general solution for the exercise is developed, which allows to code the exercise easily.

Software implementation

There are many advantages of developing the program as web software. It is multiplatform and the student can use it from any electronic device. This approach is convenient as it makes the tool accessible from virtually anywhere, in agreement with mobile learning (m-learning) statement of students using mobile devices at their time convenience. The tool allows the teacher to use it to do quizzes and workshops. Any student will be able to propose new exercises, facilitating communication with the teacher and feeding the database to increase the knowledge base.

Due to the requirements of the software, especially the one focusing on facilitating the process for users with poor programming experience, there is a need to use an expression analyzer. An expression analyzer is an algorithm that interprets a mathematical expression with pre-established functions, assigns values to each of the variables that intervene in the expression, and gives a numerical answer, see Figure 1. This condition speed up the task of introducing new exercises to the database, since instead of programming the algorithm, the existing one can be quickly modified to the required needs. Thus, the expression analyzer is used to separate the coding of the exercises from the main routines, without it, the software was conditioned to a person making a particular code for each of the exercises.

Table 1. Problem template of an exercise for axial loading.

Topic: Axial loading. Objective: Determination of normal stress in the two-forces structural elements.

Problem statement:

There is a machine consisting of two tongs that support a cylindrical element, as shown in the figure .  • Determine the normal stress developed in the link AB. • Determine the shear stress developed on pin A, consider that it is subjected to double shear. Note: Units are given in international system, SI. After performing a search on the web, we found a previously coded expression analyzer with Lesser General Public License (LGPL) [START_REF] Moreno Parra | Evaluador de Expresiones Algebraicas II[END_REF]. This analyzer allows simplifying the template definition for the proposed problems. Equations obtained during the solution steps are introduced in the template to parametrize the solution of the problem with random input.

Parameters

The software was implemented in the proprietary platform COMA (http://calumet.github.io/) [START_REF] Calumet | COMA[END_REF] used as the Learning Management System (LMS) at the School of Engineering. The software was codified in the programming languages Java and JavaScript, to follow a multiplatform (Linux, Windows, Android, Mac) approach. In addition, the HTML language was used for the presentation of the contents. The database of exercises can be constantly updated since anyone with access to the platform can propose problems. The database contains different types of questions, i.e., multiple choice, calculated, true/false, matching, and essay. As indicated before, the focus is on the development of calculated questions for problems centered in the process. One key feature that is proposed in this implementation is that values are recalculated with the input values given in previous steps during the process of solving the problem. This segmentation of the process allows assessing the competences of the student independently. For the specific case of problems related to Mechanics of Materials, it evaluates the statics, the geometry, the stress and the strains calculations, all in different steps, regardless of whether the student failed a previous step. The software offers two services: (i) Training: The teacher creates a workshop by setting the number of problems and the stipulated time to develop them. As the software has the option to review answers, the student can solve the exercise as many times as he wants, monitoring their results, until he is sure the procedure is correct. This improved feedback is required for effective formative assessment [START_REF] Sadler | Formative assessment and the design of instructional systems[END_REF], and it is extremely difficult with mainstream education, (ii) Evaluation: The teacher previously creates an exam, setting a number of problems and an adequate time to solve them, and takes the students to a computer room in order for the students to carry out the activity.

Figure 2 shows a flow diagram with the internal process that is carried out when the teacher introduces a query of the type "create questions with expressions" in the system. For a teacher to include a question in the software, he must add the statement and the problem image, followed by the data and the mathematical equations to develop the exercise, then he chooses which expressions will be used as answer. Figure 3 shows snippets for the definitions of variables and expressions in the code. When the student works on a workshop and presses the review option, the program follows a series of internal steps to show him if the answers are correct. This sequence consists in taking the randomly generated data and calculating the solution to the defined mathematical equations. After that, the software compares the values of the expressions indicated as answers with the values entered by the user, showing in red the boxes with wrong answers and in green the boxes with correct answers.

Conclusions

We presented the development of a multiplatform software tool to be used for the training and assessment of students of the subject Mechanics of Materials. The tool, implemented in Java, JavaScript and HTML, randomly generates exercises based on a parametric definition of the problems. The solution is calculated and compared with the answer given by the students and can be used to train or mark their performance. The problems are segmented into steps based on the knowledge base and the competences to be evaluated. The software is initially focused on the subject Mechanics of Materials, but it is easily extensible to other problem-based engineering subjects. By reviewing the answers, the student can check their errors and repeat the process in order to obtain effective feedback in the formative assessment process. This tool is based in innovative teaching strategies that seek to favor a teaching based on the active participation of the student. Moreover, it provides a standardized marking system that could be used across campuses and groups, reducing heterogeneities for large numbers of students due to subjective marking and non-uniform content.

Figure 1 .

 1 Figure 1. Basic algorithm of the expression analyzer.

Figure 2 .

 2 Figure 2. Basic algorithm for "create questions with expressions".

Acknowledgments

We acknowledge the support given by the group Calumet, Alejandro Riaño, and Luis Ignacio González Ramírez.

[{ "type": "variable", "variable": "b", "description": "L [m]", "min": "1.5", "max": "2.8", "value": "" }, … { "type": "variable", "variable": "z", "description": "Area bar AB [m]", "min": "0.001", "max": "0.002", "value": "" }, … { "type": "expresion", "variable": "g", "expression": "a/(2*sen(p*e/180))", "description": "R", "answer": false, "value": "" }, … { "type": "expression", "variable": "l", "expression": "k/z", "description": "Admissible shear stress AB", "answer": true, "value": "" }] Figure 3. Snippets of the array definition in JSON for the problem in Table 1, for variables and expressions.