Olgeir Sigmarsson

INFLUENCE OF TECTONISM ON THE COMPOSITION OF ACID AND BASALTIC LAVA

The Neovolcanic zones in Iceland are of threefold origin: the Mid-Atlantic rift-zone, the Snaefellsnes "leaky-transform fault", and South-Iceland transgressive volcanic zone or propagating rift segment. Basalt compositions range from alkali basalts through transitional basalts to tholeiites along the off-rift zones towards the centre of Iceland. Diminishing proportion of garnet-pyroxenite melts towards the centre of the mantle plume is the preferred explanation. The origin of silicic magmas reflects their tectonic settings and appears to be controlled by the crustal geothermal gradient.

Basaltic lavas

Mantle heterogeneity beneath oceanic islands can be inferred from basaltic bulk lava samples or glass inclusions. Lava samples represent the integrated magma formation processes (melting, mixing, contamination and fractionation) whereas glass inclusions record melts less affected by higher-level processes. Nevertheless, more precise analyses of compositional parameters can be obtained on bulk samples that permits thorough assessment of mantle heterogeneity from the final basalt composition at surface. This is especially true for well characterized sample suites from oceanic islands showing important compositional variability in their basalts. Here, we present new major-, trace element and Sr-Nd-Hf-Pb isotope results on post-glacial basalts from the Neovolcanic belt of Iceland: the Snaefellsnes Volcanic Zone (SNVZ), the South Iceland Volcanic Zone (SIVZ) and the Mid-Iceland Volcanic Zone (MIVZ). The results reveal systematic lateral variations. Along the strike of the off-rift zones (SNVZ and SIVZ), regular decrease is observed for alkalinity, La/Yb, Sm/Yb, HFSE/Y, 87 Sr/ 86 Sr and 206-207- 208 Pb/ 204 Pb towards the centre of Iceland, where olivine-tholeïtes with depleted trace element compositions, low 87 Sr/ 86 Sr and 206-207-208 Pb/ 204 Pb (and high 143 Nd/ 144 Nd and 176 Hf/ 177 Hf) are produced. These latter have similar composition as most mafic lavas along the rift zones that cut through the island from SW to NE. A clear compositional difference is thus observed for basalts erupting at the spreading axis compared to those from the non-rifting flank zones.

The systematic decrease of residual garnet signature (high La/Yb, Sm/Yb, HFSE/Y) and Sr and Pb isotope ratios (and the increase of Nd and Hf isotope ratios) in the off-rift lavas towards the assumed centre of the plume at Iceland's centre, is consistent with the occurrence of fertile garnet pyroxenites in the mantle source. Basalts from the periphery of the flank zones are probably generated by less total melting and thus sample the more enriched and fusible component(s) of the mantle beneath Iceland. Melts of this fertile mantle thus dominate the bulk melt composition at the periphery, and are progressively diluted towards the centre. In the Mid-Iceland Volcanic Zone, basalt compositions are dominated by melts derived principally from the more refractory component(s) of the mantle. The lateral variations of the basalt compositions and the proposed diminishing role of fertile pyroxenites towards the centre of the island is readily explained by increased total melting towards the hotter core of the Iceland mantle plume. Finally, the first-order observation of strong correlations between geochemical parameters and geographical location of the basaltic craters strongly suggests that mantle heterogeneity dominates the highly-incompatible trace element and radiogenic isotope ratios relative to other magma processes, such as crustal contamination and fractional crystallization forming the final basalt composition.

In-situ analysis of major-and trace elements in Holocene tephra produced during subglacial eruptions confirm the spatial variability of the basalt compositions as recorded in the lava samples . Furthermore, cross-correlations of the compositional variability permits to identify the very centre of the Iceland hot-spot, namely the most active Icelandic volcano, Grímsvötn. Moreover, eruption frequency of the principal tephra-forming volcanoes can be determined for the last 6 kyears, where robust age models for soil accumulation can be obtained from well-dated silicic key-layers.

Acid lavas

Pleistocene and Holocene peralkaline rhyolites from Torfajökull (South Iceland Volcanic Zone) and Ljósufjöll central volcanoes and trachytes from Snaefellsjökull (Snaefellsnes Volcanic Zone) allow the assessment of the mechanism for silicic magma genesis as a function of tectonic settings and variable crustal geothermal gradient. The low 18O (2.4‰) and low Sr concentration (12.2 ppm) measured in Torfajökull rhyolites are best explained by partial melting of hydrated metabasaltic crust followed by major fractionation of feldspar. In contrast, very high 87Sr/86Sr (0.70473) and low Ba (8.7 ppm) and Sr (1.2 ppm) concentrations measured in Ljósufjöll silicic lavas are best explained by fractional crystallisation and subsequent 87Rb decay. Snaefellsjökull trachytes are also principally generated by fractional crystallisation, with only small crustal contamination. The fact that silicic magmas within, or close to, the rift zone are principally generated by crustal melting whereas those from off-rift zones are better explained by fractional crystallisation clearly illustrates the controlling influence of the thermal state of the crust on silicic magma genesis in Iceland.

The origin of the Quaternary silicic rocks in Iceland thus appears to be linked to the thermal state of the crust, which in turn depends on the regional tectonic settings. This simple model can be tested on rocks from the Miocene to present, both to suggest an internally consistent model for silicic magma formation in Iceland and to constrain the link between tectonic settings and silicic magma petrogenesis. New major and trace element compositions together with O-, Sr-and Nd-isotope ratios have been obtained on silicic rocks from 19 volcanic systems ranging in age from 13 Ma to present (see sample location in Fig. 1). This allows tracing spatial and temporal evolution of both magma generation and the corresponding sources. Low δ 18 O (<5 ‰ SMOW) in most silicic rocks results from partial melting of hydrothermally altered metabasaltic crust, whereas low Ba and Sr concentrations (down to 9 and 1 ppm, respectively), in addition to Th concentrations higher than 9 ppm in the silicic rocks, indicate important role of fractional crystallisation during the final stage of magma formation. Trace element ratios, such as Th/U, Th/Zr and Th/La, record important role of accessory minerals during final differentiation. The 143 Nd/ 144 Nd proves to be an excellent marker of the silicic-magma source: high 143 Nd/ 144 Nd (0.51303 to 0.51296) characterizes a "rift-zone source", from which the silicic magmas are generated by crustal anatexis, whereas low 143 Nd/ 144 Nd (0.51290 to 0.51297) is typical of an "off-rift source", where rhyolites formed by fractional crystallisation of mantle-derived basaltic magma in a cooler environment far from the rift-zone.

The spatial and temporal distribution of samples having a "rift zone" Nd-isotope signature allows inferences to be drawn about the past tectonic regime. At Snaefellsnes Peninsula, silicic magmas younger than 5.5 Ma have an "off-rift" Nd-isotope signature, whereas those older than 6.8 Ma have typical "rift zone" characteristics. This suggests that before 6.8 Ma silicic rocks were generated in a riftzone by crustal melting caused by high geothermal gradient. But later than 5.5 Ma they were produced in a flank zone environment by fractional crystallisation alone, probably due to decreasing geothermal gradient, of basalts derived from a mantle source with lower 143 Nd/ 144 Nd. This is in agreement with an eastwards rift-jump, from Snaefellsnes towards the present Reykjanes Rift Zone , between 7 and 5.5 Ma. In the South Iceland Volcanic Zone (SIVZ), the intermediate Nd-signature observed in silicic rocks from the Torfajökull central volcano reflects the transitional character of the basalts erupted at this propagating rift segment. Therefore, the abundant evolved rocks at this major silicic complex result from partial melting of the transitional alkaline basaltic crust (<3 Ma) generated in this propagating rift segment, suggesting rapid crustal recycling there. The light and dark grey fields represent the ""rift-zone" and "off-rift" values, respectively (see text for further discussion). Note that all samples from the Snaefellsnes Peninsula (SNVZ) plot in the "off-rift" field with the exception of Drap 2. All other samples from this study plot in the "rift-zone" field with the exception of the samples from the SIVZ that fall in between. The DRAP 2 is a 6.8 Ma dacite that formed in a former rift-zone that was active from 15-6 Ma, hence its higher 143 Nd/ 144 Nd compared to other volcanics from the SNVZ . For comparison, range of 143 Nd/ 144 Nd from Öraefajökull volcanic system is also shown (Prestvik et al., 2001). The error bar (2σ) corresponds to the external error of our standard over the duration of this project.

Improved understanding of silicic magmatism in Iceland can, therefore, be used for deciphering past geodynamic settings characterized by rift-and off-rift zones resulting from interaction of a mantle plume and divergent plate boundaries.

Figure 1

 1 Figure 1 Map of Iceland showing sample locations. The Neovolcanic zones are the Snaefellsnes Volcanic Zone (SNVZ), Reykjanes Rift Zone (RRZ), Mid-Iceland Volcanic Zone (MIVZ), North-Iceland Rift Zone (NIRZ), South-Iceland Volcanic Zone (SIVZ) and the Öraefi Volcanic Belt (OVB). White areas correspond to the main glaciers.

Figure 2 .

 2 Figure2. The isotope ratio143 Nd/144 Nd vs. Th concentration in lavas from Iceland. The light and dark grey fields represent the ""rift-zone" and "off-rift" values, respectively (see text for further discussion). Note that all samples from the Snaefellsnes Peninsula (SNVZ) plot in the "off-rift" field with the exception of Drap 2. All other samples from this study plot in the "rift-zone" field with the exception of the samples from the SIVZ that fall in between. The DRAP 2 is a 6.8 Ma dacite that formed in a former rift-zone that was active from 15-6 Ma, hence its higher143 Nd/ 144 Nd compared to other volcanics from the SNVZ . For comparison, range of 143 Nd/ 144 Nd from Öraefajökull volcanic system is also shown(Prestvik et al., 2001). The error bar (2σ) corresponds to the external error of our standard over the duration of this project.


