Computation of the error functions erf and erfc in arbitrary - PowerPoint PPT Presentation

Computation of the error functions erf and erfc in arbitrary precision with correct rounding Sylvain Chevillard Arenaire, LIP, ENS-Lyon, France Sylvain.Chevillard@ens-lyon.fr Nathalie Revol INRIA, Arenaire, LIP, ENS-Lyon, France Nathalie.Revol@ens-lyon.fr ICMS, Castro Urdiales, Spain 1-3 September 2006

IEEE-754 standard for floating-point arithmetic Floating-point number : s.m.β e or 0 , a denormal, ±∞ , NaN. • s : sign • m : mantissa ∈ [2 p − 1 / 2 p , (2 p − 1) / 2 p ] , p is the precision • e : exponent ∈ [ e min , e max ] • β : basis, usually equal to 2 IEEE-754 standard for floating-point arithmetic : • fixed formats for single and double precision • specifies 4 rounding modes • specifies the arithmetic and algebraic operations + , − , × , / , √ . Advantages : • well-specified arithmetic, reproducible results • error bounds can be established and proofs can be done ICMS 2006 1 S. Chevillard & N. Revol

Desirable extensions of the IEEE-754 standard correctly rounded elementary functions : cf. current revision of the standard arbitrary precision : (software) cf. MPFR correctly rounded special functions : ? ? ? correctly rounded functions in arbitrary precision : cf. MPFR for elementary functions ? ? ? ICMS 2006 2 S. Chevillard & N. Revol

How can one return the correctly rounded evaluation of a function f ? To return the correctly rounded evaluation of f ( x ) in precision p : 1- approximate f ( x ) with larger precision q , error < ε 2- if can round ( f ( x ) , ε, p ) f(x) then return it 3- otherwise increase q , decrease ε and try again. ICMS 2006 3 S. Chevillard & N. Revol

Outline of this talk • the error function erf • algorithm to return the correctly rounded evaluation of erf ( x ) • experimental results • the complementary error function erfc • conclusion and hints for improvements • current work on interval arithmetic and algorithms : MPFI ICMS 2006 4 S. Chevillard & N. Revol

The error function erf 1 0.5 –3 –2 –1 1 2 3 x –0.5 � x 0 e − t 2 dt 2 erf ( x ) = √ π –1 The error function is very useful in statistics (cf. Gaussian distribution), diffusion equation (special configurations) and other heat transfer pro- blems. . . Goal : return the correctly rounded value of erf ( x ) in arbitrary precision. ICMS 2006 5 S. Chevillard & N. Revol

Possible formulas � x 0 e − t 2 dt 2 erf ( x ) = (1) √ π ( − 1) n x 2 n +1 � + ∞ 2 = (2) √ π n =0 n !(2 n +1) 2 n x 2 n +1 √ π e − x 2 � + ∞ 2 = (3) n =0 1 . 3 ··· (2 n +1) √ 2 � + ∞ n =0 ( − 1) n � I 2 n +1 / 2 ( x 2 ) − I 2 n +3 / 2 ( x 2 ) � = (4) e − x 2 1 / 2 3 / 2 1 1 2 = x + . . . (5) √ π x + x + x + x + ICMS 2006 6 S. Chevillard & N. Revol

Possible formulas � x 0 e − t 2 dt 2 erf ( x ) = (1) √ π ( − 1) n x 2 n +1 � + ∞ 2 = (2) √ π n =0 n !(2 n +1) 2 n x 2 n +1 √ π e − x 2 � + ∞ 2 = (3) n =0 1 . 3 ··· (2 n +1) √ 2 � + ∞ n =0 ( − 1) n � I 2 n +1 / 2 ( x 2 ) − I 2 n +3 / 2 ( x 2 ) � = (4) e − x 2 1 / 2 3 / 2 1 1 2 = x + . . . (5) √ π x + x + x + x + Discussion of the use of quadrature : • numerous evaluations of exp : costly (either in computing time or in development time) • many sources of error : evaluation of exp , quadrature, roundoff ICMS 2006 7 S. Chevillard & N. Revol

