Towards fast Puiseux series computation Adrien Poteaux , Marc - PowerPoint PPT Presentation

Towards fast Puiseux series computation Adrien Poteaux ⋆ , Marc Rybowicz † ⋆ : LIFL - Université Lille 1 † : XLIM - Université de Limoges Computer algebra and polynomials Workshop, Linz November 25th, 2013 adrien.poteaux@lifl.fr Puiseux series 1 / 16

Algebraic plane curves: a projective point of view y 3 K = Q ( α ) a number field y ′ 2 y ′′ F ( X , Y ) ∈ K [ X , Y ] 2 C = { ( x , y ) ∈ C 2 | F ( x , y ) = 0 } Let x 0 ∈ C : y 2 y ′′ Fiber at x 0 : 1 F ( x 0 ) = { roots of F ( x 0 , Y ) = 0 } . y ′ 1 y 1 Regular point : # F ( x 0 ) = d Y . Critical point : # F ( x 0 ) < d Y . = ⇒ roots of R F = Res Y ( F , F Y ) . x 0 α 1 α 2 adrien.poteaux@lifl.fr Puiseux series 1 / 16

Puiseux series: a generalization of formal power series x 0 regular ; F ( x 0 ) = { y 1 , · · · , y d Y } . Theorem (Implicit function theorem) ∞ α ik ( X − x 0 ) k s.t � There are d Y series Y i ( X ) = k = 0 F ( X , Y i ( X )) = 0 around x 0 and Y i ( x 0 ) = y i . adrien.poteaux@lifl.fr Definition 2 / 16

Puiseux series: a generalization of formal power series x 0 regular ; F ( x 0 ) = { y 1 , · · · , y d Y } . Theorem (Implicit function theorem) ∞ α ik ( X − x 0 ) k s.t � There are d Y series Y i ( X ) = k = 0 F ( X , Y i ( X )) = 0 around x 0 and Y i ( x 0 ) = y i . x 0 critical Theorem (Puiseux) ∞ k � α ik ζ jk ei s.t. There are d Y series Y ij ( X ) = e i ( X − x 0 ) k = n i F ( X , Y ij ( X )) = 0 for all 1 ≤ j ≤ e i , 1 ≤ i ≤ s, with ζ e i primitive e i -th root of unity, e i e 1 , . . . , e s partition of d Y (ramification indices). adrien.poteaux@lifl.fr Definition 2 / 16

Long term goal Puiseux series � � � Monodromy group � � � Effective Abel-Jacobi theorem ւ ց Computer Algebra Physics Integral basis computation.s KP & KdV equations Algebraic solutions of ODEs Differential Galois theory adrien.poteaux@lifl.fr Motivations 3 / 16

Difficult part is the singular part ∞ k � α ik ζ jk S ij ( X − x 0 ) = e i ( X − x 0 ) ei k = n i r ij k � α ik ζ jk ei + next terms e i ( X − x 0 ) = k = n i r ij is the regularity index ; r i = r ij for 1 ≤ j ≤ e i Next terms can be computed using quadratic Newton iterations Kung & Traub 1978, All Algebraic Functions Can Be Computed Fast adrien.poteaux@lifl.fr Singular part 4 / 16

Singulart part computation: the Newton-Puiseux algorithm Newton, 1676 → introduction of the concept. Puiseux, 1850 → rediscovers ; first procedure. Chystov, 1986 → “Newton-Puiseux bit complexity is polynomial” . Duval, 1989 → rational algorithm ; arithmetic complexity O ( D 8 ) . Walsh, 2000 → bit complexity O ˜( D 36 ) (classical algorithm). Walsh, 1999 → polynomial size for rational coefficients, no algorithm . Rybowicz & P., 2008 → improved arithmetic complexity: O ˜( D 5 ) . adrien.poteaux@lifl.fr Newton-Puiseux algorithm 5 / 16

The Newton-Puiseux algorithm: main tools F ( X , Y ) = Y 7 + Y 5 X − 2 Y 4 X + 5 Y 4 X 3 + 4 Y 2 X 2 + X 6 adrien.poteaux@lifl.fr Newton-Puiseux algorithm 6 / 16