Marc Mezzarobba
email: marc.mezzarobba@inria.fr

NumGfun: a Package for Numerical and Analytic Computation with D-finite Functions

Keywords: Descriptors: I.1.2 [Symbolic and Algebraic Manipulation]: Algorithms Algorithms, Experimentation, Theory D-finite functions, linear differential equations, certified numerical computation, bounds, Maple

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Support for computing with D-finite functions, that is, solutions of linear differential equations with polynomial coefficients, has become a common feature of computer algebra systems. For instance, Mathematica now provides a data structure called DifferentialRoot to represent arbitrary D-finite functions by differential equations they satisfy and initial values. Maple's DESol is similar but more limited.

An important source of such general D-finite functions is combinatorics, due to the fact that many combinatorial structures have D-finite generating functions. Moreover, powerful methods allow to get from a combinatorial description of a class of objects to a system of differential equations that "count" these objects, and then to extract precise asymptotic information from these equations, even when no explicit counting formula is available [START_REF] Flajolet | Analytic Combinatorics[END_REF][START_REF] Stanley | Enumerative combinatorics[END_REF]. A second major application of D-finiteness is concerned with special functions. Indeed, many classical functions of mathematical physics are D-finite (often by virtue of being defined as "interesting" solutions of simple differential equations), which allows to treat them uniformly in algorithms. This is exploited by the Encyclopedia of Special Functions [START_REF] Meunier | ESF: an automatically generated encyclopedia of special functions[END_REF] and its successor under development, the Dynamic Dictionary of Mathematical Functions [20], an interactive computergenerated handbook of special functions.

These applications require at some point the ability to perform "analytic" computations with D-finite functions, starting with their numerical evaluation. Relevant algorithms exist in the literature. In particular, D-finite functions may be computed with an absolute error bounded by 2 -n in n log O (1) n bit operations-that is, in softly linear time in the size of the result written in fixed-point notation-at any point of their Riemann surfaces [START_REF] Chudnovsky | Computer algebra in the service of mathematical physics and number theory[END_REF], the necessary error bounds also being computed from the differential equation and initial values [START_REF] Van Der Hoeven | Fast evaluation of holonomic functions[END_REF]. However, these algorithms have remained theoretical [13, §9.2.1]. The ability of computer algebra systems to work with D-finite functions is (mostly) limited to symbolic manipulations, and the above-mentioned fast evaluation algorithm has served as a recipe to write numerical evaluation routines for specific functions rather than as an algorithm for the entire class of D-finite functions.

This article introduces NumGfun, a Maple package that attempts to fill this gap, and contains, among other things, a general implementation of that algorithm. NumGfun is distributed as a subpackage of gfun [START_REF] Salvy | Gfun: A Maple package for the manipulation of generating and holonomic functions in one variable[END_REF], under the GNU LGPL. Note that it comes with help pages: the goal of the present article is not to take the place of user documentation, but rather to describe the features and implementation of the package, with supporting examples, while providing an overview of techniques relevant to the development of similar software. The following examples illustrate typical uses of NumGfun, first to compute a remote term from a combinatorial sequence, then to evaluate a special function to high precision near one of its singularities.

Example 1. The Motzkin number Mn is the number of ways of drawing non-intersecting chords between n points placed on a circle. Motzkin numbers satisfy (n + 4)Mn+2 = 3(n + 1)Mn + (2n + 5)Mn+1. Using NumGfun, the command nth_term({(n+4)*M(n+2) = 3*(n+1)*M(n)+(2*n+5)*M(n+1), M(0)=1,M(1)=1},M(n),k ) computes M 10 5 = 6187 . . . 7713 ≃ 10 47 706 in1 4.7 s and M 10 6 = 2635 . . . 9151 ≃ 10 477 112 in 1 min. Naïvely unrolling the recurrence (using Maple) takes 10.7 s for M 10 5 , and 41 s for M 2•10 5 . On this (non-generic) example, nth_term could be made competitive for smaller indices by taking advantage of the fact that the divisions that occur while unrolling the recurrence are exact.

Example 2. The double confluent Heun function

U α,β,γ,δ satisfies (z 2 -1) 3 U ′′ (z) + (2z 5 -4z 3 -αz 4 + 2z + α)U ′ (z) + (βz 2 + (2α + γ)z + δ)U (z) = 0, U (0) = 1, U ′ (0) = 0. It is singular at z = ±1.
The command evaldiffeq(eq,y(z), -0.99,1000) where eq is this differential equation yields U 1, 1 3 , 1 2 ,3 (-0.99) ≈ 4.67755...(990 digits)...05725 in 22 s.

Related work. Most algorithms implemented in NumGfun

originate in work of Chudnovsky & Chudnovsky and of van der Hoeven. Perhaps the most central of these is the "bit burst" numerical evaluation method [START_REF] Chudnovsky | Computer algebra in the service of mathematical physics and number theory[END_REF]. It belongs to the family of binary splitting algorithms for D-finite series, hints at the existence of which go back to [2, §178], and generalizes earlier work of Brent [START_REF] Brent | The complexity of multiple-precision arithmetic[END_REF] for specific elementary functions. Better known (thanks among others to [START_REF] Haible | Fast multiprecision evaluation of series of rational numbers[END_REF]) binary splitting algorithms can be seen as special cases of the bit burst algorithm. One reason why, unlike these special cases, it was not used in practice is that in [START_REF] Chudnovsky | Computer algebra in the service of mathematical physics and number theory[END_REF], none of the error control needed to ensure the accuracy of the computed result is part of the algorithm. Van der Hoeven's version [START_REF] Van Der Hoeven | Fast evaluation of holonomic functions[END_REF] addresses this issue, thus giving a full-fledged evaluation algorithm for the class of D-finite functions, as opposed to a method to compute any D-finite function given certain bounds.

These algorithms extend to the computation of limits of D-finite functions at singularities of their defining equation. The case of regular singularities is treated both in [START_REF] Chudnovsky | On expansion of algebraic functions in power and Puiseux series (I, II)[END_REF][START_REF] Chudnovsky | Computer algebra in the service of mathematical physics and number theory[END_REF], and more completely in [START_REF] Van Der Hoeven | Fast evaluation of holonomic functions near and in regular singularities[END_REF], that of irregular singular points in [START_REF] Van Der Hoeven | Efficient accelero-summation of holonomic functions[END_REF]. See [4, §12.7], [35, §1] for more history and context.

On the implementation side, routines based on binary splitting for the evaluation of various elementary and special functions are used in general-purpose libraries such as CLN [START_REF] Haible | CLN[END_REF] and MPFR [START_REF] Fousse | MPFR: A multiple-precision binary floating-point library with correct rounding[END_REF][START_REF] Jeandel | Évaluation rapide de fonctions hypergéométriques[END_REF]. Binary splitting of fast converging series is also the preferred algorithm of software dedicated to the high-precision computation of mathematical constants on standard PCs, including the current record holder for π [START_REF] Bellard | Computation of 2700 billion decimal digits of Pi using a Desktop Computer[END_REF]. Finally, even the impressive range of built-in functions of computer algebra systems is not always sufficient for applications. Works on the implementation of classes of "less common" special functions that overlap those considered in NumGfun include [START_REF] Abad | An algorithm to obtain global solutions of the double confluent Heun equation[END_REF][START_REF] Falloon | Theory and computation of spheroidal wavefunctions[END_REF].

This work is based in part on the earlier [START_REF] Mezzarobba | Génération automatique de procédures numériques pour les fonctions D-finies[END_REF].

Contribution.

The main contribution presented in this article is NumGfun itself. We recall the algorithms it uses, and discuss various implementation issues. Some of these descriptions include improvements or details that do not seem to have appeared elsewhere. Specifically: (i) we give a new variant of the analytic continuation algorithm for D-finite functions that is faster with respect to the order and degree of the equation; (ii) we improve the complexity analysis of the bit burst algorithm by a factor log log n; (iii) we point out that Poole's method to construct solutions of differential equations at regular singular points can be rephrased in a compact way in the language of noncommutative polynomials, leading to faster evaluation of D-finite functions in these points; and (iv) we describe in some detail the practical computation of all the bounds needed to obtain a provably correct result.

What NumGfun is not. Despite sharing some of the algorithms used to compute mathematical constants to billions of digits, our code aims to cover as much as possible of the class of D-finite functions, not to break records. Also, it is limited to "convergent" methods: asymptotic expansions, summation to the least term, and resummation of divergent power series are currently out of scope.

Terminology. Like the rest of gfun, NumGfun works with D-finite functions and P-recursive sequences. We recall only basic definitions here; see [30, §6.4] for further properties. A formal power series y ∈ [[z]] is D-finite over K ⊆ if it solves a non-trivial linear differential equation

y (r) (z) + ar-1(z) y (r-1) (z) + • • • + a0(z) y(z) = 0 (1)
with coefficients a k ∈ K(z). The same definition applies to analytic functions. A sequence u ∈ AE is P-recursive over K if it satisfies a non-trivial linear recurrence relation

un+s + bs-1(n) un+s-1 + • • • + b0(n) un = 0, b k ∈ K(n). ( 2
)
A sequence (un)n∈AE is P-recursive if and only if its generating series n∈AE unz n is D-finite. The poles of the coefficients a k of (1) are its singular points; nonsingular points are called ordinary. In gfun, a D-finite function is represented by a differential equation of the form [START_REF] Abad | An algorithm to obtain global solutions of the double confluent Heun equation[END_REF] and initial values at the origin, which we assume to be an ordinary point. Similarly, P-recursive sequences are encoded by a recurrence relation plus initial values, as in Ex. 1 above. If

