correct rounding of transcendental functions an approach
play

Correct rounding of transcendental functions: an approach via - PowerPoint PPT Presentation

Correct rounding of transcendental functions: an approach via Euclidean lattices and approximation theory Nicolas Brisebarre (C.N.R.S.) and Guillaume Hanrot (.N.S. Lyon) CUNY Graduate Center-Courant Seminar in Symbolic-Numeric Computing -


  1. The Table Maker’s Dilemma: an Example Consider the function 2 x and the binary64 FP number (base 2 , p = 53 ) x = 8520761231538509 2 62 We have 2 52+ x = 4509371038706515 . 4999999999999999994026 · · · (decimal) = 1 · · · . 01 · · · · · · · · · · · · 1 0 · · · (binary) . ���� � �� � 53 bits 60 consecutive 1 ′ s Hardest-to-round (HR) case for function 2 x and binary64 FP numbers. Lefèvre et al.: the value of m is 113 ( ∼ 2 p, p = 53 ) here. -10-

  2. The Table Maker’s Dilemma: an Example Consider the function 2 x and the binary64 FP number (base 2 , p = 53 ) x = 8520761231538509 2 62 We have 2 52+ x = 4509371038706515 . 4999999999999999994026 · · · (decimal) = 1 · · · . 01 · · · · · · · · · · · · 1 0 · · · (binary) . ���� � �� � 53 bits 60 consecutive 1 ′ s Hardest-to-round (HR) case for function 2 x and binary64 FP numbers. Lefèvre et al.: the value of m is 113 ( ∼ 2 p, p = 53 ) here. √ Function f : , sin , cos , arcsin , arccos , tan , arctan , exp , log , sinh , n cosh ... -10-

  3. The Table Maker’s Dilemma: Diophantine Problems Assume, w.l.o.g, that x and f ( x ) ∈ [1 , 2) . Q. (TMD) We want to determine m ∈ N , as small as possible, s.t. for all j ∈ � 2 p − 1 , 2 p − 1 � , � � j = 2 ℓ + 1 either there exists ℓ ∈ � 2 p − 1 , 2 p − 1 � s.t. f , 2 p − 1 2 p or � � � � � � j − 2 ℓ + 1 for all 2 p − 1 � ℓ � 2 p − 1 , � � � � 2 − m . � f 2 p − 1 2 p -11-

  4. The Table Maker’s Dilemma: First Challenge A breakpoint is a point where the rounding function changes. In this talk, it is the middle of two consecutive FP numbers. First challenge: Determine the set BP f of all the FP numbers x ∈ [1 , 2) such that f ( x ) is a breakpoint. In other words, determine all j, ℓ ∈ � 2 p − 1 , 2 p − 1 � s.t. � � j = 2 ℓ + 1 f . 2 p − 1 2 p -12-

  5. State of the Art A function f is algebraic if there exists P ∈ Z [ X, Y ] \ { 0 } s.t. P ( x, f ( x )) = 0 . Otherwise, f is transcendental. -13-

  6. State of the Art Algebraic functions. Works by Jeannerod, Louvet, Muller and Panhaleux. -14-

  7. State of the Art Algebraic functions. Works by Jeannerod, Louvet, Muller and Panhaleux. Transcendental elementary Functions sin , cos , arcsin , arccos , tan , arctan , exp , log , sinh , cosh . Hermite-Lindemann’s theorem: α � = 0 algebraic ⇒ e α transcendental. -14-

  8. State of the Art Algebraic functions. Works by Jeannerod, Louvet, Muller and Panhaleux. Transcendental elementary Functions sin , cos , arcsin , arccos , tan , arctan , exp , log , sinh , cosh . Hermite-Lindemann’s theorem: α � = 0 algebraic ⇒ e α transcendental. Therefore, let x a FP number, f ( x ) is not a breakpoint, except for straightforward cases ( e 0 , ln(1) , sin(0) , . . . ). -14-

  9. State of the Art Algebraic functions. Works by Jeannerod, Louvet, Muller and Panhaleux. Transcendental elementary Functions sin , cos , arcsin , arccos , tan , arctan , exp , log , sinh , cosh . Hermite-Lindemann’s theorem: α � = 0 algebraic ⇒ e α transcendental. Therefore, let x a FP number, f ( x ) is not a breakpoint, except for straightforward cases ( e 0 , ln(1) , sin(0) , . . . ). What about the Euler Gamma function? For Re( z ) > 0 , � + ∞ t z − 1 e − t d t. Γ( z ) = 0 Very little is known: Γ(1 / 2) , Γ(1 / 3) , Γ(1 / 4) , Γ(1 / 6) , Γ(2 / 3) , Γ(3 / 4) , Γ(5 / 6) are transcendental and there are some partial irrationality results. -14-

  10. Our setting Let f : [1 , 2) �→ [1 , 2) , f is transcendental and analytic in the neighborhood of [1 , 2) . -15-

  11. Our setting Let f : [1 , 2) �→ [1 , 2) , f is transcendental and analytic in the neighborhood of [1 , 2) . Let g ∈ C ([ a, b ]) , � g � ∞ , [ a,b ] := max x ∈ [ a,b ] | g ( x ) | . -15-

  12. Our Approach: Polynomial Interpolation and Lattice Basis Reduction We want to find all 2 p − 1 � i, j � 2 p − 1 s.t. � � i = 2 j + 1 f , 2 p − 1 2 p -16-

  13. Our Approach: Polynomial Interpolation and Lattice Basis Reduction We want to find all 2 p − 1 � i, j � 2 p − 1 s.t. � � i = 2 j + 1 f , 2 p − 1 2 p � � i i.e. 2 j + 1 = 2 p f . 2 p − 1 -16-

  14. Our Approach: Polynomial Interpolation and Lattice Basis Reduction We want to find all 2 p − 1 � i, j � 2 p − 1 s.t. � � i = 2 j + 1 f , 2 p − 1 2 p � � i i.e. 2 j + 1 = 2 p f . 2 p − 1 We build a trap! We find a partition of [1 , 2) = ∪ ℓ I ℓ such that, for all ℓ , we can compute P ℓ, 1 , P ℓ, 2 ∈ Z [ X, Y ] \ { 0 } with � � � P ℓ,k (2 p − 1 u, 2 p f ( u )) � < 1 , k = 1 , 2 , for all u ∈ I ℓ . -16-

  15. Our Approach: Polynomial Interpolation and Lattice Basis Reduction We build a trap! We find a partition of [1 , 2) = ∪ ℓ I ℓ such that, for all ℓ , we can compute P ℓ, 1 , P ℓ, 2 ∈ Z [ X, Y ] \ { 0 } with � � � P ℓ,k (2 p − 1 u, 2 p f ( u )) � < 1 , k = 1 , 2 for all u ∈ I ℓ . -17-

  16. Our Approach: Polynomial Interpolation and Lattice Basis Reduction We build a trap! We find a partition of [1 , 2) = ∪ ℓ I ℓ such that, for all ℓ , we can compute P ℓ, 1 , P ℓ, 2 ∈ Z [ X, Y ] \ { 0 } with � � � P ℓ,k (2 p − 1 u, 2 p f ( u )) � < 1 , k = 1 , 2 for all u ∈ I ℓ . For all ℓ , if there exist 2 p − 1 � i, j � 2 p − 1 s.t. i/ 2 p − 1 ∈ I ℓ and � � i = 2 j + 1 f , 2 p − 1 2 p -17-

  17. Our Approach: Polynomial Interpolation and Lattice Basis Reduction We build a trap! We find a partition of [1 , 2) = ∪ ℓ I ℓ such that, for all ℓ , we can compute P ℓ, 1 , P ℓ, 2 ∈ Z [ X, Y ] \ { 0 } with � � � P ℓ,k (2 p − 1 u, 2 p f ( u )) � < 1 , k = 1 , 2 for all u ∈ I ℓ . For all ℓ , if there exist 2 p − 1 � i, j � 2 p − 1 s.t. i/ 2 p − 1 ∈ I ℓ and � � i = 2 j + 1 f , 2 p − 1 2 p � �� � � �� � v u then we have, for k = 1 , 2 , P ℓ,k ( i, 2 j + 1) ∈ Z and | P ℓ,k ( i, 2 j + 1) | < 1! -17-

  18. Our Approach: Polynomial Interpolation and Lattice Basis Reduction We build a trap! We find a partition of [1 , 2) = ∪ ℓ I ℓ such that, for all ℓ , we can compute P ℓ, 1 , P ℓ, 2 ∈ Z [ X, Y ] \ { 0 } with � � � P ℓ,k (2 p − 1 u, 2 p f ( u )) � < 1 , k = 1 , 2 for all u ∈ I ℓ . For all ℓ , if there exist 2 p − 1 � i, j � 2 p − 1 s.t. i/ 2 p − 1 ∈ I ℓ and � � i = 2 j + 1 f , 2 p − 1 2 p � �� � � �� � v u then we have, for k = 1 , 2 , P ℓ,k ( i, 2 j + 1) ∈ Z and | P ℓ,k ( i, 2 j + 1) | < 1! ⇒ P ℓ,k ( i, 2 j + 1) = 0 . -17-

  19. Our Approach: Polynomial Interpolation and Lattice Basis Reduction We build a trap! We find a partition of [1 , 2) = ∪ ℓ I ℓ such that, for all ℓ , we can compute P ℓ, 1 , P ℓ, 2 ∈ Z [ X, Y ] \ { 0 } with � � � P ℓ,k (2 p − 1 u, 2 p f ( u )) � < 1 , k = 1 , 2 for all u ∈ I ℓ . For all ℓ , if there exist 2 p − 1 � i, j � 2 p − 1 s.t. i/ 2 p − 1 ∈ I ℓ and � � i = 2 j + 1 f , 2 p − 1 2 p � �� � � �� � v u then we have, for k = 1 , 2 , P ℓ,k ( i, 2 j + 1) ∈ Z and | P ℓ,k ( i, 2 j + 1) | < 1! ⇒ P ℓ,k ( i, 2 j + 1) = 0 . We eliminate (heuristic assumption!) one of the two variables and we get i and j , if they exist (Coppersmith; Boneh & Durfee; Stehlé). -17-

  20. Interpolation at Chebyshev Nodes and Uniform Approximation Definition Let n ∈ N , the Chebyshev nodes of the first kind of order n are the � � (2 k +1) π points µ k = cos , k = 0 , . . . , n . 2( n +1) -18-

  21. Interpolation at Chebyshev Nodes and Uniform Approximation Definition Let n ∈ N , the Chebyshev nodes of the first kind of order n are the � � (2 k +1) π points µ k = cos , k = 0 , . . . , n . 2( n +1) Let p n ∈ R n [ X ] that interpolates ϕ ∈ C ([ − 1 , 1]) at the ( µ k ) k =0 ,..,n , then we have � 1 � � ϕ − p n � ∞ , [ − 1 , 1] � 2 π log( n + 1) + 1 Q ∈ R n [ x ] � ϕ − Q � ∞ , [ − 1 , 1] . min -18-

  22. Interpolation at Chebyshev Nodes and Uniform Approximation Definition Let n ∈ N , the Chebyshev nodes of the first kind of order n are the � � (2 k +1) π points µ k = cos , k = 0 , . . . , n . 2( n +1) Let p n ∈ R n [ X ] that interpolates ϕ ∈ C n +1 ([ a, b ]) at the ( µ k ) k =0 ,..,n , then for any x ∈ [ − 1 , 1] , there exists ξ x ∈ ( − 1 , 1) s.t. ϕ ( x ) − p n ( x ) = ϕ ( n +1) ( ξ x ) 2 n ( n + 1)! . -18-

  23. Bernstein Ellipse � ρe iθ + ρ − 1 e − iθ � Let ρ > 1 , let E ρ := , θ ∈ [0 , 2 π ] . 2 -19-

  24. Bernstein Ellipse � ρe iθ + ρ − 1 e − iθ � Let ρ > 1 , let E ρ := , θ ∈ [0 , 2 π ] . 2 Bernstein ellipses for ρ = 1 . 05 , 1 . 25 , 1 . 45 , 1 . 65 , 1 . 85 . -19-

  25. Interpolation at Chebyshev Nodes and Uniform Approximation Definition Let n ∈ N , the Chebyshev nodes of the first kind of order n are the � � (2 k +1) π points µ k = cos , k = 0 , . . . , n . 2( n +1) Let ρ > 1 , let ϕ be a function analytic in a neighborhood of E ρ . Let p n ∈ R n [ X ] that interpolates ϕ at the ( µ k ) k =0 ,..,n , we have � ϕ − p n � ∞ , [ − 1 , 1] � 4 M ρ ( ϕ ) ρ n ( ρ − 1) . where M ρ ( ϕ ) = max z ∈E ρ | ϕ ( z ) | . -20-

  26. Interpolation at Chebyshev Nodes and Uniform Approximation Definition Let n ∈ N , the Chebyshev nodes of the first kind of order n are the � � (2 k +1) π points µ k = cos , k = 0 , . . . , n . 2( n +1) Let P ∈ R n [ X ] , we have � 2 � x ∈ [ − 1 , 1] | P ( x ) | � max π log( n + 1) + 1 k =0 ,...,n | P ( µ k ) | . max -20-

  27. Interpolation at Chebyshev Nodes and Uniform Approximation Definition Let n ∈ N , the Chebyshev nodes of the first kind of order n are the � � (2 k +1) π points µ k = cos , k = 0 , . . . , n . 2( n +1) Let ϕ ∈ C n +1 , let p n ∈ R n [ X ] that interpolates ϕ at the ( µ k ) k =0 ,..,n , let r n = � ϕ − p n � ∞ , [ − 1 , 1] , -20-

  28. Interpolation at Chebyshev Nodes and Uniform Approximation Definition Let n ∈ N , the Chebyshev nodes of the first kind of order n are the � � (2 k +1) π points µ k = cos , k = 0 , . . . , n . 2( n +1) Let ϕ ∈ C n +1 ([ − 1 , 1]) , let p n ∈ R n [ X ] that interpolates ϕ at the ( µ k ) k =0 ,..,n , let r n = � ϕ − p n � ∞ , [ − 1 , 1] , we have � 2 � � ϕ � ∞ � � p n � ∞ + r n � π log( n + 1) + 1 k =0 ,...,n | p n ( µ k ) | + r n max � 2 � π log( n + 1) + 1 k =0 ,...,n | ϕ ( µ k ) | + r n . max � -20-

  29. Interpolation at Chebyshev Nodes and Uniform Approximation Definition Let n ∈ N , the Chebyshev nodes of the first kind of order n are the � � (2 k +1) π points µ k = cos , k = 0 , . . . , n . 2( n +1) Let I = [ a, b ] , one can define scaled Chebyshev nodes of the first kind of � � (2 k +1) π order n : µ k, [ a,b ] = b − a + a + b 2 , k = 0 , . . . , n . 2 cos 2( n +1) -20-

  30. Lattice Basis Reduction -21-

  31. An Approach based on Lattice Basis Reduction Definition Let L be a nonempty subset of R d , L is a lattice iff there exists a set of vectors b 1 , . . . , b k R -linearly independent such that L = Z .b 1 ⊕ · · · ⊕ Z .b k . ( b 1 , . . . , b k ) is a basis of the lattice L . Examples. Z d , every subgroup of Z d . -22-

  32. Example: The Lattice Z (2 , 0) ⊕ Z (1 , 2) (1 , 2) (0 , 0) (2 , 0) -23-

  33. Example: The Lattice Z (2 , 0) ⊕ Z (1 , 2) (1 , 2) v − 3 u + v u (0 , 0) (2 , 0) 2 u − v -24-

  34. Example: The Lattice Z (2 , 0) ⊕ Z (1 , 2) (1 , 2) v − 3 u + v u (0 , 0) (2 , 0) 2 u − v SVP (Shortest Vector Problem) -24-

  35. Example: The Lattice Z (2 , 0) ⊕ Z (1 , 2) (1 , 2) v − 3 u + v u (0 , 0) (2 , 0) 2 u − v SVP (Shortest Vector Problem) is NP-hard. -24-

  36. Lenstra-Lenstra-Lovász Algorithm SVP (Shortest Vector Problem) is NP-hard. Factoring Polynomials with Rational Coefficients , A. K. Lenstra, H. W. Lenstra and L. Lovász, Math. Annalen 261 , 515-534, 1982. The LLL algorithm gives an approximate solution to SVP in polynomial time. -25-

  37. Lenstra-Lenstra-Lovász Algorithm Theorem Let L a lattice of rank k . LLL provides a basis ( b 1 , . . . , b k ) made of “pretty” short vectors. We have || b 1 || � 2 ( k − 1) / 2 λ 1 ( L ) where λ 1 ( L ) denotes the norm of a shortest nonzero vector of L . It terminates in at most O ( k 6 ln 3 B ) operations with B � || b i || 2 for all i . -26-

  38. Lenstra-Lenstra-Lovász Algorithm Theorem Let L a lattice of rank k . LLL provides a basis ( b 1 , . . . , b k ) made of “pretty” short vectors. We have || b 1 || � 2 ( k − 1) / 2 λ 1 ( L ) where λ 1 ( L ) denotes the norm of a shortest nonzero vector of L . It terminates in at most O ( k 6 ln 3 B ) operations with B � || b i || 2 for all i . Let ( b 1 , . . . , b k ) be an LLL-reduced basis L then || b 1 || � 2 k/ 2 ( vol L ) 1 /k || b 2 || � 2 k/ 2 ( vol L ) k − 1 . 1 and -26-

  39. Example: The Lattice Z (2 , 0) ⊕ Z (1 , 2) (1 , 2) − 3 u + v v u (0 , 0) (2 , 0) 2 u − v -27-

  40. Example: The Lattice Z (2 , 0) ⊕ Z (1 , 2) -28-

  41. Lenstra-Lenstra-Lovász Algorithm Theorem Let L a lattice of rank k . LLL provides a basis ( b 1 , . . . , b k ) made of “pretty” short vectors. We have || b 1 || � 2 ( k − 1) / 2 λ 1 ( L ) where λ 1 ( L ) denotes the norm of a shortest nonzero vector of L . It terminates in at most O ( k 6 ln 3 B ) operations with B � || b i || 2 for all i . Let ( b 1 , . . . , b k ) be an LLL-reduced basis L then 1 || b 1 || � 2 k/ 2 ( vol L ) 1 /k || b 2 || � 2 k/ 2 ( vol L ) k − 1 . and -29-

  42. How do we compute P 1 and P 2 ? Let P 1 = � u,v α u,v X u Y v and P 2 = � u,v β u,v X u Y v ∈ Z [ X, Y ] . We want to have � � � P k (2 p − 1 x, 2 p f ( x )) � < 1 , (1) k = 1 , 2 , for all x ∈ I = [ a, b ] . -30-

  43. How do we compute P 1 and P 2 ? Let P 1 = � u,v α u,v X u Y v and P 2 = � u,v β u,v X u Y v ∈ Z [ X, Y ] . We want to have � � � P k (2 p − 1 x, 2 p f ( x )) � < 1 , (1) k = 1 , 2 , for all x ∈ I = [ a, b ] . To do so, discretize (1): plug (scaled) Chebyshev nodes x k . Rough (but actually very reasonable) intuition: If P 1 and P 2 are small at these points, then they are small on the whole domain. -30-

  44. How do we compute P 1 and P 2 ? We want P 1 (2 p − 1 x k , 2 p f ( x k )) and P 2 (2 p − 1 x k , 2 p f ( x k )) small, i.e. the vectors �� � � � (2 p − 1 x k ) u (2 p f ( x k )) v � α u,v (2 p − 1 x k ) u (2 p f ( x k )) v = α u,v k � �� � u,v u,v k e u,v and �� � � � (2 p − 1 x k ) u (2 p f ( x k )) v � β u,v (2 p − 1 x k ) u (2 p f ( x k )) v = β u,v k � �� � u,v u,v k e u,v should be small. -31-

  45. How do we compute P 1 and P 2 ? We want P 1 (2 p − 1 x k , 2 p f ( x k )) and P 2 (2 p − 1 x k , 2 p f ( x k )) small, i.e. the vectors �� � � � (2 p − 1 x k ) u (2 p f ( x k )) v � α u,v (2 p − 1 x k ) u (2 p f ( x k )) v = α u,v k � �� � u,v u,v k e u,v and �� � � � (2 p − 1 x k ) u (2 p f ( x k )) v � β u,v (2 p − 1 x k ) u (2 p f ( x k )) v = β u,v k � �� � u,v u,v k e u,v should be small. � This is a shortest vector problem in the lattice Z e u,v ! u,v -31-

  46. How do we compute P 1 and P 2 ? We want P 1 (2 p − 1 x k , 2 p f ( x k )) and P 2 (2 p − 1 x k , 2 p f ( x k )) small, i.e. the vectors �� � � � (2 p − 1 x k ) u (2 p f ( x k )) v � α u,v (2 p − 1 x k ) u (2 p f ( x k )) v = α u,v k � �� � u,v u,v k e u,v and �� � � � (2 p − 1 x k ) u (2 p f ( x k )) v � β u,v (2 p − 1 x k ) u (2 p f ( x k )) v = β u,v k � �� � u,v u,v k e u,v should be small. � This is a shortest vector problem in the lattice Z e u,v ! u,v LLL yields two fairly short vectors, corresponding to the small values P 1 ( x k , f ( x k )) and P 2 ( x k , f ( x k )) we wish to obtain. -31-

  47. How do we compute P 1 and P 2 ? Dimension of the Lattice Let d ∈ N \ { 0 } and N = ( d + 1)( d + 2) / 2 . We introduce f kℓ ( x ) = (2 p − 1 x ) ℓ (2 p f ( x )) k for 0 � k � d, 0 � ℓ � d − k. -32-

  48. How do we compute P 1 and P 2 ? Dimension of the Lattice Let d ∈ N \ { 0 } and N = ( d + 1)( d + 2) / 2 . We introduce f kℓ ( x ) = (2 p − 1 x ) ℓ (2 p f ( x )) k for 0 � k � d, 0 � ℓ � d − k. Let P 1 = � α k,ℓ X ℓ Y k and 0 � k � d 0 � ℓ � d − k P 2 = � β k,ℓ X ℓ Y k ∈ Z [ X, Y ] . We want to have 0 � k � d 0 � ℓ � d − k � � � P k (2 p − 1 x, 2 p f ( x )) � < 1 , k = 1 , 2 , (2) for all x ∈ I = [ a, b ] . -32-

  49. How do we compute P 1 and P 2 ? Dimension of the Lattice Let d ∈ N \ { 0 } and N = ( d + 1)( d + 2) / 2 . We introduce f kℓ ( x ) = (2 p − 1 x ) ℓ (2 p f ( x )) k for 0 � k � d, 0 � ℓ � d − k. Let P 1 = � α k,ℓ X ℓ Y k and 0 � k � d 0 � ℓ � d − k P 2 = � β k,ℓ X ℓ Y k ∈ Z [ X, Y ] . We want to have 0 � k � d 0 � ℓ � d − k � � � P k (2 p − 1 x, 2 p f ( x )) � < 1 , k = 1 , 2 , (2) for all x ∈ I = [ a, b ] . To do so, discretize (2): plug (scaled) Chebyshev nodes x k . If P 1 and P 2 are small at these points, then they are small on the whole domain. -32-

  50. How do we compute P 1 and P 2 ? Volume of the Lattice Let d ∈ N \ { 0 } and N = ( d + 1)( d + 2) / 2 . We introduce f kℓ ( x ) = (2 p − 1 x ) ℓ (2 p f ( x )) k for 0 � k � d, 0 � ℓ � d − k. Let ( x j ) 0 � j � N − 1 denote Chebyshev nodes for the interval I = [ a, b ] , we want P 1 (2 p − 1 x j , 2 p f ( x j )) and P 2 (2 p − 1 x j , 2 p f ( x j )) to be small. -33-

  51. How do we compute P 1 and P 2 ? Volume of the Lattice Let d ∈ N \ { 0 } and N = ( d + 1)( d + 2) / 2 . We introduce f kℓ ( x ) = (2 p − 1 x ) ℓ (2 p f ( x )) k for 0 � k � d, 0 � ℓ � d − k. Let ( x j ) 0 � j � N − 1 denote Chebyshev nodes for the interval I = [ a, b ] , we want P 1 (2 p − 1 x j , 2 p f ( x j )) and P 2 (2 p − 1 x j , 2 p f ( x j )) to be small. It is related to the smallness of the volume of the lattice � Z ( f kℓ ( x j − 1 )) 1 � j � N i.e. 0 � k � d 0 � ℓ � d − k � � 1 /N det ( f kℓ ( x j − 1 )) 0 � k � d . 0 � ℓ � d − k, 1 � j � N -33-

  52. How do we compute P 1 and P 2 ? Volume of the Lattice Let d ∈ N \ { 0 } and N = ( d + 1)( d + 2) / 2 . We introduce f kℓ ( x ) = (2 p − 1 x ) ℓ (2 p f ( x )) k for 0 � k � d, 0 � ℓ � d − k. Let ( x j ) 0 � j � N − 1 denote Chebyshev nodes for the interval I = [ a, b ] . We have, for ρ > 1 , � � 1 /N det ( f kℓ ( x j − 1 )) 0 � k � d � 0 � ℓ � d − k, 1 � j � N � � ρ + ρ − 1 � � d/ 3 “negligible” (2 p 2 p − 1 ) d/ 3 � � b − a + b + a � � M ρ,a,b ( f ) d/ 3 . � � ρ ( N − 1) / 2 2 2 2 where M ρ,a,b ( f ) = max z ∈E ρ,a,b | f ( z ) | and � � ρe iθ + ρ − 1 e − iθ b − a + a + b . E ρ,a,b = 2 , θ ∈ [0 , 2 π ] 2 2 -34-

  53. How do we compute P 1 and P 2 ? Volume of the Lattice We have, for ρ > 1 , � � 1 /N det ( f kℓ ( x j − 1 )) 0 � k � d � 0 � ℓ � d − k, 1 � j � N � � ρ + ρ − 1 � � d/ 3 � � “negligible” (2 p 2 p − 1 ) d/ 3 b − a + b + a � � M ρ,a,b ( f ( z )) d/ 3 . � � ρ ( N − 1) / 2 2 2 2 where M ρ,a,b ( f ( z )) = max z ∈E ρ,a,b | f ( z ) | and � � ρe iθ + ρ − 1 e − iθ b − a + a + b . E ρ,a,b = 2 , θ ∈ [0 , 2 π ] 2 2 -35-

  54. How do we compute P 1 and P 2 ? Volume of the Lattice We have, for ρ > 1 , � � 1 /N det ( f kℓ ( x j − 1 )) 0 � k � d � 0 � ℓ � d − k, 1 � j � N � � ρ + ρ − 1 � � d/ 3 � � “negligible” (2 p 2 p − 1 ) d/ 3 b − a + b + a � � M ρ,a,b ( f ( z )) d/ 3 . � � ρ ( N − 1) / 2 2 2 2 where M ρ,a,b ( f ( z )) = max z ∈E ρ,a,b | f ( z ) | and � � ρe iθ + ρ − 1 e − iθ b − a + a + b . E ρ,a,b = 2 , θ ∈ [0 , 2 π ] 2 2 For Gamma, d = O ( p ) is enough to tackle the whole [ a, b ] -35-

  55. How do we compute P 1 and P 2 ? Volume of the Lattice We have, for ρ > 1 , � � 1 /N det ( f kℓ ( x j − 1 )) 0 � k � d � 0 � ℓ � d − k, 1 � j � N � � ρ + ρ − 1 � � d/ 3 � � “negligible” (2 p 2 p − 1 ) d/ 3 b − a + b + a � � M ρ,a,b ( f ( z )) d/ 3 . � � ρ ( N − 1) / 2 2 2 2 where M ρ,a,b ( f ( z )) = max z ∈E ρ,a,b | f ( z ) | and � � ρe iθ + ρ − 1 e − iθ b − a + a + b . E ρ,a,b = 2 , θ ∈ [0 , 2 π ] 2 2 For Gamma, d = O ( p ) is enough to tackle the whole [ a, b ] ( ρ = 2 / ( b − a ) ). -35-

  56. How do we compute P 1 and P 2 ? Volume of the Lattice We have, for ρ > 1 , � � 1 /N det ( f kℓ ( x j − 1 )) 0 � k � d � 0 � ℓ � d − k, 1 � j � N � � ρ + ρ − 1 � � d/ 3 � � “negligible” (2 p 2 p − 1 ) d/ 3 b − a + b + a � � M ρ,a,b ( f ( z )) d/ 3 . � � ρ ( N − 1) / 2 2 2 2 where M ρ,a,b ( f ( z )) = max z ∈E ρ,a,b | f ( z ) | and � � ρe iθ + ρ − 1 e − iθ b − a + a + b . E ρ,a,b = 2 , θ ∈ [0 , 2 π ] 2 2 For Gamma, d = O ( p ) is enough to tackle the whole [ a, b ] ( ρ = 2 / ( b − a ) ). If [ a, b ] = [1 , 2] , less than one hour for p = 53 and 1 CPU year for p = 113 . -35-

  57. The Table Maker’s Dilemma A breakpoint is a point where the rounding function changes. In this talk, it is the middle of two consecutive FP numbers. Two-step challenge: Determine the set BP f of all the FP numbers x such that f ( x ) is a breakpoint; -36-

  58. The Table Maker’s Dilemma A breakpoint is a point where the rounding function changes. In this talk, it is the middle of two consecutive FP numbers. Two-step challenge: Determine the set BP f of all the FP numbers x such that f ( x ) is a breakpoint; Find m ∈ N , as small as possible, s.t. for all FP number y / ∈ BP f , if we use an internal precision of m bits to evaluate the f ( y ) ’s, then we will always get RN ( f ( y )) . -36-

  59. The Table Maker’s Dilemma A breakpoint is a point where the rounding function changes. In this talk, it is the middle of two consecutive FP numbers. Two-step challenge: Determine the set BP f of all the FP numbers x such that f ( x ) is a breakpoint; Find m ∈ N , as small as possible, s.t. for all FP number y / ∈ BP f , if we use an internal precision of m bits to evaluate the f ( y ) ’s, then we will always get RN ( f ( y )) . In other words, find m ∈ N , as small as possible, such that for all j, ℓ ∈ � 2 p − 1 , 2 p − 1 � s.t. j/ 2 p − 1 / ∈ BP f and � � � � � � j − 2 ℓ + 1 � � � � 2 − m . � f 2 p − 1 2 p -36-

  60. Finding m beyond which there is no problem ? Function f : sin , cos , arcsin , arccos , tan , arctan , exp , log , sinh , cosh . -37-

  61. Finding m beyond which there is no problem ? Function f : sin , cos , arcsin , arccos , tan , arctan , exp , log , sinh , cosh . Hermite-Lindemann’s theorem ( α � = 0 algebraic ⇒ e α transcendental). If x is a FP number, there exists an m x , s.t. rounding the m x -bit approximation ⇔ rounding f ( x ) ; -37-

  62. Finding m beyond which there is no problem ? Function f : sin , cos , arcsin , arccos , tan , arctan , exp , log , sinh , cosh . Hermite-Lindemann’s theorem ( α � = 0 algebraic ⇒ e α transcendental). If x is a FP number, there exists an m x , s.t. rounding the m x -bit approximation ⇔ rounding f ( x ) ; finite number of FP numbers → ∃ m max = max x ( m x ) s.t. ∀ x , rounding the m max -bit approximation to f ( x ) is equivalent to rounding f ( x ) ; -37-

  63. Finding m beyond which there is no problem ? Function f : sin , cos , arcsin , arccos , tan , arctan , exp , log , sinh , cosh . Hermite-Lindemann’s theorem ( α � = 0 algebraic ⇒ e α transcendental). If x is a FP number, there exists an m x , s.t. rounding the m x -bit approximation ⇔ rounding f ( x ) ; finite number of FP numbers → ∃ m max = max x ( m x ) s.t. ∀ x , rounding the m max -bit approximation to f ( x ) is equivalent to rounding f ( x ) ; Questions: how large can m be? How to determine it? -37-

  64. TMD: What can Diophantine Approximation do for us? Algebraic functions. J.-M. Muller (CNRS, ENS Lyon) used an old and simple (but fundamental) result by J. Liouville to give non trivial √ estimates for : roughly, m � dp . d The exp function. Khemira and Voutier (2011): for binary64 ( p = 53 ) format, m � 16 , 490 -38-

  65. TMD: What can Diophantine Approximation do for us? Algebraic functions. J.-M. Muller (CNRS, ENS Lyon) used an old and simple (but fundamental) result by J. Liouville to give non trivial √ estimates for : roughly, m � dp . d The exp function. Khemira and Voutier (2011): for binary64 ( p = 53 ) format, m � 16 , 490( ∼ 311 p ) , whereas the hand-waving argument suggests that precision 120 should be sufficient, binary128 ( p = 113 ) format, m � 45 , 928 -38-

  66. TMD: What can Diophantine Approximation do for us? Algebraic functions. J.-M. Muller (CNRS, ENS Lyon) used an old and simple (but fundamental) result by J. Liouville to give non trivial √ estimates for : roughly, m � dp . d The exp function. Khemira and Voutier (2011): for binary64 ( p = 53 ) format, m � 16 , 490( ∼ 311 p ) , whereas the hand-waving argument suggests that precision 120 should be sufficient, binary128 ( p = 113 ) format, m � 45 , 928( ∼ 406 p ) , whereas the hand-waving argument suggests that precision 250 should be sufficient. -38-

  67. TMD: What can Diophantine Approximation do for us? Algebraic functions. J.-M. Muller (CNRS, ENS Lyon) used an old and simple (but fundamental) result by J. Liouville to give non trivial √ estimates for : roughly, m � dp . d The exp function. Khemira and Voutier (2011): for binary64 ( p = 53 ) format, m � 16 , 490( ∼ 311 p ) , whereas the hand-waving argument suggests that precision 120 should be sufficient, binary128 ( p = 113 ) format, m � 45 , 928( ∼ 406 p ) , whereas the hand-waving argument suggests that precision 250 should be sufficient. Need for an algorithmic approach! -38-

  68. TMD: Algorithmic Approaches Reminder: p = 24 , 53 , 64 , 113 . � � � � � � j − 2 ℓ + 1 1 Find all 2 p − 1 � j, ℓ � 2 p − 1 s.t. � � 2 2 p . � f � < 2 p − 1 2 p Two algorithmic approaches: Lefèvre (extension of Euclid’s algorithm) in O (2 2 p/ 3 ) : TMD solved for p = 53 . Stehlé, Lefèvre & Zimmermann (Coppersmith method, based on LLL) in O (2 p/ 2 ) : TMD solved for p = 64 . -39-

  69. TMD: Algorithmic Approaches Reminder: p = 24 , 53 , 64 , 113 . � � � � � � j − 2 ℓ + 1 1 Find all 2 p − 1 � j, ℓ � 2 p − 1 s.t. � � 2 2 p . � f � < 2 p − 1 2 p Two algorithmic approaches: Lefèvre (extension of Euclid’s algorithm) in O (2 2 p/ 3 ) : TMD solved for p = 53 . Stehlé, Lefèvre & Zimmermann (Coppersmith method, based on LLL) in O (2 p/ 2 ) : TMD solved for p = 64 . We need another approach to tackle the case of the quadruple precision ( p = 113 )! -39-

  70. Relaxed TMD � � � � � � j − 2 ℓ + 1 � < 1 2 m , 2 p − 1 � j, ℓ � 2 p − 1 , � � What about solving � f 2 p − 1 2 p when m = 6 p, 8 p, 10 p instead of 2 p ? -40-

  71. Relaxed TMD � � � � � � j − 2 ℓ + 1 � < 1 2 m , 2 p − 1 � j, ℓ � 2 p − 1 , � � What about solving � f 2 p − 1 2 p when m = 6 p, 8 p, 10 p instead of 2 p ? Lefèvre’s algorithm and Stehlé’s improvement of SLZ algorithm get unpractical for p = 113 (quadruple precision). -40-

Recommend


More recommend