Support of a polynomial F ( X , Y ) = Y 7 X 0 + Y 5 X 1 − 2 Y 4 X 1 + 5 Y 3 X 4 − Y 3 X 2 + 4 Y 2 X 2 + Y 0 X 6 × Supp(F) = { ( i , j ) ∈ N 2 | a ij � = 0 } 6 4 2 1 0 2 3 4 5 7 adrien.poteaux@lifl.fr Newton-Puiseux algorithm 6 / 16

Newton polygon � a ij Y i X j F ( X , Y ) = i , j × Supp(F) = { ( i , j ) ∈ N 2 | a ij � = 0 } 6 N ( F ) — N ( F ) : lower part of the convex hull of Supp(F). 2 1 0 2 4 7 adrien.poteaux@lifl.fr Newton-Puiseux algorithm 6 / 16

Characteristic polynomial � a ij Y i X j F ( X , Y ) = i , j × Supp(F) = { ( i , j ) ∈ N 2 | a ij � = 0 } 6 N ( F ) — N ( F ) : lower part of the convex hull of Supp(F). Characteristic polynomial : 2 ∆ , slope − m q i − i 0 � 1 φ ∆ ( T ) = a ij T q ( i , j ) ∈ ∆ 0 i 0 = 2 4 7 adrien.poteaux@lifl.fr Newton-Puiseux algorithm 6 / 16

Rational Newton-Puiseux algorithm D. Duval 1989, Rational Puiseux Expansions For each edge ∆ of N ( F ) s φ M k � – φ ∆ = k k = 1 l q – For each φ k F ( X , Y ) ← F ( ξ v k X q , X m ( ξ u k + Y )) ∆ , slope − m q X l with 0 i 0 I ( F ) · ξ k s. t. φ k ( ξ k ) = 0, · ( u , v ) such that uq − vm = 1. adrien.poteaux@lifl.fr Newton-Puiseux algorithm 7 / 16

Rational Newton-Puiseux algorithm : first turn For each edge ∆ of N 0 ( F ) s φ M k � – φ ∆ = k k = 1 – For each φ k F ( X , Y ) ← F ( ξ v k X q , X m ( ξ u k + Y )) X l with ∆ , slope − m · ξ k s. t. φ k ( ξ k ) = 0, q i 0 d Y · ( u , v ) such that uq − vm = 1. l q First turn: initial polygon N 0 ( F ) adrien.poteaux@lifl.fr Newton-Puiseux algorithm 7 / 16

Pure symbolic computation is costly H ( X , Y ) = ( Y 3 − X ) (( Y − 1 ) 2 − X ) ( Y − 2 − X 2 ) + X 2 Y 5 R H ( X ) = X 3 P ( X ) , deg X ( P ) = 23; β s.t. P ( β ) = 0 Singular parts of Puiseux series of H above β : S i ( X ) = α i , 0 , 1 ≤ i ≤ 4. 1 2 , 5 ≤ i ≤ 6. S i ( X ) = α i , 0 + α i , 1 ( X − β ) Degree of the extension field i = 5 , 6: K ( α i , 0 ) = K ( α i , 1 ) = K ( β ) → extension of degree 23, i = 1 , . . . , 4: [ K ( α i , 0 ) : K ( β )] = 4 → extension of degree 92, Coefficient growth α i , 0 → rational number with 98 digits, α i , 1 → rational number with 132 digits. adrien.poteaux@lifl.fr Newton-Puiseux algorithm 8 / 16

Numerical computations ? Direct computation: almost useless Guessing the structure ? two difficulties: Finding the correct Newton polygon Factorising “well” φ ∆ 4 x 2 − 2 . 0 x + 0 . 9999 ? = ( x − 0 . 99 )( x − 1 . 01 ) 2 ? ( x − 1 . ) 2 = 1 = ⇒ Multiplicity structure ? 0 2 4 7 = ⇒ Exact informations needed ! adrien.poteaux@lifl.fr Newton-Puiseux algorithm 9 / 16