y(z) = ∞ k=0 y k z k , we let y;n(z) = n-1 k=0 y k z k and yn;(z) = ∞ k=n y k z k .
The height of an object is the maximum bit-size of the integers appearing in its representation: the height of a rational number p/q is max(⌈log p⌉, ⌈log q⌉), and that of a complex number (we assume that elements of number fields É(ζ) are represented as i xiζ i /d with xi, d ∈ ), polynomial, matrix, or combination thereof with rational coefficients is the maximum height of its coefficients. We assume that the bit complexity M (n) of n-bit integer multiplication satisfies M (n) = n(log n) O (1) , M (n) = Ω(n log n), and M (n + m) ≥ M (n) + M (m), and that s × s matrices can be multiplied in O(s ω ) operations in their coefficient ring.

BINARY SPLITTING

"Unrolling" a recurrence relation of the form (2) to compute u0, . . . , uN takes Θ(N 2 M (log N )) bit operations, which is almost linear in the total size of u0, . . . , uN , but quadratic in that of uN . The binary splitting algorithm computes a single term uN in essentially linear time, as follows: (2) is first reduced to a matrix recurrence of the first order with a single common denominator:

q(n)Un+1 = B(n)Un, B(n) ∈ [n] s×s , q(n) ∈ [n], (3) so that UN = P (0, N ) U0/ N-1 i=0 q(i) , where P (j, i) = B(j -1) • • • B(i + 1)B(i).
One then computes P (0, N ) recursively as P (0, N ) = P (⌊N/2⌋ , N )P (0, ⌊N/2⌋), and the denominator N-1 i=0 q(i) in a similar fashion (but separately, in order to avoid expensive gcd computations).

The idea of using product trees to make the most of fast multiplication dates back at least to the seventies [4, §12.7]. The general statement below is from [12, Theorem 2.2], except that the authors seem to have overlooked the cost of evaluating the polynomials at the leaves of the tree.

Theorem 1 (Chudnovsky, Chudnovsky). Let u be a P-recursive sequence over É(i), defined by (2). Assume that the coefficients b k (n) of (2) have no poles in AE. Let d, h denote bounds on their degrees (of numerators and denominators) and heights, and d ′ , h ′ corresponding bounds for the coefficients of B(n) and q(n) in [START_REF] Bellard | Computation of 2700 billion decimal digits of Pi using a Desktop Computer[END_REF]. Assuming N ≫ s, d, the binary splitting algorithm outlined above computes one

term uN of u in O(s ω M (N (h ′ + d ′ log N )) log(N h ′ )), that is, O(s ω M (sdN (h + log N )) log(N h)), bit operations. Proof sketch. Write H = h ′ + d ′ log N .
Computing the product tree P (0, N ) takes O(s ω M (N H) log N ) bit operations [12, §2] (see also Prop. 1 below), and the evaluation of each leaf B(i) may be done in O(M (H) log d ′ ) operations [5, §3.3]. This gives uN as a fraction that is simplified in ) and get rid of some terms, so that the total complexity simplifies as stated.

O(M (N H) log(N H)) operations [8, §1.6]. Now consider how (2) is rewritten into (3). With co- efficients in [i] rather than É(i), the b k (n) have height h ′′ ≤ (d + 1)h. To get B(n) and q(n),
Since the height of uN may be as large as Ω((N +h) log N ), this result is optimal with respect to h and N , up to logarithmic factors. The same algorithm works over any algebraic number field instead of É(i). This is useful for evaluating D-finite functions "at singularities" ( §4). More generally, similar complexity results hold for product tree computations in torsion-free -algebras (or É-algebras: we then write A = É ⊗ A ′ for some -algebra A ′ and multiply in × A ′ ), keeping in mind that, without basis choice, the height of an element is defined only up to some additive constant.

Constant-factor improvements. Several techniques permit to improve the constant hidden in the O(•) of Theorem 1, by making the computation at each node of the product tree less expensive. We consider two models of computation.

In the FFT model, we assume that the complexity of long multiplication decomposes as M (n) = 3F (2n)+O(n), where F (n) is the cost of a discrete Fourier transform of size n (or of another related linear transform, depending on the algorithm). FFT-based integer multiplication algorithms adapt to reduce the multiplication of two matrices of height n in s×s to O(n) multiplications of matrices of height O(1), for a total of O(s 2 M (n) + s ω n) bit operations. This is known as "FFT addition" [START_REF] Bernstein | Fast multiplication and its applications[END_REF], "computing in the FFT mode" [START_REF] Van Der Hoeven | Fast evaluation of holonomic functions[END_REF], or "FFT invariance". A second improvement ("FFT doubling", attributed to R. Kramer in [4, §12.8]) is specific to the computation of product trees. The observation is that, at an internal node where operands of size n get multiplied using three FFTs of size 2n, every second coefficient of the two direct DFTs is already known from the level below.

The second model is black-box multiplication. There, we may use fast multiplication formulae that trade large integer multiplications for additions and multiplications by constants. The most obvious example is that the products in É(i) may be done in four integer multiplications using Karatsuba's formula instead of five with the naïve algorithm. In general, elements of height h of an algebraic number field of degree d may be multiplied in 2dM (h)+O(h) bit operations using the Toom-Cook algorithm. The same idea applies to the matrix multiplications. Most classical matrix multiplication formulae such as Strassen's are socalled bilinear algorithms. Since we are working over a commutative ring, we may use more general quadratic algorithms [9, §14.1]. In particular, for all s, Waksman's algorithm [START_REF] Waksman | On Winograd's algorithm for inner products[END_REF] multiplies s × s matrices over a commutative ring R using s 2 ⌈s/2⌉ + (2s -1)⌊s/2⌋ multiplications in R, and Makarov's [START_REF] Makarov | An algorithm for multiplying 3 × 3 matrices[END_REF] multiplies 3 × 3 matrices in 22 scalar multiplications. These formulas alone or as the base case of a Strassen scheme achieve what seems to be the best known multiplication count for matrices of size up to 20.

Exploiting these ideas leads to the following refinement of Theorem 1. Similar results can be stated for generalalgebras, using their rank and multiplicative complexity [START_REF] Bürgisser | Algebraic complexity theory[END_REF]. Proposition 1. Let d ′ and h ′ denote bounds on the degrees and heights (respectively) of B(n) and q(n) in Eq. ( 3). As N, h ′ → ∞ (s and d ′ being fixed), the number of bit operations needed to compute the product tree Proof. Each node at the level k (level 0 being the root) of the tree essentially requires multiplying s×s matrices with entries in

P (0, N ) is at most C + o(1) M N (h ′ + d ′ log N ) log(N h ′ ) , with C = (2s 2 + 1)
[i] of height H k = N (h ′ + d ′ log N )/2 k+1
, plus denominators of the same height. In the FFT model, this may be done in (2s 2 +1)M (H k ) operations. Since we assume

M (n) = Ω(n log n), we have ⌈log N⌉ k=0 2 k M (H k ) ≤ 1
2 M (H0) (a remark attributed to D. Stehlé in [START_REF] Zimmermann | The bit-burst algorithm[END_REF]). Kramer's trick saves another factor 3 2 . In the black-box model, the corresponding cost for one node is (3 MM(s) + 1)M (H k ) with Karatsuba's formula. Stehlé's argument applies again.

Note that the previous techniques save time only for dense objects. In particular, one should not use the "fast" matrix multiplication formulae in the few bottom levels of product trees associated to recurrences of the form (3), since the matrices at the leaves are companion.

Continuing on this remark, these matrices often have some structure that is preserved by successive multiplications. For instance, let sn = n-1 k=0 u k where (un) satisfies (2). It is easy to compute a recurrence and initial conditions for (sn) and go on as above. However, unrolling the recurrences (2) and sn+1 -sn = un simultaneously as

      un+1 . . . un+s-1 un+s sn+1       =       1 0 . . . . . . 1 0 * * • • • * 0 1 0 • • • 0 1             un . . . un+s-2 un+s-1 sn       (4)
is more efficient. Indeed, all matrices in the product tree for the numerator of (4) then have a rightmost column of zeros, except for the value in the lower right corner, which is precisely the denominator. With the notation MM(s) of Proposition 1, each product of these special matrices uses MM(s) + s 2 + s + 1 multiplications, vs. MM(s + 1) + 1 for the dense variant. Hence the formula ( 4) is more efficient as soon as MM(s + 1) -MM(s) ≥ s(s + 1), which is true both for the naïve multiplication algorithm and for Waksman's algorithm (compare [START_REF] Zimmermann | The bit-burst algorithm[END_REF]). In practice, on Ex. 3 below, if one puts un = sn+1 -sn in (2, 3) instead of using (4), the computation time grows from 1.7 s to 2.7 s.

The same idea applies to any recurrence operator that can be factored. Further examples of structure in product trees include even and odd D-finite series (e.g., [8, §4.9.1]). In all these cases, the naïve matrix multiplication algorithm automatically benefits from the special shape of the problem (because multiplications by constants have negligible cost), while fast methods must take it into account explicitly. Remark 1. A weakness of binary splitting is its comparatively large space complexity Ω(n log n). Techniques to reduce it are known and used by efficient implementations in the case of recurrences of the first order [START_REF] Cheng | Time-and space-efficient evaluation of some hypergeometric constants[END_REF][START_REF] Gourdon | Binary splitting method[END_REF][START_REF] Haible | CLN[END_REF][START_REF] Bellard | Computation of 2700 billion decimal digits of Pi using a Desktop Computer[END_REF].

Implementation. The implementation of binary splitting in

