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-

Download Presentation
Download Policy: The content available on the website is offered to you 'AS IS' for your personal information and use only. It cannot be commercialized, licensed, or distributed on other websites without prior consent from the author. To download a presentation, simply click this link. If you encounter any difficulties during the download process, it's possible that the publisher has removed the file from their server.

Recommend


More recommend