A symbolic-numeric approach: 1 Compute the singular part of Puiseux series modulo a well chosen prime number p This give us the structure of the Puiseux series, thus: Newton polygons, Multiplicity structures of the φ ∆ . 2 Use this information to conduct a numerical computation of the Puiseux series coefficients. adrien.poteaux@lifl.fr A modular-numeric approach 10 / 16

Modular part: main results Reduction criteria: One can reduce F mod p , 1 p > d Y 2 tc ( R F ) �≡ 0 mod p . 3 If p satisfies that: Puiseux series can be reduced modulo p , 1 The structure computed modulo p is the good one. 2 Bounds for p ; log ( p ) ≃ log ( D ) with probabilistic algorithms. References for details: Poteaux & Rybowicz 2008, On the good reduction of Puiseux series and complexity of the Newton-Puiseux algorithm over finite fields ; Poteaux & Rybowicz 2012, On the good reduction of Puiseux series and Applications ; Poteaux & Rybowicz 2011, Complexity bounds for the rational Newton-Puiseux algorithm over finite fields and related problems adrien.poteaux@lifl.fr A modular-numeric approach 11 / 16

Numerical part: following the structure If P ( X ) = ( X − α ) m ( X − β ) m = ( X − ˜ γ ) m ( X − ˜ δ ) m , γ or α ↔ ˜ α ↔ ˜ δ ? no answer = ⇒ group of roots dealt together We separate the blocs later on: filter according to the Newton polygon 1 (coefficient of F zero or not) filter according to the multiplicity structure 2 first idea = approximate gcd : singular values One example adrien.poteaux@lifl.fr A modular-numeric approach 12 / 16

Fast computation ? 1 Use only “mandatory” monomials Truncating power of X 1 already used in the D 8 of Duval, improved in our D 5 version, better ? → relaxed computations ( a-posteriori bounds) Reducing the degree in Y 2 After first turn, we only look roots above ( 0 , 0 ) Get rid of roots above ( 0 , α � = 0 ) ? Factorization 2 Fast computation in between two branch separations S ( X ) = 2 + X − X 2 + X 3 + 3 X 2 + 4 X 4 − X 2 − 2 X 5 + X 7 9 11 35 2 + 7 X 6 + . . . adrien.poteaux@lifl.fr A fast algorithm 13 / 16

A new algorithm 1 Substitutions F ( X , Y ) ← F ( X , Y + A d Y − 1 ( X )) ( compute common terms at once ; less recursive calls) 2 Factorization of the polynomial during the algorithm ( Hensel lemma → recursive calls with smaller degrees) 3 Using relax algorithms. ( no a priori bound required) 4 Improved truncation bounds (thanks to relax algorithms). leads to O ˜( D 4 ) adrien.poteaux@lifl.fr A fast algorithm 14 / 16

Better again ? The idea: all series do not need the same truncations 1 Compute half of the series in O ˜( D 3 ) (the one that requires the smallest truncations) 2 Factorize F = G H with G corresponding to the computed series, and H to the other ones. O ˜( D 3 ) ? 3 Apply the same procedure to H . Logarithmic number of recursive call: total complexity in O ˜( D 3 ) . Remaining problem: step 2 ; how to get the starting point in order to apply Hensel lemma ? or another idea ? adrien.poteaux@lifl.fr A fast algorithm 15 / 16

Conclusion A modular-numerical approach Modular part: 100% done, Numerical part: first strategy developped, Maple implementation (prototype) To be done: certified computations, error bounds, better implementation. . . A new faster algorithm D 4 strategy looks to work, Missing one idea for the D 3 algorithm, C++ implementation (based on NTL) started, To be done: finding this last missing part, writting and coding everything properly. Long term: a good effective way to compute Abel’s map ? adrien.poteaux@lifl.fr A fast algorithm 16 / 16

Annexes adrien.poteaux@lifl.fr A fast algorithm 16 / 16