NumGfun includes some of the tricks discussed in this section. FFT-based techniques are currently not used because they are not suited to implementation in the Maple language. This implementation is exposed through two user-level functions, nth_term and fnth_term, that allow to evaluate Precursive sequences (fnth_term replaces the final gcd by a less expensive division and returns a floating-point result). Additionally, gfun[rectoproc], which takes as input a recurrence and outputs a procedure that evaluates its solution, automatically calls the binary splitting code when relevant. Examples 1 and 3 illustrate the use of these functions on integer sequences and convergent series respectively.

Example 3. [7, §6.10] Repeated integration by parts of the integral representation of Γ yields for 0 < Re(z) < 1

Γ(z) = ∞ n=0 e -t t n+z z(z + 1) • • • (z + n) + ∞ t e -u u z-1 du.
Taking t = 29 3 , it follows that the sum 65000 n=0 un, where u0 = 87e -t and (3n + 4)un+1 = 3tun, is within 10 -10000 of Γ(1/3), whence Γ(1/3) ≃ 2.67893 . . . (9990 digits) . . . 99978. This computation takes 1.7 s.

HIGH-PRECISION EVALUATION OF D-FINITE FUNCTIONS

We now recall the numerical evaluation algorithms used in NumGfun, and discuss their implementation.

Let y(z) = n ynz n be a D-finite series with radius of convergence ρ at the origin. Let z ∈ É(i) be such that |z| < ρ and height(z) ≤ h. The sequence (ynz n ) is P-recursive, so that the binary splitting algorithm yields a fast method for the high-precision evaluation of y(z). Here "high-precision" means that we let the precision required for the result go to infinity in the complexity analysis of the algorithm. More precisely, (ynz n ) is canceled by a recurrence relation of height O(h). By Theorem 1, y(z) may hence be computed to the precision 10 -p in

O M N (h + log N ) log(N h) (5)
bit operations, where N is chosen so that |yN;(z)| ≤ 10 -p , i.e. N ∼ p/ log(ρ/|z|) if ρ < ∞, and N ∼ p/(τ log p) for some τ if ρ = ∞.

In practice, the numbers of digits of (absolute) precision targeted in NumGfun range from the hundreds to the millions. Accuracies of this order of magnitude are useful in some applications to number theory [START_REF] Chudnovsky | Computer algebra in the service of mathematical physics and number theory[END_REF], and in heuristic equality tests [33, §5]. It can also happen that the easiest way to obtain a moderate-precision output involves highprecision intermediate computations, especially when the correctness of the result relies on pessimistic bounds.

Analytic continuation. Solutions of the differential equation (1) defined in the neighborhood of 0 extend by analytic continuation to the universal covering of \S, where S is the (finite) set of singularities of the equation. D-finite functions may be evaluated fast at any point by a numerical version of the analytic continuation process that builds on the previous algorithm [START_REF] Chudnovsky | Computer algebra in the service of mathematical physics and number theory[END_REF]. Rewrite [START_REF] Abad | An algorithm to obtain global solutions of the double confluent Heun equation[END_REF] in matrix form

Y ′ (z) = A(z)Y (z), Y (z) = y (i) (z) i! 0≤i<r . ( 6
)
This choice of Y (z) induces, for all z0 ∈ \ S, that of a family of canonical solutions of (1) defined by

y[z0, j](z) = (z -z0) j + O (z -z0) r , 0 ≤ j < r,
that form a basis of the solutions of (1) in the neighborhood of z0. Stated otherwise, the vector y[z0] = (y[z0, j]) 0≤j<r of canonical solutions at z0 is the first row of the fundamental matrix Y [z0](z) of ( 6) such that Y [z0](z0) = Id. By linearity, for any path z0 → z1 in \ S, there exists a transition matrix Mz 0 →z 1 that "transports initial conditions" (and canonical solutions) along the path:

Y (z1) = Mz 0 →z 1 Y (z0), y[z0](z) = y[z1](z)Mz 0 →z 1 . (7)
This matrix is easy to write out explicitly:

Mz 0 →z 1 = Y [z0](z1) = 1 i! y[z0, j] (i) (z1) 0≤i,j<r , (8)
evaluations at z1 being understood to refer to the analytic continuation path z0 → z1. Transition matrices compose:

Mz 0 →z 1 →z 2 = Mz 1 →z 2 Mz 0 →z 1 , Mz 1 →z 0 = M -1 z 0 →z 1 . ( 9
) NumGfun provides functions to compute Mγ for a given path γ (transition_matrix), and to evaluate the analytic continuation of a D-finite function along a path starting at 0 (evaldiffeq). In both cases, the path provided by the user is first subdivided into a new path z0 → z1 → • • • → zm, z ℓ ∈ É(i), such that, for all ℓ, the point z ℓ+1 lies within the disk of convergence of the Taylor expansions at z ℓ of all solutions of (1). Using [START_REF] Brent | Modern Computer Arithmetic[END_REF], approximations Mℓ ∈ É(i) r×r of Mz ℓ →z ℓ+1 are computed by binary splitting.

More precisely, we compute all entries of Mℓ at once, as follows. For a generic solution y of (1), the coefficients un of u(z) = y(z ℓ + z) = ∞ n=0 unz n are canceled by a recurrence of order s. Hence the partial sums u;n(z) of u satisfy

Un+1 u;n+1(z) = C(n)z 0 K 1 B(n) Un u;n(z) , ( 10
)
where K = (1, 0, . . . , 0) and C(n) ∈ É(n) s×s . Introducing an indeterminate δ, we let z = z ℓ+1 -z ℓ + δ ∈ É(i)[δ] and compute B(N -1) . . . B(0) mod δ r by binary splitting (an idea already used in [START_REF] Van Der Hoeven | Fast evaluation of holonomic functions[END_REF]), for some suitable N . The upper left blocks of the subproducts are kept factored as a power of z times a matrix independent on z. In other words, clearing denominators, the computation of each subproduct

P = 1 d D • p 0 R d = P high P low (p = numer(z) m ) is performed as D ← D high D low ; p ← p high p low ; R ← p low (R high C low ) + d high R low ; d ← d high d low .
The powers of z can be computed faster, but the saving is negligible.

Applying the row submatrix R of the full product to the matrix U0 = ( 1 i! y[z ℓ , j] (i) (z ℓ )) 0 i<s,0 j<r yields a row vector equal to (y[z ℓ , j];N (z ℓ+1 + δ)) 0 j<r + O(δ r ), each entry of which is a truncated power series whose coefficients are the entries of the corresponding column of Mℓ . This way of computing Mℓ is roughly min(r, s ω-1 ) times more efficient than the fastest of the variants from [START_REF] Chudnovsky | Computer algebra in the service of mathematical physics and number theory[END_REF][START_REF] Van Der Hoeven | Fast evaluation of holonomic functions[END_REF].

In the function transition_matrix, we then form the product Mm-1 • • • M0, again by binary splitting. In the case of evaldiffeq, we compute only the first row R of Mm-1, and we form the product R Mm-2 • • • M0 Ĩ, where Ĩ is an approximation with coefficients in É(i) of the vector Y (z0) of initial conditions (or this vector itself, if it has symbolic entries). The whole computation is done on unsimplified fractions, controlling all errors to guarantee the final result. Let us stress that no numerical instability occurs since all numerical operations are performed on rational numbers. We postpone the discussion of approximation errors (including the choice of N ) to §5.

Example 4. Transition matrices corresponding to paths that "go round" exactly one singularity once are known as local monodromy matrices. A simple example is that of the equation (1 + z 2 )y ′′ + 2zy ′ = 0, whose solution space is generated by 1 and arctan z. Rather unsurprisingly, around i:

> transition_matrix((1+z^2)*diff(y(z),z,y(z),z)
+2*z*diff(y(z),z), y(z), [0,I+1,2*I,I-1,0], 20);

1.00000000000000000000 3.14159265358979323846 0 1.00000000000000000000

More generally, expressing them as entries of monodromy matrices is a way to compute many mathematical constants.

Another use of analytic continuation appears in Ex. 2: there, despite the evaluation point lying within the disk of convergence, NumGfun performs analytic continuation along the path 0 → -1 2 → -3 4 → -22 25 → -47 50 → -99 100 to approach the singularity more easily (cf. [START_REF] Chudnovsky | Computer algebra in the service of mathematical physics and number theory[END_REF]Prop. 3.3]).

The "bit burst" algorithm. The complexity result (5) is quasi-optimal for h = O(log p), but becomes quadratic in p for h = Θ(p), which is the size of the approximation z ∈ É(i) needed to compute y(z) for arbitrary z ∈ . This issue can be avoided using analytic continuation to approach z along a path z0 → z1 → • • • → zm = z made of approximations of z whose precision grow geometrically:

|z ℓ+1 -z ℓ | ≤ 2 -2 ℓ , height(z ℓ ) = O(2 ℓ ), (11)
thus balancing h and |z|. This idea is due to Chudnovsky and Chudnovsky [START_REF] Chudnovsky | Computer algebra in the service of mathematical physics and number theory[END_REF], who called it the bit burst algorithm, and independently to van der Hoeven with a tighter complexity analysis [START_REF] Van Der Hoeven | Fast evaluation of holonomic functions[END_REF][START_REF] Van Der Hoeven | Transséries et analyse complexe effective[END_REF]. The following proposition improves this analysis by a factor log log p in the case where y has a finite radius of convergence. No similar improvement seems to apply to entire functions.

y(0) = 1 24 , y ′ (0) = 1 12 , y ′′ (0) = 5 24 , y ′′′ (0) = 5 24 .
The singular points are z ≈ 3.62 and z ≈ 0.09 ± 0.73i. By analytic continuation (along a path 0 → πi homotopic to a segment) followed by bit-burst evaluation, we obtain y(πi) ≈ -0.52299...(990 digits)...53279 -1.50272...90608i after about 5 min. This example is among the "most general" that NumGfun can handle. Intermediate computations involve recurrences of order 8 and degree 17.

