Alin Bostan

Bruno Salvy

Computer algebra for Combinatorics

> value(subs(_alpha^(-20*k)=1,%)): > simplify(%) assuming k::posint: > factor(%); 2 (5 k + 1) * on a PC, (Intel Xeon X5160, 3GHz processor, with 8GB RAM), running Magma V2.16-7

Discrete Fourier Transform

[Gauss 1866, Cooley-Tukey 1965] DFT Problem: Given n = 2 k , f ∈ K[x] <n , and ω ∈ K a primitive n-th root of unity, compute (f (1), f (ω), . . . , f (ω n-1 ))

Then f (ω j ) = f even (ω 2j ) + ω j f odd (ω 2j ), and (ω 2j ) 0≤j<n = n 2 -roots of 1.

are they, other than just names for mathematical objects that are useful only in situations of contrived simplicity? Why are we so pleased when a complicated calculation "comes out" as a Bessel function, or a Laguerre polynomial? What determines which functions are "special"? -Sir Michael Berry [Ber01] People who like this sort of thing will find this the sort of thing they like.

-Barry Hughes, quoting Abraham Lincoln at the beginning of an appendix on "Special Functions for Random Walk Problems" [Hug95, p. 569] Problem 6

A flea starts at (0, 0) on the infinite two-dimensional integer lattice and executes a biased random walk: At each step it hops north or south with probability 1/4, east with probability 1/4 + , and west with probability 1/4 -. The probability that the flea returns to (0, 0) sometime during its wanderings is 1/2. What is ?

Asking for the that gives a certain probability p of return yields a problem hardly any more difficult than calculating the probability for a given : it just adds the need to use a numerical root-finder. But the problem looks more interesting the way it is stated. In §6.1 we give a short argument, why the problem is solvable.

We will discuss several methods for calculating the probability of return. In §6.2, using virtually no probability theory, we transform the problem to one of lin-Computer algebra conjectures and proves

p( ) = 1 - A 2 • 2 F 1   1 2 , 1 2 1 2 √ 1 -16 2 A   -1
, with A = 1 + 8 2 + 1 -16 2 .

Algebraic balanced urns

Computer algebra conjectures and proves larger classes of algebraic balanced urns.

More in Basile Morcrette's talk!

Gessel's conjecture

• Gessel walks: walks in N 2 using only steps in S = { , , ←, →}

• g(i, j, n) = number of walks from (0, 0) to (i, j) with n steps in S Question: Nature of the generating function

G(x, y, t) = ∞ i,j,n=0 g(i, j, n) x i y j t n ∈ Q[[x, y, t]]
Computer algebra conjectures and proves:

Theorem [B. & Kauers 2010] G(x, y, t) is an algebraic function † and G(1, 1, t) = 1 2t • 2 F 1 -1/12 1/4 2/3 - 64t(4t + 1) 2 (4t -1) 4 - 1 2t .
This is equivalent to finding the number M 20k of solutions (a, b, c) ∈ N 3 of a + 2b + 4c = 20k.

Euler-Comtet's denumerants:

n≥0 M n x n = 1 (1 -x)(1 -x 2 )(1 -x 4 ) . > f:=1/(1-x)/(1-x^2)/(1-x^4): > S:=series(f,x,201): > [seq(coeff(S,x,20*k),k=1..10)];
[ 36,121,256,441,676,961,1296,1681,2116,2601] > subs(n=20*k,gfun[ratpolytocoeff](f,x,n)):

17 32 + (20k + 1)(20k + 2) 16 + 5k + (-1) -20k (20k + 1) 16 + 5(-1) -20k 32 + α 2 +1=0 - ( 1 16 -1 16 α)α -20k α INTRODUCTION 2.

Computer Algebra

General framework

Computeralgebra = effectivemathematicsandalgebraiccomplexity

• Effective mathematics: what can we compute?

• their complexity: how fast?

Computer algebra books

Mathematical Objects

• Main objects

