Marc Mezzarobba
email: marc@mezzarobba.net

A NOTE ON THE SPACE COMPLEXITY OF FAST D-FINITE FUNCTION EVALUATION

) ) and space O(p). The standard fast algorithm for this task, due to Chudnovsky and Chudnovsky, achieves the same time complexity bound but requires Θ(p lg p) bits of memory.

This work is in the public domain. As such, it is not subject to copyright. Where it is not legally possible to consider this work as released into the public domain, any entity is

Introduction

Binary splitting is a well-known and widely applicable technique for the fast multiple precision numerical evaluation of rational series. For any series n s n with lim sup n |s n | 1/n < 1 whose terms s n obey a linear recurrence relation with polynomial coefficients, e.g., ln 2 = ∞ n=0 s n , s n = 1 (n + 1)2 n+1 , 2(n + 2)s n+1 -(n + 1)s n = 0, the binary splitting algorithm allows one to compute the partial sum N -1 n=0 s n in O(M(N (lg N ) 2 )) bit operations [START_REF] Chudnovsky | Computer algebra in the service of mathematical physics and number theory[END_REF][START_REF] Brent | Modern Computer Arithmetic[END_REF]. Here M(n) stands for the complexity of multiple precision integer multiplication, and lg denotes the binary logarithm. As N = O(p) terms of the series are enough to make the approximation error less than 2 -p , the complexity of the algorithm is softly linear in the precision p, assuming M(n) = O(n(lg n) O(1) ).

Methods based on binary splitting tend to be favored in practice even in cases when asymptotically faster algorithms (typically AGM iterations [START_REF] Borwein | Pi and the AGM[END_REF]) would apply. One high-profile example is the computation of billions of digits of classical constants such as π, ζ [START_REF] Brent | Modern Computer Arithmetic[END_REF] or γ. Basically all record computation in recent years were achieved by evaluating suitable series using variants of binary splitting [START_REF] Gourdon | Constants and records of computation[END_REF][START_REF] Yee | Mathematical constants -billions of digits[END_REF].

A drawback of the classical binary splitting algorithm, both from the complexity point of view and in practice, is its comparatively large memory usage. Indeed, the algorithm amounts to the computation of a product tree of matrices derived from the recurrence-see Sect. 3 below for details. The intermediate results are matrices of rational numbers whose bit sizes roughly double from one level to the next. Near the root, their sizes can (and in general do) reach Θ(p lg p), even though the output has size Θ(p).

However, the space complexity can be lowered to O(p) using a slight variation of the classical algorithm. The basic idea is to truncate the intermediate results to a precision O(p) when they start taking up more space than the final result.

Of course, these truncations introduce errors. To make the trick into a genuine algorithm, we need to analyze the errors, add a suitable number of "guard digits" at each step and check that the space and time complexities of the resulting process stay within the expected bounds.

The opportunity to improve the practical behavior of binary splitting using truncations has been noticed by authors of implementations on several occasions over the last decade or so. Gourdon and Sebah [START_REF] Gourdon | Binary splitting method[END_REF] describe truncation as a "crucial" optimization. Besides the expected drop of memory usage, they report running time improvements by an "appreciable" constant factor. Cheng et al. [START_REF] Cheng | Time-and space-efficient evaluation of some hypergeometric constants[END_REF] compare truncation with alternative (less widely applicable but sometimes more efficient) approaches. Most recently, Kreckel [START_REF] Kreckel | decimal(γ)[END_REF] explicitly asks how to make sure that the new roundoff errors do not affect the correctness of the result. Indeed, the above-mentioned error analysis did not appear in the literature until very recently. An article by Yakhontov [START_REF] Yakhontov | Calculation of hypergeometric series with quasi-linear time and linear space complexity[END_REF]27] now provides the required bounds in the case of the generalized hypergeometric series p F q , which covers all examples where the truncation trick had been used before. But the applicability of the method is actually much wider.

The purpose of this note is to present a more general and arguably simpler analysis. Our version is more general in two main respects. First, besides hypergeometric series, it applies to the solutions of linear ordinary differential equations with rational coefficients, also known as D-finite (or holonomic) series [START_REF] Stanley | Differentiably finite power series[END_REF]. D-finite series are exactly those whose coefficients obey a linear recurrence relation with rational coefficients, while hypergeometric series correspond to recurrences of the first order. Second, we take into account the coefficient size of the recurrence that generates the series to be computed. Allowing the size of the coefficients to vary with the target precision p makes it possible to use the modified binary splitting procedure as part of the "bit burst" algorithm [START_REF] Chudnovsky | Computer algebra in the service of mathematical physics and number theory[END_REF] to handle evaluations at general real or complex points approximated by rationals of size Θ(p).

Additionally, our analysis readily adapts to other applications of binary splitting. The simplicity and generality of the proof are direct consequences of viewing the algorithm primarily as the computation of a product tree. See Gosper [START_REF] Gosper | Strip mining in the abandoned orefields of nineteenth century mathematics[END_REF] and Bernstein [1, §12-16] for further comments on this point of view.

The remainder of this note is organized as follows. Section 2 contains some notations and assumptions. In Sect. 3, we recall the standard binary splitting algorithm, which will serve as a subroutine in the linear-space version. Then, in Sect. 4, we state and analyze the "truncated" variant that achieves the linear space complexity for general D-finite functions. Finally, Sect. 5 offers a few comments on other variants of the binary splitting method and possible extensions of the analysis.

Setting

The performance of the binary splitting algorithm crucially depends on that of integer multiplication. Following common usage, we denote by M(n) a bound on the time needed to multiply two integers of at most n bits. Currently the best theoretical bound [START_REF] Fürer | Faster integer multiplication[END_REF]

is M(n) = O(n(lg n) exp O(lg * n)), where lg * n = min{k lg •k n 1}.
In practice, implementations such as GMP [START_REF] Granlund | GNU Multiple Precision Arithmetic Library[END_REF] use variants of the Schönhage-Strassen algorithm of complexity O(n(lg n)(lg lg n)). We make the usual assumption [START_REF] Zur Gathen | Modern Computer Algebra[END_REF] that the function n → M(n)/n is nondecreasing. It follows that M(n) + M(m) M(n + m). We also assume that the space complexity of integer multiplication is linear, which is true for the standard algorithms. Write à = É(i), and define the bit size of a number (x + iy)/w ∈ à (where w, x, y ∈ ) as ⌈lg w⌉ + ⌈lg x⌉ + ⌈lg y⌉ + 1. Consider a linear differential equation with coefficients in Ã(z). It will prove convenient to clear all denominators (both polynomial and integer) and multiply the equation by a power of z to write it as [START_REF] Bernstein | Fast multiplication and its applications[END_REF] a r (z) z d dz

r + • • • + a 1 (z)z d dz + a 0 (z) • y(z) = 0, a k ∈ [i][z].
Let s = max k deg a k , and let h 1 denote the maximum bit size of the coefficients of the a k . Although our complexity estimates depend on r and h 1 , we do not consider more general dependencies on the equation. Thus, the a k are assumed to vary only in ways that can be described in terms of these two parameters. Specifically, we assume that s = O(1) and that the coefficients of a k (z)/a r (0) are all restricted to some bounded domain.

We also assume that 0 is an ordinary (i.e. nonsingular) point of (1). This implies that a r (0) = 0 and s r. The case of regular singular points (those for which we still have a r (0) = 0 but possibly s < r [13, Chap. 9]) is actually similar [START_REF] Van Der Hoeven | Fast evaluation of holonomic functions near and in regular singularities[END_REF][START_REF] Mezzarobba | Autour de l'évaluation numérique des fonctions Dfinies[END_REF]; we focus on ordinary points to avoid cumbersome notations.

Let ρ = min{|z| : a r (z) = 0} ∈ (0, ∞]. Then any formal series solution y(z) = n 0 y n z n of (1) converges on the disk |z| < ρ. We select a particular solution (say, by specifying initial values y(0), . . . , y (r-1) (0) in some fixed, bounded domain), and an evaluation point ζ ∈ Ã with |ζ| < ρ. Let h 2 denote the bit size of ζ, and

let h = h 1 + h 2 .
Again, h 2 is allowed to grow to infinity, but we assume that |ζ| is bounded away from ρ.

Given p 0, our goal is to compute a complex number ω ∈ Ã such that

|ω -y(ζ)| 2 -p
. By a classical argument, which can be reconstructed by substituting a series with indeterminate coefficients into (1), the sequence (y n ) obeys a recurrence relation of the form

(2) b 0 (n)y n+r + b 1 (n)y n+r-1 + • • • + b s (n)y n+r-s = 0, b j ∈ Ã [n] . Writing a k (z) = a k,0 + a k,1 z + • • • + a k,s z s , the b j are given explicitly by (3) b j (n) = r k=0 a k,j (n + r -j) k .
Based on the matrix form of the recurrence (2), set ( 4)

B(n) = ζC(n) 0 R 1 ∈ Ã(n) (s+1)×(s+1)
where

C(n) =      1 . . . 1 -bs(n) b0(n) • • • • • • -b1(n) b0(n)      , R = 0 . . . 0 s-r zeroes 1 0 . . . 0 r-1 zeroes . Let P (a, b) = B(b -1) • • • B(a + 1)B(a)
for all a b. (In particular, P (a, a) is the identity matrix.)

Algorithm 1. BinSplit(a, b) 1 If b -a (some threshold) 2 Return B(b -1) • • • B(a)
where B is defined by ( 5)

3 else 4 Return BinSplit(⌊ a+b 2 ⌋, b) • BinSplit(a, ⌊ a+b 2 ⌋)
One may check that b 0 (n) = 0 for n 0, due to the fact that 0 is an ordinary point of [START_REF] Bernstein | Fast multiplication and its applications[END_REF]. Thus the computation of a partial sum S N = N -1 n=0 y n ζ n reduces to that of the matrix product P (0, N ). Indeed, we have

(y n+r-s ζ n , . . . , y n+r-1 ζ n , S n ) T = P (0, n) (y r-s , . . . , y r-1 , 0) T
where y r-s = 0, . . . , y -1 = 0, y 0 , . . . , y r-1 are easily determined from the initial values of the differential equation.

Review of the Classical Binary Splitting Algorithm

Since the entries of the matrix B(n) are rational functions of n, the bit size of A slight complication stems from the fact that removing common divisors between the numerators and denominators of the fractions appearing in the intermediate P (a, b) ∈ Ã r×r would in general be too expensive. Multiplying the numer- ators and denominators separately and doing a single final division yields better complexity bounds. Let

P (a, b) grows as O((b -a) lg b) when b, (b -a) → ∞. This bound is sharp in the sense that it is reached for some (in fact, most) differential equations. Computing P (a, b) as B(b -1) • [B(b -2) • [• • • B(a)]]
(5) B(n) = b 0 (n) ζB(n) ∈ [i][n] (s+1)×(s+1) , ζ = ζ/ ζ ( ζ ∈ [i] , ζ ∈ ).
The entries of B(n) are polynomials of degree at most r and bit size O(h). To compute P (a, b) by binary splitting, we multiply the B(n) for a n < b using Algorithm 1, and then divide the resulting matrix by its bottom right entry. The general algorithm considered here was first published by Chudnovsky and Chudnovsky [START_REF] Chudnovsky | Computer algebra in the service of mathematical physics and number theory[END_REF], with (up to minor details) the analysis summarized in Prop. 1. The idea of binary splitting was known long before [START_REF] Gosper | Strip mining in the abandoned orefields of nineteenth century mathematics[END_REF][START_REF] Bernstein | Fast multiplication and its applications[END_REF]. 1) , both bounds are sharp.

Proposition 1. [5] As b, N = b -a, h, r → ∞ with r = O(N ), Algorithm 1 com- putes an unreduced fraction equal to P (a, b) in O(M N (h+r lg b) lg N ) operations, using O N (h + r lg b) bits of memory. Assuming M(n) = n(lg n)(lg lg n) O(
sketch. The bit sizes of the matrices that get multiplied together at any given depth 0 δ < ⌈lg N ⌉ in the recursive calls are at most C2 -δ N (h+d lg b) for some C. Since there are at most 2 δ such products and the multiplication function M(•) was assumed to be subadditive, the contribution of each level is bounded by M(C(ba)(h + d lg b)), whence the total time complexity. See [START_REF] Chudnovsky | Computer algebra in the service of mathematical physics and number theory[END_REF][START_REF] Mezzarobba | Autour de l'évaluation numérique des fonctions Dfinies[END_REF] for details. The intermediate results stored or multiplied together at any stage of the computation are disjoint subproducts of B(b -1) • • • B(a), and we assumed the space complexity of n-bit integer multiplication to be O(n), so the space required by the algorithm is linear in the combined size of the B(n). Finally, it is not hard to construct examples of differential equations that reach these bounds.

Remark 1. The link between our setting and the more common description of the algorithm for hypergeometric series is as follows. In the notation of Haible and Pananikolaou [START_REF] Haible | Fast multiprecision evaluation of series of rational numbers[END_REF] also used in Yakhontov's article, the partial sums of the hypergeometric series are related to its defining parameters a, b, p, q by s(i + 1)

S(i) = p(i) q(i) 0 a(i) b(i) p(i) q(i) b(i)q(i) s(i) S(i -1) , s(i) = b(i) a(i) s(i).
This equation becomes

B(i) T (0,i) = b(i)p(i) 0 a(i)p(i) b(i)q(i) B(i-1)P (i-1) T (0,i-1)
upon clearing denominators. The standard recursive algorithm for hypergeometric series may be seen an "inlined" computation of the associated product tree. Each recursive step is equivalent to the computation of the matrix product Br Pr 0 Tr BrQr

B l P l 0 T l B l Q l .
We return to the evaluation of a D-finite power series within its disk of convergence. From the differential equation ( 1), suitable initial conditions, the evaluation point ζ and a target precision p, one can compute [START_REF] Mezzarobba | Effective bounds for P-recursive sequences[END_REF] a truncation order N such that |S Ny(ζ)| 2 -p and ( 6)

N ∼ Kp = lg(|ζ| /ρ) -1 p, if ρ < ∞ N = Θ(p/ lg p), if ρ = ∞.
Combined with these estimates, Proposition 1 implies the following. This result is the basis of more general evaluation algorithms for D-finite functions [START_REF] Chudnovsky | Computer algebra in the service of mathematical physics and number theory[END_REF]. Indeed, binary splitting can be used to compute the required series sums at each step when solving a differential equation of the form (1) by the so-called method of Taylor series [START_REF] Mathews | Bibliography for Taylor series method for D.E.'s[END_REF]. Corollary 1 thus extends to the evaluation of y outside the disk |z| < ρ. Chudnovsky and Chudnovsky further showed how to reduce the cost of evaluation from Ω(hp) = Ω(p 2 ) to softly linear in p when h = Θ(p). This last situation is very natural since it covers the case where the point ζ is itself a O(p)-digits approximation resulting from a previous computation. The method, known as the bit burst algorithm, consists in solving the differential equation along a path made of approximations of ζ of exponentially increasing precision. Its time complexity is O(M(p(lg p) 2 )) [START_REF] Mezzarobba | NumGfun: a package for numerical and analytic computation with Dfinite functions[END_REF]. The improvements from the next section apply to all these settings. See also [START_REF] Van Der Hoeven | Transséries et analyse complexe effective[END_REF] for an overview of more sophisticated applications.

"Truncated" Binary Splitting

The superiority of binary splitting over alternatives like summing the series in floating-point arithmetic results from the controlled growth of intermediate results. Indeed, in the product tree computed by Algorithm 1, the exact representations of most subproducts P (a, b) are much more compact than Θ(p)-digits approximations would be. However, as already mentioned, the bit sizes of the P (a, b) also grow larger than p near the root of the tree. The size of a subproduct appearing at depth δ is roughly 2 -δ N (h + r lg N ). Assuming N = Θ(p), this means that the intermediate results get significantly larger than the output in the top Θ(lg lg p) levels of the tree.

A natural remedy is to use a hybrid of binary splitting and naive summation. More precisely, we split the full product P (0, N ) into ∆ = Θ(ln N ) subproducts of O(p) bits each, which are computed by binary splitting. The results are accumulated by successive multiplications at precision O(p).

We make use of the following notations to state and analyze the algorithm. In Equations ( 7) to [START_REF] Granlund | GNU Multiple Precision Arithmetic Library[END_REF] below, the coefficients of a general matrix A ∈ k×k are denoted a p,q = x p,q + iy p,q (1 p, q k) with x p,q , y p,q ∈ Ê. Let • be a submultiplicative norm on k×k , and let β k > 0 be such that

(7) A β k N (A), N (A) = max{|x i,j | , |y i,j |} 1 i,j k .
For definiteness, assume for now that • = • 1 is the matrix norm induced by the vector 1-norm. (We will discuss this choice later.) Then it holds that ( 8)

N (A) A 1 = k max j=1 k i=1 |a i,j | √ 2kN (A)
and ( 9)

P (a, b) b-1 n=a B(n) b-1 n=a 1 + |ζ| + |ζ| s max k=1 b k (n) b 0 (n) .
Observe that, since 1 is an eigenvalue of B(n) and the norm • is assumed to be submultiplicative, we have B(n) 1 for all n. Besides, it is clear from (3) that B(n) is bounded.

Given a ∈ É and ε < 1, let [START_REF] Gourdon | Binary splitting method[END_REF] Trunc(a, ε) = sgn(a) ⌊2 e |a|⌋ 2 -e , e = lg ε -1 .

We have |Trunc(a, ε) -a| ε; the size of Trunc(a, ε) is O(lg ε -1 ) for bounded a; and Trunc(a, ε) may be computed in O(M(h + e)) bit operations where h is the bit size of a. We extend the definition to matrices A ∈ Ã k×k by [START_REF] Granlund | GNU Multiple Precision Arithmetic Library[END_REF] Trunc(A, ε) = Trunc(x p,q , β -1 k ε) + i Trunc(y p,q , β -1 k ε) 1 p,q k , so that again Trunc(A, ε) -A ε. Note that we often write expressions of the form Trunc(a⋆b, ε) for some operator ⋆. Though this does not affect our complexity bounds, it is usually better to compute the approximate value of a⋆b directly instead of starting with an exact computation and truncating the result. See Brent and Zimmermann [START_REF] Brent | Modern Computer Arithmetic[END_REF] for some relevant algorithms.

The complete binary splitting algorithm with truncations is stated as Algorithm 2. Its key properties are summarized in the following propositions. Proposition 2. The output P = TruncBinSplit(p) of Algorithm 2 is such that P -P (0, N ) 2 -p .

Proof. Set P (q) = P (0, ⌊ q ∆ N ⌋) and Q (q) = P (⌊ q ∆ N ⌋, ⌊ q+1 ∆ N ⌋). Then, for 0 q ∆, it holds that ( 12) P (q) -P (q) q ∆ ε M ∆-q .

Algorithm 2. TruncBinSplit(p)

The notation X (q) , q = 0, 1, . . . refers to a single memory location X at different points q of the computation. [START_REF] Van Der Hoeven | Fast evaluation of holonomic functions[END_REF][START_REF] Mezzarobba | Effective bounds for P-recursive sequences[END_REF] 3 Set ∆ = ⌈ N p (h + r lg N )⌉, where h and r are given following Eq. ( 1)

1 Set ε = 2 -p 2 Compute N such that |S N -y(ζ)| ε
4 Compute M such that max ∆-1 q=0 P (⌊ q ∆ N ⌋, ⌊ q+1 ∆ N ⌋) + ε M C N/∆
, where C does not depend on p, h, r [say, by approximating the right-hand side of ( 9) from above with O(lg p) bits of precision]

5 Initialize P (0) := id ∈ Ã (s+1)×(s+1) 6 For q = 0, 1, . . . , ∆ -1 7 Q = ( Qi,j ) := BinSplit( q ∆ N , q+1 ∆ N ) (Algorithm 1) 8 Q(q) := Trunc( Q-1 s+1,s+1 • Q, 1 2∆ M -∆+1 ε) 9 P (q+1) := Trunc( Q(q) • P (q) , 1 2∆ M -∆+q+1 ε) 10 Return P (∆)
Indeed, this is true for q = 0. After Step 8 of each loop iteration, we have the bound Q(q) -Q

(q) 1 2∆ M -∆+1 ε ε since B(n)
1 for all n. Using ( 12) and the inequality Q(q) M from Step 3, it follows that Q(q) P (q) -Q (q) P (q) Q(q) -Q (q) P (q) + Q(q) P (q) -P (q) 2q + 1 2∆ ε M ∆-q-1 . After taking into account the truncation error from Step 9, we obtain P (q+1) -P (q+1) = P (q+1) -Q (q) P (q) q + 1 ∆ ε M ∆-q-1 which concludes the induction. Proposition 3. Not counting the cost of Step 2, Algorithm 2 runs in time

(13) O M(p)(h + r lg p) lg p , if ρ < ∞, O M(p)(h + r lg p) , if ρ = ∞,
as p, h, r → ∞ with r = O(lg p) and h = O(p). In both cases, it uses O(p) bits of memory (where the hidden constant is independent of h and r, under the same growth assumptions).

We neglect the cost of finding N to avoid a lengthy discussion of the complexity of the corresponding bound computation algorithms. It could actually be checked to be polynomial in r and lg p.

Time

Space (classical) Space (trunc.) 1. Complexity of some D-finite function evaluation algorithms based on binary splitting. The rows labeled "BinSplit" summarize the cost of computing a single sum by binary splitting, with or without truncations. Those labeled "BitBurst" refer to the computation of y(ζ) by the "bit burst" method, using either of Algorithm 1 and Algorithm 2 at each step. All entries are asymptotic bounds as p, h → ∞ with h = O(p). In the "BinSplit" case, we also let r tend to infinity under the assumption that r = O(lg p). The whole point of the "bit burst" method is to get rid the dependency on h.

ρ < ∞ BinSplit O(M(p(h + r lg p) lg p)) O(p(h + r lg p)) O(p) BitBurst O(M(p(lg p) 2 )) O(p lg p) O(p) ρ = ∞ BinSplit O(M(p(h + r lg p))) O(p(r + h/ lg p)) O(p) BitBurst O(M(p(lg p) 2 )) O(p) O(p) Table
The truncations in Steps 8 and 9 ensure that the bit sizes of P and Q are always at most [START_REF] Kreckel | decimal(γ)[END_REF] lg The required memory comprises space for the current values of P (q) and Q (q) , any temporary storage used by the operations from Steps 7 to 9, and an additional O(lg p) bits to manipulate auxiliary variables such as M and q. We have seen that P (q) and Q (q) have bit size O(p). Besides, our assumption that fast integer multiplication could be performed in linear space implies the same property for division by Newton's method. Thus, Steps 8 and 9 use O(p) bits of auxiliary storage. Finally, again by Proposition 1, the calls to Algorithm 1 use O((N/∆)(h+ r lg p)) = O(p) bits of memory.

ε -1 + ∆ lg M + lg ∆ + O(1) = O(p).
Plugging Algorithm 2 into the numerical evaluation algorithms mentioned at the end of Sect. 3 yields corresponding improvements for the evaluation of D-finite functions at more general points. Table 1 summarizes the complexity bounds we obtain. The omitted proofs are direct adaptations of those that apply without truncations [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][START_REF] Mezzarobba | Autour de l'évaluation numérique des fonctions Dfinies[END_REF]. There would be much to say on the hidden constant factors. The main result may be stated more precisely as follows.

Theorem 1. Let U ⊂ be a simply connected domain such that 0 ∈ U and a r (z) = 0 for all z ∈ U . Fix ℓ 0 , . . . , ℓ r-1 ∈ and ζ ∈ U . Assume that 0 is an ordinary point of (1), and let y be the unique solution of (1) defined on U and such that y (k) (0) = ℓ k , 0 k < r. Then, the value y(ζ) may be computed with error bounded by 2 -p in time O(M(p)(lg p) 2 ) and space O(p), not counting the resources needed to approximate the ℓ k or ζ to precision O(p) or to find suitable truncation orders for the Taylor series involved.

Finally, some comments are in order regarding the "working precision", that is, the size p ′ of the entries of P and Q in Algorithm 2. Equation ( 14) suggests a number of "guard digits" p ′p = Θ(p). Moreover, if the bound M is computed using [START_REF] Gourdon | Constants and records of computation[END_REF], the hidden constant depends on the choice of • .

Let B ∞ = lim n→∞ B(n). For the norm • opt given by Lemma 1 below, we have Proof. We mimic the classical proof of Householder's theorem [START_REF] Serre | Matrices[END_REF]Sect. 4.2]. By (3), the limit C ∞ = lim n→∞ C(n) is the companion matrix of the polynomial z s a r (1/z). The eigenvalues of ζC ∞ are strictly smaller than 1 in absolute value since |ζ| < ρ.

Let Γ be such that Γ -1 C ∞ Γ is in (lower) Jordan normal form. Let λ > 0, and set Π = diag(Γ, 1) • diag(1, λ, . . . , λ s ). Then Π -1 B ∞ Π is lower triangular, with off-diagonal entries tending to zero as λ → 0. Hence we have Π -1 B ∞ Π 1 = 1 for λ small enough. We choose such a λ (e.g., λ = 1-|ζ|/ρ 2 max(1,|ζ|) ) and set A opt = Π -1 AΠ 1 .

One way to eliminate the overestimation in the algorithm is to compute approximations of the matrices P (⌊ q ∆ N ⌋, ⌊ q+1 ∆ N ⌋) with O(lg p) digits of precision before doing the computation at full precision. One then uses the norms of these approximate products instead of those of the individual B(n) to determine M . We can also explicitly construct an approximation Π of the matrix Π from the proof of Lemma 1 precise enough that Π-1 B ∞ Π 1 = 1, and use the corresponding norm instead of • 1 in (9). (Compare [22, Algorithm B].) Other options include computing symbolic bounds on the coefficients of P (a, b) as a function of a and b [START_REF] Mezzarobba | Effective bounds for P-recursive sequences[END_REF] or finding an explicit integer n 0 such that n n 0 ⇒ B(n) opt = 1 based on the symbolic expression of n. Which variant to use in practice depends on the features of the implementation platform.

In any case, replacing the O(•) in the space complexity bound by an explicit constant would also require more specific assumptions on the memory representation of the objects we work with, as well as finer control on the space complexity of integer multiplication and division (see, e.g., Roche [START_REF] Roche | Efficient Computation with Sparse and Dense Polynomials[END_REF]).

Final Remarks

What we lose and what we retain. The price we pay for the reduced memory usage is the ability to easily extend the computation to higher precision. Indeed, the classical algorithm computes the exact value of the matrix P (0, N ), from which we can deduce P (0, N ′ ) for any N ′ > N in time roughly proportional to N ′ -N . This is no longer true with the linear-space variant. In some "lucky" cases where P (0, N ) can be represented exactly in linear space, it is possible to get the memory usage down to O(N ) while preserving restartability: see Cheng et al. [START_REF] Cheng | Time-and space-efficient evaluation of some hypergeometric constants[END_REF] and the references therein. Additionally, the resulting running time is reportedly lower than using truncations, probably owing to the fact that the size of the subproducts in the lg(N/∆) lower levels of the tree is reduced as well. Unfortunately, the applicability of the technique is limited to very special cases.

Two other traditional selling points of the binary splitting method are its easy parallelization and good memory locality. Nothing is lost in this respect, except that the memory bound grows to Θ(t • p) when using t = o(lg N ) parallel tasks in the approximate part of the computation.

Generalizations. The idea of binary splitting "with truncations" and the outline of its analysis adapt to various settings not covered here. For instance, we may consider systems of linear differential equations instead of scalar equations [START_REF] Chudnovsky | Computer algebra in the service of mathematical physics and number theory[END_REF]. Product trees of matrices over number fields à ′ = É(α) other than É(i) or over rings of truncated power series à ′ [[ε]] / ε k are also useful, respectively, to evaluate limits of D-finite functions at regular singular points of their defining equations, and to make the analytic continuation process more efficient for equations of large order [START_REF] Van Der Hoeven | Fast evaluation of holonomic functions[END_REF][START_REF] Mezzarobba | Autour de l'évaluation numérique des fonctions Dfinies[END_REF]. It is not essential either that the coefficients of the recurrence relation satisfied by the y n are rational functions of n: all we really ask is that they have suitable growth properties and can be computed fast.

Implementation. We are working on an implementation of the algorithm from Sect. 4 in an experimental branch of the software package NumGfun [START_REF] Mezzarobba | NumGfun: a package for numerical and analytic computation with Dfinite functions[END_REF]. The current state of the code is available from http://marc.mezzarobba.net/supplementary-material/trunc-CASC2012/. A comparison (updated periodically) with the implementation of binary splitting without truncations used in previous releases of NumGfun is also included.

  then takes time at least quadratic in ba, as can be seen from the combined size of the intermediate results. The term "binary splitting" refers to the technique of reorganizing the product into a balanced tree of subproducts, using the relation P (a, b) = P (m, b) • P (a, m) with m = ⌊ 1 2 (a + b)⌋, and so on recursively.

Corollary 1 .

 1 Write ℓ = h + r lg p. Under the assumptions of Proposition 1, one can compute y(ζ) in O(M(ℓp lg p)) bit operations, using O(ℓp) bits of memory. The complexity goes down to O(M(ℓp)) operations and O(ℓp/ lg p) bits of memory when a r (z) is a constant.

Proof.

  Computing the bound M using Equation (9) as suggested is more than enough to ensure that lg M = O(N/∆). It requires O(N ) arithmetic operations on O(lg p)-bit numbers, that is, o(N (lg N ) 2 ) bit operations. By Proposition 1, each of the ∆ calls to BinSplit requires O M( N ∆ (h + r lg N )) lg N = O M(p) lg p bit operations. The resulting matrices Q (p) all have size O(p), hence the divisions from Step 8 can be done in O(M(p)) operations using Newton's method [25, Chap. 9].

  It follows that the matrix multiplications from Step 9 take O(M(p)) operations each. Summing up, each iteration of the loop from Step 6 can be performed in O(M(p) lg p) operations, for a total of O(∆M(p) lg p). Equation (13) follows upon setting N = O(p) or N = O(p/ lg p) according to (6).

lg P (a, b) opt b- 1 n=aLemma 1 .

 11 lg B ∞ opt + O(n -1 ) = O lg(ba) , and hence lg P (a, b) = O(lg(ba)) for any norm • . There exists a matrix norm • opt such that B ∞ opt = 1.

Acknowledgments. I would like to thank Nicolas Brisebarre and Bruno Salvy for encouraging me to write this note and offering useful comments, and Anne Vaugon for proofreading parts of it.


Charles Sanson

Odile Diou

Julie Thévenot

Emmanuel Ibarboure

Alain Soum

Annie Brûlet

Sylvain Miraux

Eric Thiaudière

Sisareuth Tan

Alain Brisson

Vincent Dupuis

Olivier Sandre
email: olivier.sandre@ipb.fr

Sébastien Lecommandoux
email: lecommandoux@enscbp.fr

Doxorubicin Loaded Magnetic Polymersomes: Theranostic Nanocarriers for MR Imaging and Magneto-Chemotherapy

Hydrophobically modified magnetic nanoparticles (MNPs) were encapsulated within the membrane of poly(trimethylene carbonate)-b-poly(L-glutamic acid) (PTMC-b-PGA) block copolymer vesicles using a nanoprecipitation process. This formulation method provides a high loading of MNPs (up to 70 wt %) together with a good control over the sizes of the vesicles (100 -400 nm). The deformation of the vesicle membrane under an applied magnetic field was evidenced by anisotropic SANS. These hybrid objects display contrast enhancement properties in Magnetic Resonance Imaging, a diagnostic method routinely used for three-dimensional and noninvasive scans of the human body. 1 They can also be guided in a magnetic field gradient. The feasibility of drug release triggered by magnetic induction was evidenced using the anticancer drug doxorubicin (DOX), which is co-encapsulated in the membrane. Magnetic polymersomes are thus proposed as multimodal drug nanocarriers for bioimaging and magneto-chemotherapy. In the new field of nanomedicine, a great variety of polymeric drug delivery nanocarriers have shown efficient entrapment and controlled release of drugs in vitro. However, evaluating the fate of the carrier after injection in vivo is a delicate task: as a consequence, an imaging probe is usually coencapsulated with the drugs in the polymer nanoparticles. Such multifunctional nanocarriers for cancer diagnostics and treatment open the field of "theranostics", i.e. combination of imaging and therapeutic functions in "all-in-one" nanoparticles. 2 Among the different nanoparticles used in biotechnological applications, superparamagnetic nanocrystals made of γ-Fe 2 O 3 around 10 nm in size are of particular interest since they are biocompatible and offer contrast enhancement in MRI. In addition, MNPs can also be implemented to induce cancer cells apoptosis. Local heating is obtained by exposure to radio-frequency magnetic fields in the frequency range 100 kHz -1MHz, due to thermal dissipation by the magnetic moments forced to oscillate at the frequency of the applied magnetic field. Hyperthermia therapy (i.e. treating by heat) has been recognized as a promising form of cancer therapy, particularly in synergy with chemo-or radio-therapy. 3 For a solid tumor located not too deep under the skin (such as melanoma), a conventional hyperthermia treatment is applied, the cancer cells being destroyed by the heat shock. We proposed a totally different strategy with a polymer vesicle delivering a drug through the temperature increase induced by radiofrequency magnetic fields not macroscopically but locally, at the scale of the membrane: the MNPs being confined inside the polymer block, the heat dissipates in the vicinity of the polymer membrane, which strongly impacts its permeability. Therefore we called this strategy "magneto-chemotherapy" as depicted on Fig. 1. 4 Fig. 1: Dually-loaded vesicles prepared by addition at controlled speed of an aqueous buffer into a mixture in DMSO of PTMC-b-PGA (m=19 and n=24), doxorucin and magnetic nanoparticles coated by the BEY surfactant of given formula (bottom sketch). The cryo-TEM image (top right panel) shows a vesicle wrapped by a mantle of MNPs, which excitation by a RF magnetic field heats up the membrane locally and accelerates the DOX release.

In this highlight paper, we describe a convenient way to prepare fully biodegradable polymersomes featuring a dual magnetic and thermo-responsive membrane. The hydrophobic core of membranes is made of the PTMC semi-crystalline blocks and the hydrophilic leaflets arise from the PGA polypeptide blocks of the copolymer. One pot nanoprecipitation (Fig. 1) leads to dual-loaded vesicles with iron oxide MNPs of appropriate size and coating (BEY) embedded within the membranes together with the antitumor drug doxorubicin. The amphiphilic diblock copolymer was synthesized by ringopening polymerization of γ-benzyl-L-glutamate Ncarboxyanhydride initiated by amino functionalized PTMC macroinitiator upon a previously published method. 5 PTMC 24 -b-PGA 19 presents a hydrophilic weight fraction of 50 wt % for an overall molar mass M n = 4900 g/mol with a dispersity index 1.15. On the other hand, magnetic iron oxides MNPs were synthesized in water by alkaline coprecipitation of ferrous and ferric salts followed by coating with a surfactant called Beycostat (bottom sketch on Fig 1). The gyration and hydrodynamic radii of these MNPs were respectively R G MNP =3.05±0.06 nm and R H MNP =4.70±0.07 nm as measured by SANS and DLS.

Nanoprecipitation allowed reaching quantitative loading contents (LC) with controlled final sizes and low polydispersity indexes, as reported on this table for two durations of the buffer addition (PBS pH=7.4), with and without DOX.

The two-dimensional confinement of the MNPs inside the vesicular membranes was evidenced by small angle neutron and light scattering techniques and observed by transmission electron microscopy. The magnetic membranes exhibited a reversible deformation under a permanent magnetic field of moderate intensity (0.1 Tesla) as seen on Fig. 2. Fig. 2 Top: SANS patterns of magnetic vesicles loaded at 50 wt. % MNPs in H 2 O in the q range 3×10 -3 -3×10 -2 Å -1 at zero field and under B=0.1 Tesla. Colors correspond to iso-intensity ranges. Patterns were averaged in angular sectors (-30° 30°) along the field direction or (-15° 15°) around the perpendicular direction. Bottom: Kratky-Porod fits Ln[q 2 I(q)] ~ -q 2 δ 2 /12 in the low q regime lead to membrane thicknesses in the two directions δ // and δ ┴ and to a scheme of the overall deformation.

With pure H 2 O as solvent, the nuclear scattering of the inorganic component of membranes was much greater than that of the polymer leaflets. Thus a variation in intensity in both directions can be ascribed to a reorganization of the MNPs within the membrane. The variation of δ // and δ ┴ as a function of the field intensity is interpreted by the stretching of the membrane near the magnetic poles due to an overall deformation of the vesicle into an elongated ellipsoid. The apparent increase of thickness near the equator and concomitant decrease near the magnetic poles (δ ┴ = 18 nm, δ // = 8 nm at B=0.6 T) indicates that MNPs move away from the poles where dipolar repulsions are strong and accumulate in the other portions of membrane where dipolar interactions are expected to be attractive Finally, the release of the anticancer drug was assayed under several heating conditions. A global heating at 37°C (melting temperature of PTMC 6 ) caused a significantly higher release than the one observed at 23°C, as explained by a more permeable membrane. But the application of an oscillating magnetic field at 23°C induced hyperthermia at a local (nanometric) scale only in the vicinity of the membrane that provoked an increase in DOX release (by 2 times), thereby evidencing the validity of the concept of magnetochemotherapy with those magnetic polymersomes.

Funding

Financial support was provided by the European Commission (FP7 CP-IP 213631-2 NANOTHER).

Fig. 3 :

 3 Fig. 3: In vitro release assays of doxorubicin at two temperatures (23 or 37°C) and with or without an applied radio-frequency magnetic field B=2.65 mT at 500 kHz.