REGULAR SINGULAR POINTS

The algorithms of the previous section extend to the case where the analytic continuation path passes through regular singular points of the differential equation [START_REF] Chudnovsky | On expansion of algebraic functions in power and Puiseux series (I, II)[END_REF][START_REF] Chudnovsky | Computer algebra in the service of mathematical physics and number theory[END_REF][START_REF] Van Der Hoeven | Fast evaluation of holonomic functions near and in regular singularities[END_REF]. Work is in progress to support this in NumGfun, with two main applications in view, namely special functions (such as Bessel functions) defined using their local behaviour at a regular singular point, and "connection coefficients" arising in analytic combinatorics [15, § VII.9]. We limit ourselves to a sketchy (albeit technical) discussion.

Recall that a finite singular point z0 of a linear differential equation with analytic coefficients is called regular when all solutions y(z) have moderate growth y(z) = 1/(z-z0) O (1) as z → z0 in a sector with vertex at z0, or equivalently when the equation satisfies a formal property called Fuchs' criterion. The solutions in the neighborhood of z0 then have a simple structure: for some t ∈ AE, there exist linearly independent formal solutions of the form (with z0 = 0)

y(z) = z λ t-1 k=0 φ k (z) log k z, λ ∈ , φ k ∈ [[z]] ( 12
)
in number equal to the order r of the differential equation. Additionally, the φ k converge on the open disk centered at 0 extending to the nearest singular point, so that the solutions [START_REF] Chudnovsky | Computer algebra in the service of mathematical physics and number theory[END_REF] in fact form a basis of analytic solutions on any open sector with vertex at the origin contained in this disk.

(See for instance [START_REF] Poole | Introduction to the theory of linear differential equations[END_REF] for proofs and references.) Several extensions of the method of indeterminate coefficients used to obtain power series solutions at ordinary points allow to determine the coefficients of the series φ k . We proceed to revisit Poole's variant [28, §16] of Heffter's method [START_REF] Heffter | Einleitung in die Theorie der linearen Differentialgleichungen[END_REF]Kap. 8] using the "operator algebra point of view" on indeterminate coefficients: a recursive formula for the coefficients of the series expansion of a solution is obtained by applying simple rewriting rules to the equation. This formulation makes the algorithm simpler for our purposes (compare [START_REF] Tournier | Solutions formelles d'équations différentielles[END_REF][START_REF] Chudnovsky | On expansion of algebraic functions in power and Puiseux series (I, II)[END_REF][START_REF] Van Der Hoeven | Fast evaluation of holonomic functions near and in regular singularities[END_REF]) and leads to a small complexity improvement. It also proves helpful in error control ( §6).

Since our interest lies in the D-finite case, we assume that 0 is a regular singular point of Equation [START_REF] Abad | An algorithm to obtain global solutions of the double confluent Heun equation[END_REF]. Letting θ = z d dz , the equation rewrites as L(z, θ) • y = 0 where L is a polynomial in two noncommutative indeterminates (and L(z, θ) has no nontrivial left divisor in K[z]). Based on [START_REF] Chudnovsky | Computer algebra in the service of mathematical physics and number theory[END_REF], we seek solutions as formal logarithmic sums

y(z) = n∈λ+ k≥0 y n,k log k z k! z n , λ ∈ .
(z) = n∈λ+ k≥0 y n-1,k log k z k! z n and θ • y(z) = n∈λ+ k≥0 (ny n,k + y n,k+1 ) log k z k! z n . Thus the coefficient sequence of L(z, θ) • y is L(S -1 n , n + S k ) • y.
Assume that y(z) satisfies (1). Then L(S -1 n , n+S k )•y = 0; additionally, [START_REF] Chudnovsky | Computer algebra in the service of mathematical physics and number theory[END_REF]

implies that y n,k = 0 for k ≥ t, which translates into S t k • y = 0. Set L(z, θ) = Q0(θ) + R(z, θ)z, and fix n ∈ λ + . If Q0(n) = 0, then Q0(n + X) is invertible in K(λ)[[X]],
and

y = -Q0(n + S k ) -1 mod S t k R(S -1 n , n + S k )S -1 n • y.
In general, let µ(n) be the multiplicity of n as a root of Q0.

Take Tn ∈ K(λ)[X] such that Tn(X) X -µ(n) Q0(n + X) = 1 + O(X t ) (explicitly, Tn(S k ) = t-1 v=0 ∂ v ∂X v X µ(n) Q 0 (X) X=n S v k ). Then, it holds that S µ(n) k • y = -Tn(S k )R(S -1 n , n + S k )S -1 n • y, (13)
hence the y n ′ ,k with n ′ < n determine (y n,k ) k≥µ(n) , while the y n,k , k < µ(n) remain free.

Coupled with the condition y n,k = 0 for n-λ < 0 following from [START_REF] Chudnovsky | Computer algebra in the service of mathematical physics and number theory[END_REF], this implies that a solution y is characterized by

(y n,k ) (n,k)∈E , where E = {(n, k) | k < µ(n)}. As in §3, this choice of initial values induces that of canonical solutions (at 0) y[n, k] defined by y[n, k] n,k = 1 and y[n, k] n ′ ,k ′ = 0 for (n ′ , k ′ ) ∈ E \ {(n, k)}.
The notion of transition matrix extends, see [33, §4] for a detailed discussion.

Equation ( 13) is a "recurrence relation" that lends itself to binary splitting. The main difference with the setting of §2 is that the companion matrix of the "recurrence operator" now contains truncated power series, i.e.,

B ∈ K(n)[[S k ]]/(S t k ). The coefficients y[u, v]n = t-1 k=0 y[u, v] n,k log k z of canoni- cal solutions may be computed fast by forming the product B(n -1) • • • B(u) ∈ K(λ)[[S k ]]/(S t
k ) and applying it to the initial values Yu = (0, . . . , 0, log v z) T . Compared to the algorithm of [START_REF] Van Der Hoeven | Fast evaluation of holonomic functions near and in regular singularities[END_REF], multiplications of power series truncated to the order t replace multiplications of t × t submatrices, so that our method is slightly faster. As in §3, all entries of Mz 0 →z 1 may (and should) be computed at once.

ERROR CONTROL

Performing numerical analytic continuation rigorously requires a number of bounds that control the behaviour of the function under consideration. We now describe how error control is done in NumGfun. Some of the ideas used in this section appear scattered in [START_REF] Van Der Hoeven | Fast evaluation of holonomic functions[END_REF]- [START_REF] Van Der Hoeven | On effective analytic continuation[END_REF]. NumGfun relies almost exclusively on a priori bounds; see [27, §5.2] for pointers to alternative approaches, and [START_REF] Van Der Hoeven | On effective analytic continuation[END_REF] for further useful techniques, including how to refine rough a priori bounds into better a posteriori bounds.

We start by discussing the computation of majorant series for the canonical solutions of the differential equation. Then we describe how these are used to determine, on the one hand, at which precision each transition matrix should be computed for the final result to be correct, and on the other hand, where to truncate the series expansions to achieve this precision. Finally, we propose a way to limit the cost of computing bounds in "bit burst" phases.

Remark 2. In practice, numerical errors that happen while computing the error bounds themselves are not always controlled, due to limited support for interval arithmetic in Maple. However, we have taken some care to ensure that crucial steps rely on rigorous methods.

Majorant series. A formal power series

g ∈ Ê+[[z]] is a majorant series of y ∈ [[z]
], which we denote by y g, if |yn| ≤ gn for all n. If y(z) g(z) and ŷ(z) ĝ(z), then

y(z) ≤ g(|z|), d dz y(z) d dz g(z), y + ŷ g + ĝ, y ŷ gĝ, y • ŷ g • ĝ. (14)
We shall allow majorant series to be formal logarithmic sums or matrices. The relation extends in a natural way: we write n,k y n,k z n log z k n,k g n,k z n log z k iff |y n,k | ≤ g n,k for all n and k, and Y = (yi,j) G = (gi,j) iff Y and G have the same format and yi,j gi,j for all i, j. In particular, for scalar matrices, Y G if |yi,j| ≤ gi,j for all i, j. The inequalities [START_REF] Falloon | Theory and computation of spheroidal wavefunctions[END_REF] still hold.

Bounds on canonical solutions. The core of the error control is a function that computes g(z) such that ∀j, y[z0, j](z0 + z) g(z) [START_REF] Flajolet | Analytic Combinatorics[END_REF] (in the notations of §3) given Equation ( 1) and a point z0. This function is called at each step of analytic continuation.

The algorithm is based on that of [START_REF] Mezzarobba | Effective bounds for P-recursive sequences[END_REF] ("SB" in what follows), designed to compute "asymptotically tight" symbolic bounds. Run over É(i) instead of É, SB returns (in the case of convergent series) a power series satisfying [START_REF] Flajolet | Analytic Combinatorics[END_REF] n) . The series g is specified by τ , α, and other parameters of no interest here that define φ. The tightness property means that τ and α are the best possible.