-integers Z -polynomials K[x] -rational functions K(x) -power series K[[x]] -matrices M r (K)
linear recurrences with constant, or polynomial, coefficients K[n] S n linear differential equations with polynomial coefficients K[x] ∂ x where K is a field (generally supposed of characteristic 0 or large)

• Secondary/auxiliary objects

-polynomial matrices M r (K[x]) -power series matrices M r (K[[x]])

Overview

Complexity yardsticks

Important features:

• addition is easy: naive algorithm already optimal

• multiplication is the most basic (non-trivial) problem

• almost all problems can be reduced to multiplication

Are there quasi-optimal algorithms for:

• integer/polynomial/power series multiplication? Yes!

• What can be computed in 1 minute with a CA system * polynomial product † in degree 14,000,000 (>1 year with schoolbook)

product of two integers with 500,000,000 binary digits factorial of N = 20, 000, 000 (output of 140,000,000 digits)

gcd of two polynomials of degree 600,000 resultant of two polynomials of degree 40,000 factorization of a univariate polynomial of degree 4,000 factorization of a bivariate polynomial of total degree 500 resultant of two bivariate polynomials of total degree 100 (output 10,000)

product/sum of two algebraic numbers of degree 450 (output 200,000) determinant (char. polynomial) of a matrix with 4,500 (2,000) rows determinant of an integer matrix with 32-bit entries and 700 rows

Inverse DFT

IDFT Problem:

Given n = 2 k , v 0 , . . . , v n-1 ∈ K and ω ∈ K a primitive n-th root of unity, compute f ∈ K[x] <n such that f (1) = v 0 , . . . , f (ω n-1 ) = v n-1 • V ω • V ω -1 = n • I n → performing the inverse DFT in size n amounts to: -performing a DFT at 1 1 , 1 ω , • • • , 1 ω n-1
-dividing the results by n.

• this new DFT is the same as before:

1 ω i = ω n-i ,
so the outputs are just shuffled.

Consequence: the cost of the inverse DFT is O(n log(n))

FFT polynomial multiplication

Suppose the basefield K contains enough roots of unity To multiply two polynomials f, g in K[x], of degrees < n:

• find N = 2 k such that h = f g has degree less than N N ≤ 4n • compute DFT(f, N ) and DFT(g, N ) O(N log(N ))
• multiply pointwise these values to get DFT(h, N ) O(N )

• recover h by inverse DFT O(N log(N ))

Complexity: O(N log(N )) = O(n log(n))
General case: Create artificial roots of unity O(n log(n) log log n)

Example: fast factorial

Problem: Compute N ! = 1 × • • • × N Naive iterative way: unbalanced multiplicands Õ(N 2 )
• Binary Splitting: balance computation sequence so as to take advantage of fast multiplication (operands of same sizes):

N ! = (1 × • • • × N/2 ) size 1 2 N log N × (( N/2 + 1) × • • • × N ) size 1 2 N log N
and recurse. Complexity Õ(N ).

• Extends to matrix factorials

A(N )A(N -1) • • • A(1) Õ(N )
-→ recurrences of arbitrary order.

Application to recurrences

Problem: Compute the N -th term u N of a P -recursive sequence

p r (n)u n+r + • • • + p 0 (n)u n = 0, (n ∈ N)
Naive algorithm: unroll the recurrence Õ(N 2 ) bit ops.

Binary splitting: U n = (u n , . . . , u n+r-1 ) T satisfies the 1st order recurrence

U n+1 = 1 p r (n) A(n)U n with A(n) =         p r (n) . . . p r (n) -p 0 (n) -p 1 (n) . . . -p r-1 (n)         . =⇒ u N reads off the matrix factorial A(N -1) • • • A(0)
[Chudnovsky-Chudnovsky, 1987]: Binary splitting strategy Õ(N ) bit ops.

Application: fast computation of e = exp(1)

[Brent 1976] e n = n k=0 1 k! -→ exp(1) = 2.7182818284590452 . . . Recurrence e n -e n-1 = 1/n! ⇐⇒ n(e n -e n-1 ) = e n-1 -e n-2 rewrites   e N -1 e N   = 1 N   0 N -1 N + 1   C(N )   e N -2 e N -1   = 1 N ! C(N )C(N -1) • • • C(1)   0 1   .
x κ+1 = N (x κ ) = x κ -(x 2 κ -2)/(2x κ ), x 0 = 1
x 1 = 1.5000000000000000000000000000000