Possible formulas � x 0 e − t 2 dt 2 erf ( x ) = (1) √ π ( − 1) n x 2 n +1 � + ∞ 2 = (2) √ π n =0 n !(2 n +1) 2 n x 2 n +1 √ π e − x 2 � + ∞ 2 = (3) n =0 1 . 3 ··· (2 n +1) √ 2 � + ∞ n =0 ( − 1) n � I 2 n +1 / 2 ( x 2 ) − I 2 n +3 / 2 ( x 2 ) � = (4) e − x 2 1 / 2 3 / 2 1 1 2 = x + . . . (5) √ π x + x + x + x + Discussion of the use of alternate power series : • the remainder is easy to bound • sum of terms of alternate signs : cancellation ICMS 2006 8 S. Chevillard & N. Revol

Possible formulas � x 0 e − t 2 dt 2 erf ( x ) = (1) √ π ( − 1) n x 2 n +1 � + ∞ 2 = (2) √ π n =0 n !(2 n +1) 2 n x 2 n +1 √ π e − x 2 � + ∞ 2 = (3) n =0 1 . 3 ··· (2 n +1) √ 2 � + ∞ n =0 ( − 1) n � I 2 n +1 / 2 ( x 2 ) − I 2 n +3 / 2 ( x 2 ) � = (4) e − x 2 1 / 2 3 / 2 1 1 2 = x + . . . (5) √ π x + x + x + x + Discussion of the use of this power series : • sum of positive terms : numerical stability • the remainder is less easy to bound ICMS 2006 9 S. Chevillard & N. Revol

Possible formulas � x 0 e − t 2 dt 2 erf ( x ) = (1) √ π ( − 1) n x 2 n +1 � + ∞ 2 = (2) √ π n =0 n !(2 n +1) 2 n x 2 n +1 √ π e − x 2 � + ∞ 2 = (3) n =0 1 . 3 ··· (2 n +1) √ 2 � + ∞ n =0 ( − 1) n � I 2 n +1 / 2 ( x 2 ) − I 2 n +3 / 2 ( x 2 ) � = (4) e − x 2 1 / 2 3 / 2 1 1 2 = x + . . . (5) √ π x + x + x + x + Discussion of the use of Bessel functions of fractional order : • the problem is now to evaluate the Bessel functions of fractional order ICMS 2006 10 S. Chevillard & N. Revol

Possible formulas � x 0 e − t 2 dt 2 erf ( x ) = (1) √ π ( − 1) n x 2 n +1 � + ∞ 2 = (2) √ π n =0 n !(2 n +1) 2 n x 2 n +1 √ π e − x 2 � + ∞ 2 = (3) n =0 1 . 3 ··· (2 n +1) √ 2 � + ∞ n =0 ( − 1) n � I 2 n +1 / 2 ( x 2 ) − I 2 n +3 / 2 ( x 2 ) � = (4) e − x 2 1 / 2 3 / 2 1 1 2 = x + . . . (5) √ π x + x + x + x + Discussion of the use of this continued fraction : • the remainder is less easy to bound • numerical stability is not ensured ICMS 2006 11 S. Chevillard & N. Revol

Chosen formula � x 0 e − t 2 dt 2 erf ( x ) = (1) √ π ( − 1) n x 2 n +1 � + ∞ 2 = (2) √ π n =0 n !(2 n +1) 2 n x 2 n +1 √ π e − x 2 � + ∞ 2 = (3) n =0 1 . 3 ··· (2 n +1) √ 2 � + ∞ n =0 ( − 1) n � I 2 n +1 / 2 ( x 2 ) − I 2 n +3 / 2 ( x 2 ) � = (4) e − x 2 1 / 2 3 / 2 1 1 2 = x + . . . (5) √ π x + x + x + x + Motivation for the choice of this alternate power series : • the remainder is easy to bound • special care to avoid cancellation in the sum of terms of alternate signs ICMS 2006 12 S. Chevillard & N. Revol