of the form g(z) = ∞ n=0 n! -τ α n φ(n)z n , where τ ∈ É+, α ∈ É, and φ(n) = e o(
However, the setting of numerical evaluation differs from that of symbolic bounds in several respects: (i) the issue is no more to obtain human-readable formulae, but bounds that are easy to evaluate; (ii) bound computations should be fast; (iii) at regular singular points, the fundamental solutions are logarithmic sums, not power series. For space reasons, we only sketch the differences between SB and the variant we use for error control.

First, we take advantage of (i) to solve (ii). The most important change is that the parameter α is replaced by an approximation α ≥ α. This avoids computations with algebraic numbers, the bottleneck of SB. Strictly speaking, it also means that we are giving up the tightness of the bound. However, in constrast with the "plain" majorant series method [START_REF] Van Der Hoeven | Fast evaluation of holonomic functions near and in regular singularities[END_REF][START_REF] Van Der Hoeven | Majorants for formal power series[END_REF][START_REF] Van Der Hoeven | On effective analytic continuation[END_REF], we are free to take α arbitrarily close to α, since we do not use the ratio (α/ α) n to mask subexponential factors. The algorithms from [START_REF] Mezzarobba | Effective bounds for P-recursive sequences[END_REF] adapt without difficulty. Specifically, Algorithms 3 and 4 are modified to work with α instead of α. Equality tests between "dominant roots" (Line 9 of Algorithm 3, Line 5 of Algorithm 4) can now be done numerically and are hence much less expensive. As a consequence, the parameter T from §3.3 is also replaced by an upper bound. This affects only the parameter φ(n) of the bound, which is already pessimistic.

The issue (iii) can be dealt with too. One may compute a series g such that y(z) g(z)

t-1 k=0

log k z k!
(with the notations of §4) by modifying the choice of K in [27, §3.3] so that

t-1 k=0 ∂ k ∂X k X µ(n)
Q 0 (X) X=n ≤ K/n r , and replacing [27, Eq. ( 14)] by Eq. ( 13) above in the reasoning.

Global error control. We now consider the issue of controlling how the approximations at each step of analytic continuation accumulate. Starting with a user-specified error bound ǫ, we first set ǫ ′ so that an ǫ ′ -approximation r ∈ É(i) of the exact result r rounds to ⋄(r) ∈ n∈AE 10 -n [i] with |⋄(r) -r| < ǫ. No other rounding error occur, since the rest of the computation is done on objects with coefficients in É(i). However, we have to choose the precision at which to compute each factor of the product Π = R Mm-1 • • • M0 Ĩ (in evaldiffeq, and of similar products for other analytic continuation functions) so that |r -r| < ǫ ′ .

As is usual for this kind of applications, we use the Frobenius norm, defined for any (not necessarily square) matrix by (ai,j) F = ( i,j |ai,j| 2 ) 1/2 . The Frobenius norm satisfies AB F ≤ A F B F for any matrices A, B of compatible dimensions. If A is a square r × r matrix,

A ∞ ≤ A 2 ≤ A F ≤ r A ∞ ( 16
)
where • 2 is the matrix norm induced by the Euclidean norm while • ∞ denotes the entrywise uniform norm. Most importantly, the computation of A F is numerically stable, and if A has entries in É(i), it is easy to compute a good upper bound on A F in floating-point arithmetic.

We bound the total error ǫ Π on Π using repeatedly the rule

à B -AB F ≤ à F B -B F + Ã-A F B F .
From there, we compute precisions ǫA such that having à -A F < ǫA for each factor A of Π ensures that ǫ Π < ǫ ′ . Upper bounds on the norms à F and A F appearing in the error estimate are computed either from an approximate value of A (usually à itself) if one is available at the time, or from a matrix G such that A G given by [START_REF] Falloon | Theory and computation of spheroidal wavefunctions[END_REF][START_REF] Flajolet | Analytic Combinatorics[END_REF].

Local error control.

Let us turn to the computation of individual transition matrices. We restrict to the case of ordinary points. Given a precision ǫ determined by the "global error control" step, our task is to choose N such that M -Mz 0 →z 1 F ≤ ǫ if M is computed by truncating the entries of (8) to the order N . If g satisfies [START_REF] Flajolet | Analytic Combinatorics[END_REF], it suffices that g

(i) N; (|z1 -z0|) ≤ i!
r ǫ for all i, as can be seen from [START_REF] Falloon | Theory and computation of spheroidal wavefunctions[END_REF][START_REF] Fousse | MPFR: A multiple-precision binary floating-point library with correct rounding[END_REF]. We find a suitable N by dichotomic search. For the family of majorant series g used in NumGfun, the g (i)

N; (x) are not always easy to evaluate, so we bound them by expressions involving only elementary functions [27, §4.2]. The basic idea, related to the saddle-point method from asymptotic analysis, is that if x ≥ 0, inequalities like gn;(x) ≤ (x/tn) n g(tn)(1 -x/tn) -1 are asymptotically tight for suitably chosen tn.

Points of large bit-size. Computing the majorants ( 15) is expensive when the point z0 has large height. This can be fixed by working with an approximate value c of z0 to obtain a bound valid in a neighborhood of c that contains z0. This technique is useful mainly in "bit burst" phases (where, additionally, we can reuse the same c at each step).

Assume that

Y [c](c + z) G(z) for some c ≈ z0. Since Y [z0](z0 + z) = Y [c](c + (z0 -c) + z)M -1 c→z 0 by (7), it follows that Y [z0](z0 + z) G(|z0 -c| + z) M -1 c→z 0 1 where 1 is a square matrix filled with ones. Now Mc→c+z = Id +O(z), whence Mc→c+z -Id G(z) -G(0). For |z0 -c| < η, this implies Mc→z 0 -Id F ≤ δ := G(η) -G(0) F . Choosing η small enough that δ < 1, we get M -1 c→z 0 F ≤ (1 -δ) -1 . If G was computed from (15), i.e., for G(z) = 1 j! g (j) (z) i,j , this finally gives the bound y[z0, j](z0 + z) r j! (1 -δ) g (j) (η + z),
valid for all z0 in the disk |z0 -c| < η. Note however that simple differential equations have solutions like exp(K/(1 -z)) with large K. The condition δ < 1 then forces us to take c of size Ω(K). Our strategy to deal with this issue in NumGfun is to estimate K using a point c of size O(1) and then to choose a more precise c (as a last resort, c = z0) based on this value if necessary. Remark 3. If the evaluation point z is given as a program, a similar reasoning allows to choose automatically an approximation z ∈ É(i) of z such that |y(z) -y(z)| is below a given error bound [32, §4.3]. In other applications, it is useful to have bounds on transition matrices that hold uniformly for all small enough steps in a given domain. Such bounds may be computed from a majorant differential equation with constant coefficients [35, §5].

REPEATED EVALUATIONS

Drawing plots or computing integrals numerically requires to evaluate the same function at many points, often to a few digits of precision only. NumGfun provides limited support for this through the function diffeqtoproc, which takes as input a D-finite function y, a target precision ǫ, and a list of disks, and returns a Maple procedure that performs the numerical evaluation of y. For each disk D = {z, |z -c| < ρ}, diffeqtoproc computes a polynomial p ∈ É(i) [z] such that |⋄(p(z)) -y(z)| < ǫ for z ∈ D ∩ É(i), where ⋄ again denotes rounding to complex decimal. The procedure returned by diffeqtoproc uses the precomputed p when possible, and falls back on calling evaldiffeq otherwise.

The approximation polynomial p is constructed as a linear combination of truncated Taylor series of fundamental solutions y[c, j], with coefficients obtained by numerical analytic continuation from 0 to c. The way we choose expansion orders is similar to the error control techniques of §5: we first compute a bound Bcan on the fundamental solutions and their first derivatives on the disk D. The vector Y (c) of "initial values" at c is computed to the precision ǫ ′ /Bcan where ǫ ′ ≤ ǫ/(2r). We also compute Bini ≥ Y (c) F . Each y[c, j] is expanded to an order N such that y[c, j]N; ∞,D ≤ ǫ ′ /Bini, so that finally |p(z) -y(z)| ≤ ǫ for z ∈ D.

The most important feature of diffeqtoproc is that it produces certified results. At low precisions and in the absence of singularities, we expect that interval-based numerical solvers will perform better while still providing (a posteriori) guarantees. Also note that our simple choice of p is far from optimal. If approximations of smaller degree or height are required, a natural approach is to aim for a slightly smaller error y -p ∞,D above, and then replace p by a polynomial p for which we can bound pp ∞ [36, §6.2].

Example 6. The following plot of the function y defined by (z-1) y ′′′ -z(2z-5) y ′′ -(4z-6) y ′ +z 2 (z-1) y with the initial values y(0) = 2, y ′ (0) = 1, y ′′ (0) = 0 was obtained using polynomial approximations on several disks that cover the domain of the plot while avoiding the pole z = 1. The whole computation takes about 9 s. Simple numerical integrators typically fail to evaluate y beyond z = 1.

FINAL REMARKS

Not all of NumGfun was described in this article. The symbolic bounds mentioned in §5 are also implemented, with functions that compute majorant series or other kinds of bounds on rational functions (bound_ratpoly), D-finite functions (bound_diffeq and bound_diffeq_tail) and Precursive sequences (bound_rec and bound_rec_tail). This implementation was already presented in [START_REF] Mezzarobba | Effective bounds for P-recursive sequences[END_REF].

Current work focuses on adding support for evaluation "at regular singular points" (as outlined in §4), and improving performance. The development version of NumGfun already contains a second implementation of binary splitting, written in C and called from the Maple code. In the longer term, I plan to rewrite other parts of the package "from the bottom up", both for efficiency reasons and to make useful subroutines independant of Maple.

Proposition 3 .

 3 Let us call the double sequence y = (y n,k ) n∈λ+ ,k∈AE the coefficient sequence of y. The shift operators Sn and S k on such double sequences are defined by Sn • y = (y n+1,k ), and S k • y = (y n,k+1 ). The heart of the method lies in the following observation. Let y(z) be a formal logarithmic sum. The relation L(z, θ) • y = 0 holds (formally) iff the coefficient sequence y of y satisfies L(S -1 n , n + S k ) • y = 0. Proof. The operators z and θ act on logarithmic sums by z • y

  it remains to reduce to common denominator the whole equation; hence d ′ ≤ sd and h

′ ≤ s(h ′′ + log s + log d). These two conversion steps take O(M (sdh log 2 d)) and O(M (d ′ h ′ log s)) operations respectively, using product trees. The assumption d, s = o(N ) allows to write H = O(sd(h + log N )

  /6 in the FFT model, and C = (3 MM(s) + 1)/4 in the black-box model. Here MM(s) ≤ (s 3 + 3s 2 )/2 is the algebraic complexity of s × s matrix multiplication over .

  Proof. By[START_REF] Bostan | Fast algorithms for polynomial solutions of linear differential equations[END_REF] and[START_REF] Chudnovsky | On expansion of algebraic functions in power and Puiseux series (I, II)[END_REF], the step z ℓ → z ℓ+1 of the bitburst algorithm takes O(M (p(2 ℓ + log p) log p/2 ℓ )) bit operations. Now

	m ℓ=0	M		p(2 ℓ + log p) 2 ℓ	log p ≤ M mn log p +	∞ ℓ=0	p log 2 p 2 ℓ
	and m = O(log p), hence the total complexity.
	Example 5. Consider the D-finite function y defined by the following equation, picked at random:
	5 12	-	1 4	z +	19 24	z 2 -	5 24	z 3	y (4) + -	7 24	+	2 3	z +	13 24	z 2 +	1 12	z 3	y ′′′
		+	7 12	-	19 24	z +	1 8	z 2 +	1 3	z 3	y ′′ + -	3 4	+	5 12	z +	5 6	z 2 +	1 2	z 3	y ′
														+	5 24	+	23 24	z +	7 8	z 2 +	1 3	z 3	y = 0,

Proposition 2. Let y be a D-finite function. One may compute a 2 -p -approximation of y at a point z ∈ É(i) of height O(p) in O(M (p log 2 p)) bit operations.

All timings reported in this article were obtained with the following configuration: Intel T7250 CPU at

GHz, 1 GB of RAM, Linux 2.6.32, Maple 13, GMP 4.2.1.

Acknowledgements. I am grateful to my advisor, B. Salvy, for encouraging me to conduct this work and offering many useful comments. Thanks also to A. Benoit, F. Chyzak, P. Giorgi, J. van der Hoeven and A. Vaugon for stimulating conversations, bug reports, and/or comments on drafts of this article, and to the anonymous referees for helping make this paper more readable.

This research was supported in part by the MSR-INRIA joint research center.


Youssef Elfoutayeni
email: youssef.elfoutayeni@campusmarrakech.com

Mohamed Khaladi
email: khaladi@ucam.ac.ma

Using vector divisions in solving linear complementarity problem

Keywords: Linear complementarity problem, Vector division, Global convergence, Newton's method, secant method

, he uses vector divisions with the secant method; but for using this method we must have a function continuous with partial derivatives on an open set of IR n ; so we built a sequence of functions F (p, x) ∈ C ∞ which converges uniformly to the function F (x); and we show that finding the zero of the function F is completely equivalent to finding the zero of the sequence of the functions F (p, x). We close our paper with some numerical simulation examples to illustrate our theoretical results.

Introduction:

The complementarity problem noted (CP ) is a classical problem of optimization theory of finding (z, w) ∈ IR n × IR n such that:

   < z, w >= 0 w -f (z) = 0 z, w 0 (1)
where f , a continuous operator from IR n into itself, is given data.

The constraint < z, w >= 0 is called the complementarity condition since for any i, 1 i n, z i = 0 if w i > 0, and vice versa; it may be the case that z i = w i = 0.

In the case that the function f is a nonlinear continuous operator from IR n into itself, so the problem is called a NonLinear Complementarity Problem associated with the function f and noted (N LCP ).

In the case that the function f is affine, i.e f (z) = q + M z, where q is an element of IR n and M is a real square matrix of order n, so the problem is called a Linear Complementarity Problem associated with the matrix M and the vector q and noted (LCP ).

Solving the (LCP ) in general, however, appears to be difficult. One simple reason is that there is no known good characterization of the nonexistence of a solution to the system for any given function f . The linear complementarity problem plays an important role in several fields (game theory, operational research ...); moreover, Cottle et Dantzig [START_REF] Cottle | A life in mathematical programming[END_REF] et Lemke [START_REF] Lemke | Bimatrix equilibrium points and mathematical programming[END_REF] have proved that all problems of linear programming (LP ), convex quadratic programming (CQP ), and also the problems of Nash equilibrium of a game bi-matrix can be written as a linear complementarity problem.

The question is precisely under what conditions on the matrix M and the vector q this problem admits one and only one solution, if this is the case, how can we express this solution as a function of the matrix and vector mentioned above. This question has not been completely resolved yet. However, many results already exist, for instance Lemke [START_REF] Lemke | Bimatrix equilibrium points and mathematical programming[END_REF] who gave sufficient conditions on the matrix M and the vector q under which the number of solutions of LCP (M, q) is finite. Samelson [START_REF] Samelson | A partition theorem for Euclidean n-space[END_REF], Ingeton [START_REF] Ingleton | Aproblem in Linear Inequalites[END_REF], Murty [START_REF] Murty | On a characterization of P-matrices[END_REF], Watson [START_REF] Watson | A Variational Approach to the Linear Complementaruty Problem[END_REF], Kelly [START_REF] Kelly | Q-Matrices and Spherical Geometry[END_REF] and Cottle [START_REF] Cottle | The Lineair Complementarity Problem[END_REF] have by contrast shown that the matrix M is a P -matrix if and only if the linear complementarity problem associated with a matrix M and a vector q has a unique solution for all q ∈ IR n (We should remind that a matrix M is called a P -matrix if all principal minors are strictly positive (see [START_REF] Fiedler | On matrices with non-positive off-diagonal elements and positive principal minors[END_REF]), and we should note Y. ELFOUTAYENI and M. KHALADI 2/11 that any symmetric and positive definite matrix is a P -matrix, but not vice-versa).

Preliminaries:

In this section, we summarize some notations which will be used in this paper.

In particular, IR n denotes the space of real n-dimensional vectors,

IR n + := {x ∈ IR n : x i 0, i = 1.
.n} is the nonnegative orthant and its interior is IR n

++ := {x ∈ IR n : x i > 0, i = 1..n}. With x ∈ IR n we define |x| = (|x 1 |, .., |x n |) T ∈ IR n .
We denote by I the identity matrix. Let x, y ∈ IR n , x T y or < x, y > is the inner product of the x and y; ||x|| is the Euclidean norm.

For x ∈ IR n and k a nonnegative integer, x (k) refers to the vector obtained after k iterations; for 1 i n, x i refers to the i th element of x, and x (k) i refers to the i th element of the vector obtained after k iterations. Let x, y ∈ IR n , the expression x y (respectively x < y) meaning that x i y i (respectively x i < y i ) for each i = 1..n.

For x ∈ IR n we denote by e x = (e x1 , .., e xn ) T ∈ IR n and for x ∈ IR n ++ we denote by ln(x) = (ln(x 1 ), .., ln(x n )) T . The transpose of a vector is denoted by super script T , such as the transpose of the vector x is given by x T .

Remember that the spectrum σ(A) of the matrix A is the set of its eigenvalues and its spectral radius ρ is given by: ρ(A) := sup{|λ| such that λ ∈ σ(A)}.

Equivalent reformulation of LCP:

It is known in [START_REF] Murty | Linear Complementarity, Linear and Nonlinear Programming[END_REF] that the linear complementarity problem LCP (M, q) is completely equivalent to solving nonlinear equation

F (x) = 0 with F is a function from IR n into itself defined by F (x) := (M + I)x + (M -I)|x| + q;
More precisely (see [START_REF] Murty | Linear Complementarity, Linear and Nonlinear Programming[END_REF]), on the one hand, if x * is a zero of the function F , then

z * := |x * | + x * w * := |x * | -x * (2)
define a solution of LCP (M, q).

Y. ELFOUTAYENI and M. KHALADI 3/11

On the other hand, if (z * , w * ) is a solution of LCP (M, q), then

x * := z * -w * 2
is a zero of the function F .

We mention that this equation is solved by the fixed point algorithm (see [START_REF] Schafer | On the moduls algorithm for the linear complementarity problem[END_REF]), this algorithm is defined by:

x (0) ∈ IR n arbitrary, x (k+1) = (I + M ) -1 (I -M )|x (k) | -(I + M ) -1 q (3)
For the case that M is symetric and positive definite, it was shown in [START_REF] Van Bokhoven | Piecewise-linear modelling and analysis[END_REF] (see also Section 9.2 in [START_REF] Murty | Linear Complementarity, Linear and Nonlinear Programming[END_REF]) that

D := (I + M ) -1 (I -M ) it holds ||D|| 2 = ρ(D T D) = ρ(D 2 ) = ρ(D) < 1,
where ρ(.) denotes the spectral radius of a matrix; hence (3) converges by the contraction-mapping theorem (see Theorem 5.1.3 in [START_REF] Ortega | Iterative solution of nonlinear equations in serval variables[END_REF]) and

x * = lim k→+∞ x (k)
is the unique solution of the F (x) = 0. We mention also that the convergence of algorithm (3) is only linear; in this paper, we consider the use of vector divisions with the secant method (see [START_REF] Shi | Using vector divisions in solving nonlinear systems of equations[END_REF]) in instead of the algorithm (3) decribed above, this algorithm has global convergence (see [START_REF] Shi | Using vector divisions in solving nonlinear systems of equations[END_REF]); but for using this algorithm we must have a function F a continuous with partial derivatives on a set of IR n ; the next section can answer this problematic.

The main result

We consider again the function F : IR n → IR n defined by Remark 3 : The uniqueness of the root of the function F results from the uniqueness of the solution of linear complementarity problem LCP (M, q), in fact, supposing that x * 1 and x * 2 , two distinct roots of the function F ,exist, then

z * 1 := |x * 1 | + x * 1 . z * 2 := |x * 2 | + x * 2 .
Since z * 1 = z * 2 (uniqueness of the solution of LCP (M, q)) then

|x * 1 | + x * 1 = |x * 2 | + x * 2 . ( 4
)
Simularly, we use the same method for w * 1 and w * 2 we have

|x * 1 | -x * 1 = |x * 2 | -x * 2 .
(5) so (4)-( 5)

means that x * 1 = x * 2 .
Now, we give the following algorithm for solving F (p, x) = 0 (see [START_REF] Shi | Using vector divisions in solving nonlinear systems of equations[END_REF]):

Algorithm:

Step 0: Determine ǫ, p, k * , ρ, σ such that k * is a positive integer, 0 < ρ < 1/2, and ρ < σ < 1;

Step 1: Select two points x (0) and x (1) ∈ IR n ;

Step 2: for k = 1, 2, ... until termination, do the following:

1-Compute the steepset descent direction

d (k) := -J(p, x (k) ) T F (p, x (k) ),
δ ij = 0 if i = j. 2-if k equals a multiple of k * , then insert a steepset descent direc- tion step, that is, let s (k) := d (k)
and go to step 2.7; 3-Compute:

u (k) := ξ 1 F (p, x (k) ) v (k) := ξ 2 (x (k) -x (k-1) )
Y. ELFOUTAYENI and M. KHALADI 6/11

with

ξ 1 = - ||x (k) -x (k-1) || 2 < x (k) -x (k-1) , F (p, x (k) ) -F (p, x (k-1) ) > and ξ 2 = - < F (p, x (k) ) -F (p, x (k-1) ), F (p, x (k) ) > || F (p, x (k) ) -F (p, x (k-1) )|| 2 . 4-if (u (k) -v (k) ) T d (k) = 0, then choose α (k) > -< v (k) , d (k) > < u (k) -v (k) , d (k) > such that α (k) maximizes the value of Cos[s (k) , d (k) ] = < s (k) , d (k) > ||s (k) ||.||d (k) || ; Set s (k) := α (k) u (k) + (1 -α (k) )v (k)
and go to step 2.7;

5-if (u (k) -v (k) ) T d (k) = 0 and < v (k) , d (k) >> 0, then set s (k) := u (k) + v (k) 2
and go to step 2.7;

6-if (u (k) -v (k) ) T d (k) = 0 and < v (k) , d (k) > 0, then set s (k) := d (k)
and go to step 2.7; 7-take a line search along the direction s (k) to determine the step length γ suth that

f (p, x (k) + γs (k) ) f (p, x (k) ) -γρ < d (k) , s (k) > and < ▽f (p, x (k) + γs (k) ), s (k) > -σ < d (k) , s (k) >; where f (p, x) = 1 2 || F (p, x)|| 2 8-set x (k+1) := x (k) + γs (k)
and go to next iteration.

Numerical examples

In this part, we consider some examples to test our method. The results of Fixed Point and using vector divisions methods for these examples are presented here for comparison purposes. The results and expected solutions for each example have been presented on the Tables 1 and2.

Example 4 : Consider the following linear complementarity problem:

Find vector z in IR n satisfying z T (M z + q) = 0, M z + q 0, z 0,

where M =     4 -1 0 0 -1 4 -1 0 0 -1 4 -1 0 0 -1 4     and q =     -4 3 -4 2     .
The exact solutions is x * = (1, 0, 1, 0) T .

We apply the fixed point and using vector divisions methods to solve this example with the initial approximation

x (0) = (1.1, 0.1, 1.2, 0.2) T .
The solution of this problem with six significant digits is presented in Ta-ble1.

Iteration z 1 z 2 z 3 z 4 k=01
0,0000000 0,0000000 0,0000000 0,0000000 k=05 0,7883251 0,0000000 1,3448593 0,0000000 k=10 1,0197946 0,0000000 0,9884737 0,0000000 Fixed k=15 0,9985643 0,0000000 1,0012907 0,0000000 point k=20 1,0001030 0,0000000 0,9999377 0,0000000 method k=25 0,9999918 0,0000000 1,0000058 0,0000000 k=30 1,0000005 0,0000000 0,9999997 0,0000000 k=33 0,9999999 0,0000000 1,0000001 0,0000000 k=34 1,0000001 0,0000000 1,0000000 0,0000000 k=35 1,0000000 0,0000000 1,0000000 0,0000000 Using vector k=01 0,0000000 0,0000000 0,0000000 0,0000000 divisions k=02 4,0000000 0,0000000 4,0000000 0,0000000 method k=03 1,0000000 0,0000000 1,0000000 0,0000000 Table 1: The results of different methods for example1.

Example 5 : Let's solve the following linear complementarity problem Find vector z in IR n satisfying z T (M z + q) = 0, M z + q 0, z 0,

where M =     8 -1 0 -5 1 5 -1 0 2 -1 6 -1 6 0 -1 7     and q =     1 -2 -3 4     .
The exact solutions is x * = (0, 15 29 , 17 29 , 0) T . We apply the fixed point and using vector divisions methods to solve this example with the initial approximation x (0) = (-1, -2, -3, -4) T .

The solution of this problem with six significant digits is presented in Ta-Y. ELFOUTAYENI and M. KHALADI ble2.

Iteration

z 1 z 2 z 3 z 4 k=01 0,0000000 0,0000000 0,0000000 0,0000000 k=05 0,1115059 0,6300237 0,3146258 0,0000000 k=10 0,0000000 0,4906849 0,6581657 0,0000000 k=15 0,0000000 0,5309163 0,5863067 0,0000000 k=20 0,0000000 0,5192743 0,5936312 0,0000000 k=25 0,0000000 0,5188339 0,5869623 0,0000000 k=30 0,0000000 0,5173795 0,5862876 0,0000000 Fixed k=35 0,0000000 0,5171712 0,5859293 0,0000000 point k=40 0,0000000 0,5171411 0,5860432 0,0000000 method k=45 0,0000000 0,5171911 0,5861363 0,0000000 k=50 0,0000000 0,5172239 0,5861902 0,0000000 k=55 0,0000000 0,5172388 0,5862079 0,0000000 k=60 0,0000000 0,5172428 0,5862107 0,0000000 k=65 0,0000000 0,5172428 0,5860930 0,0000000 k=70 0,0000000 0,5172421 0,5862079 0,0000000 k=75 0,0000000 0,5172416 0,5862071 0,0000000 k=79 0,0000000 0,5172424 0,5862069 0,0000000 k=01 0,0000000 2,0000000 3,0000000 0,0000000 k=02 0,0000000 0,0000000 0,4367816 0,0000000 Using vector k=03 0,0000000 2,6264367 0,5000000 0,0000000 divisions k=04 0,0000000 0,6983749 0,6163958 0,0000000 method k=05 0,0000000 0,5172414 0,5862069 0,0000000 k=06 0,0000000 0,5172414 0,5862069 0,0000000

Conclusion:

In this paper we have used that the linear complementarity problem is completely equivalent to solving nonlinear equation F (x) = 0 with F is a function from IR n into itself defined by F (x) = (M + I)x + (M -I)|x| + q; for solving this equation we have used the Shi's method, this method uses vector divisions with the secant method; based on that, we have proposed a globally convergent hybrid algorithm for solving this equation; to do so, we had to build a sequence of functions F (p, x) ∈ C ∞ which converges uniformly to the function F (x); and we have shown that finding the zero of the function F is completely equivalent to finding the zero of the sequence of the functions F .

  Therefore, w * := |x * | -x * and z * := |x * | + x * define the unique solution of the LCP (M, q).

FF

  (x) := (M + I)x + (M -I)|x| + q and let's consider the sequence of functions F : IN * ×IR n → IR n defined by F (p, x) := (M + I)x + 1 p (M -I) ln(e 0 + e px + e -px ) + q Proposition 1 : F (p, x) converges uniformly to F (x) when p → +∞. Proof. : To show that we can start by 1 p ln(e 0 + e px + e -px ) -|x| = 1 p [ln(e 0 + e px + e -px ) + ln(e -p|x| )] = 1 p ln[e -p|x| * (1 + e px + e -px )] = 1 p ln(e -p|x| + e p(x-|x|) + e -p(x+|x|) ) then we have 0 1 p ln(e 0 + e px + e -px ) -|x| 1 p ln(3) so 1 p ln(e 0 + e px + e -px ) is uniform convergence to |x| when p → +∞; Moreover, the operator (M -I) is linear then we have 1 p (M -I) ln(e 0 +e px +e -px ) converges uniformly to (M -I)|x| as p → +∞ and from the expression of the sequence of the functions F and the function F we have F (p, x) converges uniformly to F (x) when p → +∞. Theorem 2 : Let x * (p) be a solution of the equation F (p, x) = 0, then x * (p) is an approximation solution of F (x) = 0 for p is large enough. Proof. :To show that, we use proposition(1) wich we can interpret as ∀ǫ > 0, ∃ p * > 0 such that for all p > p * we have ||F (x * (p))|| = ||F (x * (p)) -F (p, x * (p))|| ǫ then we have for any ǫ > 0, x * (p) is the approximation solution of

whereF

  (p, x) := (M + I)x + 1 p (M -I) ln(e 0 + e px + e -px ) + q; J(p, x) := (M + I) + (M -I)E(p, x) and E ij (p, x) := δ ij e pxi -e -pxi 1 + e pxi + e -pxi ; we racall that δ ii = 1 and

Table 2 :

 2 The results of different methods for example2.

Y. ELFOUTAYENI and M. KHALADI