Semi-Automatic Floating-Point Implementation of Special Functions - PowerPoint PPT Presentation

Semi-Automatic Floating-Point Implementation of Special Functions Christoph Lauter 1 Marc Mezzarobba 1,2 Pequan group 1 Université Paris 6 2 CNRS ARITH 22, Lyon, 2015-06-23

Code Generation for Mathematical Functions expression explicit function equation black box f = e sin ( x ) f = log ( x ) Φ ( f ) = 0 x �→ f ( x ) }main() ❞♦✉❜❧❡ ❢✉♥✭❞♦✉❜❧❡ ①✮ ④ { int temp; float celsius; char repeat; char flag; do { flag='n"; do { if(flag=='n') printf("Input a valid temperature :"); else printf("input a valid temperature,stupid:"); scanf("%d",&temp); flag='y'; } while (temp<0||temp >100); ✳✳✳ celsius=(5.0/9.0)*(temp-32); printf("%d degrees F is %6.2f degrees celsius\n",temp,cels printf("Do you have another temperature?"); repeat=getchar(); putchar('\n'); ⑥

Code Generation for Mathematical Functions expression explicit function equation black box f = e sin ( x ) f = log ( x ) Φ ( f ) = 0 x �→ f ( x ) }main() ❞♦✉❜❧❡ ❢✉♥✭❞♦✉❜❧❡ ①✮ ④ { int temp; float celsius; char repeat; char flag; do { flag='n"; do { if(flag=='n') printf("Input a valid temperature :"); else printf("input a valid temperature,stupid:"); scanf("%d",&temp); flag='y'; } while (temp<0||temp >100); ✳✳✳ celsius=(5.0/9.0)*(temp-32); printf("%d degrees F is %6.2f degrees celsius\n",temp,cels printf("Do you have another temperature?"); repeat=getchar(); putchar('\n'); ⑥

Our Focus: Special Functions 2 1 0 –6 –4 –2 2 4 6 –1 Ai ( x ) D-finite Functions –2 = sol ns of linear ODEs with polynomial coeff ts 2 1 p r ( x ) · f ( r ) ( x ) + · · · + p 1 ( x ) · f ′ ( x ) + p 0 ( x ) · f ( x ) = 0 0 –6 –4 –2 2 4 6 –1 Ci ( x ) p 0 , p 1 , . . . , p r ∈ R [ x ] –2 2 1 Example: f ( x ) = arctan ( x ) 0 –6 –4 –2 2 4 6 –1 erf ( x ) –2 � f ( 0 ) � � 0 � ( x 2 + 1 ) · f ′′ ( x ) + 2 x · f ′ ( x ) = 0, = 2 f ′ ( 0 ) 1 1 0 –4 –2 2 4 –6 6 –1 J 0 ( x ) –2

Problem Statement }main() { int temp; float celsius; char repeat; char flag; do { flag='n"; do { if(flag=='n') printf("Input a valid temperature :"); else printf("input a valid temperature,stupid:"); scanf("%d",&temp); flag='y'; } while (temp<0||temp >100); celsius=(5.0/9.0)*(temp-32); printf("%d degrees F is %6.2f degrees celsius\n",temp,cels printf("Do you have another temperature?"); repeat=getchar(); putchar('\n'); � p r f ( r ) + · · · + p 0 f = 0, ❞♦✉❜❧❡ ❢✉♥✭❞♦✉❜❧❡ ①✮❀ f ( 0 ) , . . . , f ( r − 1 ) ( 0 ) ∀ x ∈ D ∩ ❞♦✉❜❧❡ , � � ❢✉♥ ( x ) − f ( x ) D ⊆ R , � � ε > 0 � ≤ ε � � f ( x ) � ◮ “Semi-automatic”: may require implementation hints Future goal: full automation ◮ Rigorous in principle, some shortcuts in current prototype

Overview differential equation Frontend sets of polynomial approximations truncated Taylor series ◮ NumGfun (Maple) ◮ ≈ rigorous multiple compact rough approx. precision solver for ODE Backend tight approximation ◮ Sollya ◮ ≈ Metalibm-lutetia FP representable poly. (cf. previous talk by O. Kupriianova) evaluation scheme code

The Frontend: Taylor Series differential equation Frontend p r ( x ) · f ( r ) ( x ) + · · · p 0 ( x ) · f ( x ) = 0 truncated Taylor series ∞ ∑ c n x n Ansatz: f ( x ) = n = 0 compact rough approx. ∞ ◮ f ′ ( x ) = ( n + 1 ) c n + 1 x n ∑ Backend n = 0 ∞ c n − 1 x n tight approximation ∑ ◮ x · f ( x ) = n = 1 c n = b 0 ( n ) c n − 1 + b 1 ( n ) c n − 2 + · · · FP representable poly. evaluation scheme f ( x ) = c 0 + c 1 x + · · · + c d − 1 x d − 1 + · · · code

The Frontend: Error Bounds differential equation d − 1 ∞ Frontend c n x n + ∑ ∑ c n x n f ( x ) = n = 0 n = d truncated Taylor series � �� � ? | x | ≤ ρ compact rough approx. “Majorizing series”: Backend using the ODE, tight approximation f ( x ) = ∑ ˆ find a simple ˆ c n x n such that | c n | ≤ ˆ c n for all n . FP representable poly. d − 1 ∞ � c n x n � ! c n ρ n evaluation scheme � ∑ � ∑ � f ( x ) − � ≤ ˆ ≤ ε . n = 0 n = d code

The Frontend: Analytic Continuation differential equation ◮ What if D �⊆ [ − ρ , ρ ] ? Frontend i truncated Taylor series compact rough approx. Backend tight approximation − i FP representable poly. evaluation scheme code

The Frontend: Analytic Continuation differential equation ◮ What if D �⊆ [ − ρ , ρ ] ? Frontend i truncated Taylor series compact rough approx. x 1 Backend tight approximation − i � f ( x 1 ) � � f ( 0 ) � � � 1.0000... 0.5404... FP representable poly. = f ′ ( x 1 ) f ′ ( 0 ) 0.0000... 0.7352... evaluation scheme code

The Frontend: Analytic Continuation differential equation ◮ What if D �⊆ [ − ρ , ρ ] ? Frontend i truncated Taylor series compact rough approx. x 1 Backend tight approximation − i � f ( x 1 ) � � f ( 0 ) � � � 1.0000... 0.5404... FP representable poly. = f ′ ( x 1 ) f ′ ( 0 ) 0.0000... 0.7352... evaluation scheme code

The Frontend: Analytic Continuation differential equation ◮ What if D �⊆ [ − ρ , ρ ] ? Frontend i truncated Taylor series compact rough approx. x 1 x 2 Backend tight approximation − i � f ( x 1 ) � � f ( 0 ) � � � 1.0000... 0.5404... FP representable poly. = f ′ ( x 1 ) f ′ ( 0 ) 0.0000... 0.7352... � f ( x 2 ) � � f ( x 1 ) � � 1.0000... 0.5095... � evaluation scheme = f ′ ( x 2 ) f ′ ( x 1 ) 0.0000... 0.5055... code

The Frontend: Analytic Continuation differential equation ◮ What if D �⊆ [ − ρ , ρ ] ? Frontend i truncated Taylor series compact rough approx. x 1 x 2 Backend tight approximation − i � f ( x 1 ) � � f ( 0 ) � � � 1.0000... 0.5404... FP representable poly. = f ′ ( x 1 ) f ′ ( 0 ) 0.0000... 0.7352... � f ( x 2 ) � � f ( x 1 ) � � 1.0000... 0.5095... � evaluation scheme = f ′ ( x 2 ) f ′ ( x 1 ) 0.0000... 0.5055... . . . code

The Frontend: Economization differential equation � � � ≤ ε � f ( x ) − ( c 0 + · · · + c d − 1 x d − 1 ) Have: Frontend 2 x ∈ C , | x | ≤ 1 truncated Taylor series + · · · + c 7 x 7 + c 8 x 8 + c 9 x 9 c 0 compact rough approx. − c 9 7 x 7 + c ( 1 ) c ( 1 ) + · · · + c ( 1 ) 8 x 8 2 9 T 9 ( x ) 0 Backend − c 8 c ( 2 ) + · · · + c ( 2 ) 7 x 7 tight approximation 2 8 T 8 ( x ) 0 . . FP representable poly. . | T n | ≤ 1 over [ − 1, 1 ] � � d − 1 − k x k ) � f ( x ) − ( c ( k ) + · · · + c ( k ) � � evaluation scheme � ≤ ε 0 x ∈ [ − 1, 1 ] code

The Backend: Tight Approximation differential equation Frontend ◮ Frontend: truncated Taylor series � D = D i compact rough approx. x ∈ D i | p i ( x ) − f ( x ) | ≤ ε deg p i = O ( 500 ) Backend ◮ Backend: tight approximation � D ⋆ D = i FP representable poly. � � p ⋆ i ( x ) � � x ∈ D ⋆ f ( x ) − 1 � ≤ ε � � i � evaluation scheme deg p ⋆ i = O ( 10 ) code

The Backend: Tight Approximation differential equation Frontend ◮ Backend: leveraging Metalibm truncated Taylor series ◮ Reuse domain splitting algorithm sketched in last talk compact rough approx. � D = D ⋆ i ◮ Reuse minimax implementation Backend with relative error bounds tight approximation � � p ⋆ i ( x ) � � x ∈ D ⋆ f ( x ) − 1 � ≤ ε � � i � FP representable poly. ◮ Addressing issues ◮ Pure black-box interface too slow evaluation scheme ◮ Zeros of f in domain D ? code

The Backend: FP Polynomials differential equation Frontend truncated Taylor series ◮ Implementation need: p ⋆ i ∈ F [ x ] compact rough approx. Backend ◮ Leverage existing FP minimax heuristics tight approximation FP representable poly. evaluation scheme code

The Backend: FP Polynomials differential equation Frontend truncated Taylor series ◮ Implementation need: p ⋆ i ∈ F [ x ] compact rough approx. Backend ◮ Leverage existing FP minimax heuristics tight approximation FP representable poly. ◮ unless f has a zero in the domain evaluation scheme code

The Backend: Zeros of f differential equation Frontend ◮ Suppose f ( c ) = 0 for c ∈ D truncated Taylor series ◮ We want � � p ⋆ i ( x ) � � compact rough approx. f ( x ) − 1 � ≤ ε � � � Backend ◮ We need p ⋆ i ( c ) = 0 tight approximation ◮ No if p ⋆ i ∈ F [ x ] but c ∈ R \ F FP representable poly. evaluation scheme code

The Backend: Zeros of f differential equation Frontend ◮ Suppose f ( c ) = 0 for c ∈ D truncated Taylor series ◮ Actually we want � � p ⋆ i ( x ) compact rough approx. � � ∀ x ∈ F . f ( x ) − 1 � ≤ ε � � � Backend ◮ p ⋆ i ( x ) never needs to be 0 exactly tight approximation ◮ Okay but we need to compute p ⋆ i ◮ Express c as a symbol, evaluated FP representable poly. with Newton-Raphson on f ◮ Give f ( x − c ) / x k to minimax algorithm, yielding q evaluation scheme i from x k q ( x ) ◮ Deduce p ⋆ code