Other useful formulas erf ( − x ) = − erf ( x ) No argument reduction possible (cf. sin( x + 2 π ) = sin x or exp(2 x ) = (exp x ) 2 ). e − x 2 e − x 2 2 < erfc ( x ) ≤ 2 √ π · √ √ π · for x ≥ 0 . x 2 + 2 � x + x 2 + 4 x + π ICMS 2006 13 S. Chevillard & N. Revol

Outline • the error function erf • algorithm to return the correctly rounded evaluation of erf ( x ) • experimental results • the complementary error function erfc • conclusion and hints for improvements • current work on interval arithmetic and algorithms : MPFI ICMS 2006 14 S. Chevillard & N. Revol

Algorithm Input : x, p Output : correctly rounded value of erf ( x ) with p bits 1. handle special cases : x = 0 , x = ±∞ , x < 0 2. check whether the last enclosure gives rapidly the answer 3. determine the truncation rank N needed to have an absolute error ≤ ε 4. determine the computing precision q to have roundoff error ≤ ε 5. evaluate y that approximates erf ( x ) using the alternate power series ; bound from above the roundoff error ε ′ ≤ ε , on the fly 6. if can round ( y, ε + ε ′ , p ) then return y else increase N and q and do steps (5) and (6) again ICMS 2006 15 S. Chevillard & N. Revol

Algorithm : step (2) Input : x, p Output : correctly rounded value of erf ( x ) with p bits 1. handle special cases : x = 0 , x = ±∞ , x < 0 2. check whether the last enclosure gives rapidly the answer 3. determine the truncation rank N needed to have an absolute error ≤ ε 4. determine the computing precision q to have roundoff error ≤ ε 5. evaluate y that approximates erf ( x ) using the alternate power series ; bound from above the roundoff error ε ′ ≤ ε , on the fly 6. if can round ( y, ε + ε ′ , p ) then return y else increase N and q and do steps (5) and (6) again ICMS 2006 16 S. Chevillard & N. Revol

Algorithm : step (2) Reminder : e − x 2 e − x 2 2 < erfc ( x ) ≤ 2 √ π · √ √ π · for x ≥ 0 . x 2 + 2 � x + x 2 + 4 x + π Step (2) : compute both sides of this enclosure : y L , y R if can round ( y L , y R − y L , p ) then return it. ICMS 2006 17 S. Chevillard & N. Revol

Algorithm : step (3) Input : x, p Output : correctly rounded value of erf ( x ) with p bits 1. handle special cases : x = 0 , x = ±∞ , x < 0 2. check whether the last enclosure gives rapidly the answer 3. determine the truncation rank N needed to have an absolute error ≤ ε 4. determine the computing precision q to have roundoff error ≤ ε 5. evaluate y that approximates erf ( x ) using the alternate power series ; bound from above the roundoff error ε ′ ≤ ε , on the fly 6. if can round ( y, ε + ε ′ , p ) then return y else increase N and q and do steps (5) and (6) again ICMS 2006 18 S. Chevillard & N. Revol

Algorithm : step (3) Reminder : power series + ∞ + ∞ ( − 1) n x 2 n +1 ( − 1) n x 2 n +1 erf ( x ) = 2 a n with a n = 2 � � √ π n !(2 n + 1) = √ π n !(2 n + 1) n =0 n =0 Property : alternate power series � + ∞ n =0 a n with non-increasing term a n (for n large enough) ⇒ remainder | � + ∞ n = N a n | = | erf ( x ) − � N − 1 n =0 a n | ≤ | a N | . Step (3) : ε = 2 − p − 8 : 8 extra bits evaluate a n until a N ≤ ε : N is the truncation rank. ICMS 2006 19 S. Chevillard & N. Revol