x 2 = 1.4166666666666666666666666666667

x 3 = 1.4142156862745098039215686274510

x 4 = 1.4142135623746899106262955788901

x 5 = 1.4142135623730950488016896235025

Newton's tangent method: power series case Newton's tangent method: power series case

x κ+1 = N (x κ ) = x κ -(x 2 κ -(1 -t))/(2x κ ), x 0 = 1 x 1 = 1 - 1 2 t x 2 = 1 - 1 2 t - 1 8 t 2 - 1 16 t 3 - 1 32 t 4 - 1 64 t 5 - 1 128 t 6 - 1 256 t 7 - 1 512 t 8 - 1 1024 t 9 + • • • x 3 = 1 - 1 2 t - 1 8 t 2 -
In order to solve ϕ(x, g) = 0 in K [[x]] (where ϕ ∈ K[[x, y]], ϕ(0, 0) = 0 and ϕ y (0, 0) = 0), iterate

g κ+1 = g κ - ϕ(g κ ) ϕ y (g κ ) mod x 2 κ+1 g -g κ+1 = g -g κ + ϕ(g) + (g κ -g)ϕ y (g) + O((g -g κ ) 2 ) ϕ y (g) + O(g -g κ ) = O((g -g κ ) 2 ).
The number of correct coefficients doubles after each iteration 1966, Sieveking 1972& Kung 1974, Brent 1975] Division, logarithm and exponential of power series in K [[x]] can be computed at precision N using O(M(N )) operations in K Division, logarithm and exponential of power series

Total cost = 2 × the cost of the last iteration Theorem [Cook
[Sieveking1972, Kung1974, Brent1975]
To compute the reciprocal of f ∈ K[[x]] with f (0) = 0, choose ϕ(g) = 1/gf :

g 0 = 1/f 0 and g κ+1 = g κ + g κ (1 -f g κ ) mod x 2 κ+1 for κ ≥ 0. Complexity: C(N ) = C(N/2) + O(M(N )) =⇒ C(N ) = O(M(N ))
Corollary: division of power series at precision N in O(M(N ))

Division, logarithm and exponential of power series

[Sieveking1972, Kung1974, Brent1975]
To compute the reciprocal of f ∈ K[[x]], choose ϕ(g) = 1/gf :

g 0 = 1/f 0 and g κ+1 = g κ + g κ (1 -f g κ ) mod x 2 κ+1 for κ ≥ 0. Complexity: C(N ) = C(N/2) + O(M(N )) =⇒ C(N ) = O(M(N ))
Corollary: division of power series at precision N in O(M(N ))

Corollary: Logarithm log(f ) = - i≥1 (1 -f ) i i of f ∈ 1 + xK[[x]] in O(M(N )): • compute the Taylor expansion of h = f /f modulo x N -1 O(M(N ))
• take the antiderivative of h O(N )

Division, logarithm and exponential of power series

[Sieveking1972, Kung1974, Brent1975]
To compute the reciprocal of f ∈ K[[x]], choose ϕ(g) = 1/gf :

g 0 = 1/f 0 and g κ+1 = g κ + g κ (1 -f g κ ) mod x 2 κ+1 for κ ≥ 0. Complexity: C(N ) = C(N/2) + O(M(N )) =⇒ C(N ) = O(M(N ))
Corollary: division of power series at precision N in O(M(N ))

Corollary: Logarithm log(f ) = - i≥1 (1 -f ) i i of f ∈ 1 + xK[[x]] in O(M(N )): • compute the Taylor expansion of h = f /f modulo x N -1 O(M(N ))
• take the antiderivative of h O(N )

Corollary: Exponential exp(f ) = i≥0 f i i! of f ∈ xK[[x]]. Use φ(g) = log(g) -f : g 0 = 1 and g κ+1 = g κ -g κ (log(g κ ) -f ) mod x 2 κ+1 for κ ≥ 0. Complexity: C(N ) = C(N/2) + O(M(N )) =⇒ C(N ) = O(M(N ))
Application: Euclidean division for polynomials [Strassen, 1973] Pb:

Given F, G ∈ K[x] ≤N , compute (Q, R) in Euclidean division F = QG + R Naive algorithm: O(N 2 ) Idea: look at F = QG + R from infinity: Q ∼ +∞ F/G Let N = deg(F ) and n = deg(G). Then deg(Q) = N -n, deg(R) < n and F (1/x)x N rev(F ) = G(1/x)x n rev(G) • Q(1/x)x N -n rev(Q) + R(1/x)x deg(R) rev(R) •x N -deg(R)
Algorithm:

• Compute rev(Q) = rev(F )/rev(G) mod x N -n+1 O(M(N )) • Recover Q O(N ) • Deduce R = F -QG O(M(N ))
Application: conversion coefficients ↔ power sums [Schönhage, 1982] Any polynomial

F = x n + a 1 x n-1 + • • • + a n in K[x]
can be represented by its first n power sums

S i = F (α)=0 α i
Conversions coefficients ↔ power sums can be performed

• either in O(n 2 ) using Newton identities (naive way):

ia i + S 1 a i-1 + • • • + S i = 0, 1 ≤ i ≤ n • or in O(M(n)) using generating series rev(F ) rev(F ) = - i≥0 S i+1 x i ⇐⇒ rev(F ) = exp   - i≥1 S i i x i   Application: special bivariate resultants [B-Flajolet-S-Schost, 2006]
Composed products and sums: manipulation of algebraic numbers

F ⊗ G = F (α)=0,G(β)=0 (x -αβ), F ⊕ G = F (α)=0,G(β)=0 (x -(α + β))
Output size:

N = deg(F ) deg(G)
α s β s = (αβ) s and (α + β) s s! x s = α s s! x s β s s! x s Corollary: Fast polynomial shift P (x + a) = P (x) ⊕ (x + a) O(M(deg(P )))
Newton iteration on power series: operators and systems

In order to solve an equation φ(Y ) = 0, with φ

: (K[[x]]) r → (K[[x]]) r , 1. Linearize: φ(Y κ -U ) = φ(Y κ ) -Dφ| Y κ • U + O(U 2 ),
where Dφ| Y is the differential of φ at Y .

Iterate:

Y κ+1 = Y κ -U κ+1
, where U κ+1 is found by 3. Solve linear equation:

Dφ| Y k • U = φ(Y κ ) with val U > 0.
Then, the sequence Y κ converges quadratically to the solution Y .

Application: inversion of power series matrices [Schulz, 1933] To compute the inverse Z of a matrix of power series Y ∈ M r (K[[x]]):

• Choose the map φ :

Z → I -Y Z with differential Dφ| Y : U → -Y U • Equation for U : -Y U = I -Y Z κ mod x 2 κ+1 • Solution: U = -Y -1 (I -Y Z κ ) = -Z κ (I -Y Z κ ) mod x 2 κ+1
This yields the following Newton-type iteration for Y -1

Z κ+1 = Z κ + Z κ (I r -Y Z κ ) mod x 2 κ+1 Complexity: C inv (N ) = C inv (N/2) + O(MM(r, N )) =⇒ C inv (N ) = O(MM(r, N ))
Application: non-linear systems

In order to solve a system

Y = H(Y ) = φ(Y ) + Y , with H : (K[[x]]) r → (K[[x]]) r , such that I r -∂H/∂Y is invertible at 0. 1. Linearize: φ(Y κ -U ) -φ(Y κ ) = U -∂H/∂Y (Y κ ) • U + O(U 2 ). 2. Iterate Y κ+1 = Y κ -U κ+1
, where U κ+1 is found by 3. Solve linear equation:

(I r -∂H/∂Y (Y κ )) • U = H(Y κ ) -Y κ with val U > 0.
This yields the following Newton-type iteration:

   Z κ+1 = Z κ + Z κ (I r -(I r -∂H/∂Y (Y κ ))Z κ ) mod x 2 κ+1 Y κ+1 = Y κ -Z κ+1 (H(Y κ ) -Y κ ) mod x 2 κ+1
computing simultaneously a matrix and a vector.

Example: Mappings

> mappings:={M=Set(Cycle(Tree)),Tree=Prod(Z,Set(Tree))}: > combstruct[gfeqns](mappings,labeled,x);

[M (x) = 1 1 -Tree(x) , Tree(x) = x exp(Tree(x))]
> countmappings:=SeriesNewtonIteration(mappings,labelled,x): > countmappings(10); To compute the solution

M = 1 + x + 2
Y ∈ M r (K[[x]]) of the system Y = AY • choose the map φ : Y → Y -AY , with differential φ. • the equation for U is U -AU = Y κ -AY κ mod x 2 κ+1
• the method of variation of constants yields the solution

U = Y κ V κ mod x 2 κ+1 , Y κ -AY κ = Y κ V κ mod x 2 κ+1
This yields the following Newton-type iteration for Y :

Y κ+1 = Y κ -Y κ Y -1 κ (Y κ -AY κ ) mod x 2 κ+1 Complexity: C solve (N ) = C solve (N/2) + O(MM(r, N )) =⇒ C solve (N ) = O(MM(r, N ))
TOOLS FOR CONJECTURES 1. Hermite-Padé Approximants

Definition

Definition: Given a column vector

F = (f 1 , . . . , f n ) T ∈ K[[x]] n and an n-tuple d = (d 1 , . . . , d n ) ∈ N n , a
Hermite-Padé approximant of type d for F is a row vector P = (P 1 , . . . , P n ) ∈ K[x] n , (P = 0), such that:

(1)

P • F = P 1 f 1 + • • • + P n f n = O(x σ ) with σ = i (d i + 1) -1, (2) deg(
P i ) ≤ d i for all i.
σ is called the order of the approximant P.

Very useful concept in number theory (transcendence theory):

• [Hermite, 1873]: e is transcendent.

• [Lindemann, 1882]: π is transcendent, and so does e α for any α ∈ Q \ {0}.

• [Beukers, 1981]: reformulate Apéry's proof that Rivoal, 2000]: there exist an infinite number of k such that ζ(2k + 1) / ∈ Q.

ζ(3) = n 1 n 3 is irrational. • [

Worked example

Let us compute a Hermite-Padé approximant of type (1, 1, 1) for (1, C, C 2 ), where

C(x) = 1 + x + 2x 2 + 5x 3 + 14x 4 + 42x 5 + O(x 6 ).
This boils down to finding α 0 , α 1 , β 0 , β 1 , γ 0 , γ 1 such that

α 0 +α 1 x+(β 0 +β 1 x)(1+x+2x 2 +5x 3 +14x 4 )+(γ 0 +γ 1 x)(1+2x+5x 2 +14x 3 +42x 4 ) = O(x 5 ).
By identifying coefficients, this is equivalent to a homogeneous linear system:

         
1 0 1 0 1 0 0 1 1 1 2 1 0 0 2 1 5 2 0 0 5 2 14 5 0 0 14 5 42 14

          ×              α 0 α 1 β 0 β 1 γ 0 γ 1              = 0 ⇐⇒          
1 0 1 0 1 0 1 1 1 2 0 0 2 1 5 0 0 5 2 14 0 0 14 5 42

          ×           α 0 α 1 β 0 β 1 γ 0           = -γ 1           0 1 2 5 14          
.

By homogeneity, one can choose γ 1 = 1. Then, the violet minor shows that one can take (β 0 , β 1 , γ 0 ) = (-1, 0, 0). The other values are α 0 = 1, α 1 = 0.

Thus the approximant is (1, -1, x), which corresponds to P = 1y + xy 2 such that P (x, C(x)) = 0 mod x 5 .

Algebraic and differential approximation = guessing

• Hermite-Padé approximants of n = 2 power series are related to Padé approximants, i.e. to approximation of series by rational functions

• algebraic approximants = Hermite-Padé approximants for f = A -1 , where

A ∈ K[[x]]
seriestoalgeq, listtoalgeq

• differential approximants = Hermite-Padé approximants for f = A ( -1) , where

A ∈ K[[x]]
seriestodiffeq, listtodiffeq > listtoalgeq ([1,1,2,5,14,42,132,429],y(x)); 2 [1 -y(x) + x y(x) , ogf] > listtodiffeq ([1,1,2,5,14,42,132,429] Ideas:

• Compute a whole matrix of approximants

• Exploit divide-and-conquer Algorithm: > proba:=proc(i,j,n,c) option remember; if abs(i)+abs(j)>n then 0 elif n=0 then 1 else expand(proba(i-1,j,n-1,c)*(1/4+c)+proba(i+1,j,n-1,c)*(1/4-c) +proba(i,j+1,n-1,c)*1/4+proba(i,j-1,n-1,c)*1/4) fi end: > seq(proba(0,0,k,c),k=0..6);

1. If σ = n(d + 1) -1 ≤ threshold, call the naive algorithm 2. Else: (a) recursively compute P 1 ∈ K[x] n×n s.t. P 1 • F = O(x σ/2 ), deg(P 1 ) ≈ d 2 (b) compute "residue" R such that P 1 • F = x σ/2 • R + O(x σ/2 ) (c) recursively compute P 2 ∈ K[x] n×n s.t. P 2 • R = O(x σ/2 ), deg(P 2 ) ≈ d 2 (
1, 0, 1 4 -2c 2 , 0, 9 64 - 9 4 c 2 + 6c 4 , 0, 25 256 - 75 32 c 2 + 15c 4 -20c 6
> gfun:-listtodiffeq([seq(proba(0,0,2*k,c),k=0..20)],y(x));

[ -1 + 8 c 2 + 48 xc 4 y (x) + 4 -8 x + 64 xc 2 + 192 x 2 c 4 d dx y (x) + 4 x + 64 x 3 c 4 -4 x 2 + 32 x 2 c 2 d 2 dx 2 y (x) , y (0) = 1, D (y) (0) = 1/4 -2 c 2 , ogf ]

Linear differential operators

Definition: If K is a field, K x, ∂; ∂x = x∂ + 1 , or simply K(x) ∂ , denotes the associative algebra of linear differential operators with coefficients in K(x).

  multiplication in K[x] <n , and of n-bit integers = O(n 2 ) by the naive algorithm = O n 1.58 by Karatsuba's algorithm = O n log α (2α-1) by the Toom-Cook algorithm (α ≥ 3) = O n log n loglog n by the Schönhage-Strassen algorithm MM(r) = complexity of matrix product in M r (K) = O(r 3 ) by the naive algorithm = O(r 2.81 ) by Strassen's algorithm = O(r 2.38 ) by the Coppersmith-Winograd algorithm MM(r, n) = complexity of polynomial matrix product in M r (K[x] <n ) = O(r 3 M(n)) by the naive algorithm = O(MM(r) n log(n) + r 2 n log n loglog n) by the Cantor-Kaltofen algo = O(MM(r) n + r 2 M(n)) by the B-Schost algorithm Fast polynomial multiplication in practice Practical complexity of Magma's multiplication in F p [x], for p = 29 × 2 57 + 1.

e

  N in Õ(N ) bit operations[Brent 1976] generalizes to the evaluation of any D-finite series at an algebraic number[Chudnovsky-Chudnovsky 1987] Õ(N ) bit ops.Implementation in gfun[Mezzarobba, S. 2010] > rec:={n*(e(n) -e(n-1)) = e(n-1) -e(n-2), e(0)=1, e(1)=2}; > pro:=rectoproc(rec,e(n)); pro := proc(n::nonnegint) local i1, loc0, loc1, loc2, tmp2, tmp1, i2; if n <= 22 then loc0 := 1; loc1 := 2; if n = 0 then return loc0 else for i1 to n -1 do loc2 := (-loc0 + loc1 + loc1*(i1 + 1))/(i1 + 1); loc0 := loc1; loc1 := loc2 end do end if; loc1 else tmp1 := 'gfun/rectoproc/binsplit'([ 'ndmatrix'(Matrix([[0, i2 + 2], [-1, i2 + 3]]), i2 + 2), i2, 0, n, matrix_ring(ad, pr, ze, ndmatrix(Matrix(2, 2, [[...],[...]], datatype = anything, storage = empty, shape = [identity]), 1)), expected_entry_size], Vector(2, [...], datatype = anything)); tmp1 := subs({e(0) = 1, e(1) = 2}, tmp1); tmp1 end if end proc > tt:=time(): x:=pro(50000): time()-tt, evalf(x-exp(1), 200000); 1.827, 0. Application: record computation of π [Chudnovsky-Chudnovsky 1987] fast convergence hypergeometric identity n (6n)!(13591409 + 545140134n) n! 3 (3n)!(8 • 100100025 • 327843840) n . Used in Maple & Mathematica: 1st order recurrence, yields 14 correct digits per iteration -→ 4 billion digits [Chudnovsky-Chudnovsky 1994] Current record on a PC: 10000 billion digits [Kondo & Yee 2011]HIGH PRECISION3. Newton IterationNewton's tangent method: real case[Newton, 1671]

  Linear algebra: χ xy , χ x+y in K[x, y]/(F (x), G(y)) O(MM(N )) Resultants: Res y F (y), y deg(G) G(x/y) , Res y (F (y), G(xy)) O(N 1.5 ) Better: ⊗ and ⊕ are easy in Newton representation O(M(N ))

Theorem

  For any vector F = (f 1 , . . . , f n ) T ∈ K[[x]] n and for any n-tuple d = (d 1 , . . . , d n ) ∈ N n , there exists a Hermite-Padé approximant of type d for F.Proof: The undetermined coefficients of P i = d i j=0 p i,j x j satisfy a linear homogeneous system with σ = i (d i + 1) -1 equations and σ + 1 unknowns.Corollary Computation in O(MM(σ)) = O(σ θ ), for 2 ≤ θ ≤ 3.There are better algorithms:• The linear system is structured (Sylvester-like / quasi-Toeplitz)• Derksen's algorithm (Gaussian-like elimination) O(σ 2 ) • Beckermann-Labahn's algorithm (DAC) Õ(σ) = O(σ log 2 σ)Quasi-optimal computationTheorem[Beckermann-Labahn, 1994] One can compute a Hermite-Padé approximant of type(d, . . . , d) for F = (f 1 , . . . , f n ) in O(MM(n, d) log(nd)).

  d) return P := P 2 • P 1 The precise choices of degrees is a delicate issue Gcd, extended gcd, Padé approximants in O(M(n) log n) Example: Flea from SIAM 100-Digit Challenge

  Next steps: dsolve (+ help) and evaluation at x = 1.TOOLS FOR CONJECTURES2. p-Curvature of Differential OperatorsImportant classes of power series algebraic hypergeom D-finite seriesAlgebraic: S(x) ∈ K[[x]] root of a polynomial P ∈ K[x, y]. D-finite: S(x) ∈ K[[x]] satisfying a linear differential equation with polynomial (or rational function) coefficients c r (x)S (r) (x) + • • • + c 0 (x)S(x) = 0. Hypergeometric: S(x) = n s n x n such that s n+1 s n ∈ K(n). E.g. n (b) n (c) n x n n!, (a) n = a(a + 1) • • • (a + n -1).

K

  [x] ∂ is called the (rational) Weyl algebra. It is the algebraic formalization of the notion of linear differential equation with rational function coefficients:a r (x)y (r) (x) + • • • + a 1 (x)y (x) + a 0 (x)y(x) = 0 ⇐⇒ L(y) = 0, where L = a r (x)∂ r + • • • + a 1 (x)∂ + a 0 (x) The commutation rule ∂x = x∂ + 1 formalizes Leibniz's rule (f g) = f g + f g .Implementation in the DEtools package: diffop2de, de2diffop, mult DEtools[mult](Dx,x,[Dx,x]); x Dx + 1 Weyl algebra is Euclidean Theorem [Libri 1833, Brassinne 1864, Wedderburn 1932, Ore 1932] K(x) ∂ is a non-commutative (left and right) Euclidean domain: for any A, B ∈ K(x) ∂ , there exist unique operators Q, R ∈ K(x) ∂ such that A = QB + R, and deg ∂ (R) < deg ∂ (B). This is called the Euclidean right division of A by B. Moreover, any A, B ∈ K(x) ∂ admit a greatest common right divisor (GCRD) and a least common left multiple (LCLM). They can be computed by a non-commutative version of the extended Euclidean algorithm. rightdivision, GCRD, LCLM from the DEtools package > rightdivision(Dx^10,Dx^2-x,[Dx,x])+ 80) Dx + 100 x + xproves that Ai (10) (x) = (20x 3 + 80)Ai (x) + (100x 2 + x 5 )Ai(x) of a series are enough to guess candidate differential equations below the red curve. GCRD of candidates could jump above the red curve.

The Grothendieck-Katz p-curvatures conjecture Q: when does a differential equation possess a basis of algebraic solutions? E.g. for the Gauss hypergeometric equation x(1 -x)∂ 2 + (γ -(α + β + 1)x)∂ -αβx, Schwarz's list (1873) classifies algebraic 2 F 1 's in terms of α, β, γ

Conjecture [Grothendieck, 1960's, unpublished;Katz, 1972] Let A ∈ Q(x) n×n . The system (S) : y = Ay has a full set of algebraic solutions if and only if, for almost all prime numbers p, the system (S p ) defined by reduction of (S) modulo p has a full set of algebraic solutions over F p (x).

Definition:

The p-curvature of (S) is the matrix A p , where

Theorem [Cartier, 1957] The sufficient condition of the G.-K. Conjecture is equivalent to A p = 0 mod p.

For each p, this can be checked algorithmically.

Grothendieck's conjecture

Q: when does a differential equation possess a basis of algebraic solutions?

For a scalar differential equation, the G.-K. Conjecture can be reformulated:

The equation (E) : L(y) = 0 has a basis of algebraic solutions if and only if, for almost all prime numbers p, the operator L right-divides ∂ p modulo p.

For each p, this can be checked algorithmically.

Conjecture is proved for

Grothendieck's conjecture for combinatorics

Suppose that we have guessed a linear differential equation L(f ) = 0 (by differential Hermite-Padé approximation) for some power series

and that we want to recognize whether f is algebraic or not.

Recipe 1: try algebraic guessing.

Recipe 2: For several primes p, compute p-curvatures mod p, and check whether they are zero; equivalently, test if ∂ p mod L = 0 (mod p).

For many power series coming from counting problems (diagonals, constant terms, integrals of algebraic functions, . . . ) Grothendieck's conjecture is true.

Grothendieck's conjecture at work

(excerpt from Rodriguez-Villegas's "Integral ratios of factorials") Algebraicity of u can be however guessed using any prior knowledge, by computing p-curvatures of the (minimal) order-8 operator L s.t. L(u) = 0 For p < 300, they are all zero, except when p ∈ {11, 13, 17, 19, 23}

G-series and global nilpotence

(4) OGF of any P -recursive, integer-valued, exponentially bounded, sequence

Theorem [Chudnovsky 1985] The minimal-order linear differential operator annihilating a G-series is globally nilpotent: for almost all prime numbers p, it right-divides ∂ pµ modulo p, for some µ ≤ deg ∂ L.

(this condition is equivalent to the nilpotence mod p of the p-curvature matrix)

Global nilpotence for combinatorics

Suppose we have guessed (by differential approximation) a linear differential equation L(f ) = 0 for a power series f ∈ Q [[x]] which is a G-series (typically, the OGF of a P -recursive, integer-valued, exponentially bounded, sequence).

A way to empirically certify that L is very plausible: