Reliable multiprecision implementation of a class of special - PowerPoint PPT Presentation

Reliable multiprecision implementation of a class of special functions Team: A. Cuyt, V.B. Petersen, B. Verdonk, H. Waadeland, F. Backeljauw, S. Becuwe, M. Colman Universiteit Antwerpen Sør-Trøndelag University College Norwegian University of Science and Technology 1 / 34

Special functions ⊲ Special functions Type I Type II Type III Building blocks ◮ Incomplete gamma and related Implementation ◮ Error function and related ◮ Exponential integrals ◮ Hypergeometric 2 F 1 ( a , n ; c ; x ) , n ∈ Z ◮ Confluent hypergeometric 1 F 1 ( n ; b ; x ) , n ∈ Z ◮ Legendre functions ◮ Bessel functions of integer and fractional order ◮ Coulomb wave functions 2 / 34

Special functions Special functions ⊲ Type I Type II Type III Building blocks ◮ Incomplete gamma and related Implementation ◮ Error function and related ◮ Exponential integrals ◮ Hypergeometric 2 F 1 ( a , n ; c ; x ) , n ∈ Z ◮ Confluent hypergeometric 1 F 1 ( n ; b ; x ) , n ∈ Z ◮ Legendre functions ◮ Bessel functions of integer and fractional order ◮ Coulomb wave functions 3 / 34

Type I Special functions ⊲ Type I Type II Type III Building blocks monotone behaviour of a m ( x ) , b m ( x ) Implementation ∞ a m ( x ) K f ( x ) = b 0 ( x ) + b m ( x ) m = 1 a 1 ( x ) = b 0 ( x ) + a 2 ( x ) b 1 ( x ) + a 3 ( x ) b 2 ( x ) + b 3 ( x ) + · · · 4 / 34

Example: Special functions ⊲ Type I Type II γ ( a , x ) = x a e − x ∞ ( m − 1 ) x Type III K Building blocks a − x + a + ( m − 1 ) − x m = 2 Implementation x � = 0, a > 0 a 1 ( x ) = x a e − x b 1 ( x ) = a − x a m ( x ) = ( m − 1 ) x b m ( x ) = a + ( m − 1 ) − x a m ( x ) / ( b m ( x ) b m − 1 ( x )) ց 0 5 / 34

Special functions Special functions Type I ⊲ Type II Type III Building blocks ◮ Incomplete gamma and related Implementation ◮ Error function and related ◮ Exponential integrals ◮ Hypergeometric 2 F 1 ( a , n ; c ; x ) , n ∈ Z ◮ Confluent hypergeometric 1 F 1 ( n ; b ; x ) , n ∈ Z ◮ Legendre functions ◮ Bessel functions of integer and fractional order ◮ Coulomb wave functions 6 / 34

Type II Special functions Type I ⊲ Type II Type III Building blocks Implementation bound a ( k ) m ( x ) , b ( k ) m ( x ) , for m � n k , 1 � k � l l � f ( x ) = f k ( x ) k = 1 a ( k ) ∞ m ( x ) K f k ( x ) = b 0 ( x ) + b ( k ) m ( x ) m = 1 7 / 34

Special functions Type I Example: ⊲ Type II Type III Building blocks 2 F 1 ( a , b ; c ; x ) ∞ a m x Implementation K f k ( x ) = 2 F 1 ( a , b + 1; c + 1; x ) = 1 + 1 m = 1 x < 1 a 2 m + 1 = −( a + m )( c − b + m ) ( c + 2 m )( c + 2 m + 1 ) , m � 0 a 2 m = −( b + m )( c − a + m ) ( c + 2 m − 1 )( c + 2 m ) , m � 1 8 / 34

Special functions Special functions Type I Type II ⊲ Type III Building blocks ◮ Incomplete gamma and related Implementation ◮ Error function and related ◮ Exponential integrals ◮ Hypergeometric 2 F 1 ( a , n ; c ; x ) , n ∈ Z ◮ Confluent hypergeometric 1 F 1 ( n ; b ; x ) , n ∈ Z ◮ Legendre functions ◮ Bessel functions of integer and fractional order ◮ Coulomb wave functions 9 / 34

Type III Special functions Type I Type II ⊲ Type III Building blocks each f k ( x ) either type I or type II Implementation l � f ( x ) = f 0 ( x ) f k ( x ) k = 1 a ( k ) ∞ m ( x ) K f k ( x ) = b 0 ( x ) + b ( k ) m ( x ) m = 1 f 0 ( x ) obtained separately 10 / 34

Example: Special functions Type I Type II ⊲ Type III J 5 / 2 ( x ) = J 1 / 2 ( x ) J 3 / 2 ( x ) J 5 / 2 ( x ) Building blocks J 1 / 2 ( x ) J 3 / 2 ( x ) Implementation � 2 J 1 / 2 ( x ) = πx sin x − x 2 / ( 4 ( ν + m − 1 )( ν + m )) = x/ ( 2 ν + 2 ) ∞ J ν + 1 ( x ) K J ν ( x ) 1 + 1 m = 2 x � 0, v � 0 a m ( x ) ր 0 11 / 34

Building blocks Special functions ⊲ Building blocks Equivalence and notation Interval Sequence Theorem Reliability Implementation ◮ equivalence transformation ◮ continued fraction approximant ◮ continued fraction tail ◮ truncation error bound ◮ roundoff error bound 12 / 34

Equivalence and notation Special functions Building blocks ⊲ Equivalence and notation Interval Sequence Theorem Reliability b m ( x ) = 1 Implementation non-terminating fraction: a m ( x ) � = 0, b m ( x ) � = 0 a 1 ( x ) a 1 ( x ) /b 1 ( x ) = a 2 ( x ) a 2 ( x ) /b 1 ( x ) b 2 ( x ) b 1 ( x ) + 1 + a 3 ( x ) 1 + a 3 ( x ) /b 2 ( x ) b 3 ( x ) b 2 ( x ) + b 3 ( x ) + · · · 1 + · · · 13 / 34

Equivalence and notation Special functions Building blocks approximant ⊲ Equivalence and notation Interval Sequence Theorem Reliability N − 1 a m ( x ) a N ( x ) Implementation K f N ( x ; w N ) = 1 + 1 + w N m = 1 N → ∞ f N ( x ; w N ) = f ( x ) lim tail ∞ a m ( x ) f ( N ) = K 1 m = N + 1 N → ∞ f ( N ) ? lim = 0 14 / 34

Example: Special functions 1 = − 1 + √ 1 + 4 x ∞ x Building blocks K ⊲ Equivalence and notation 2 Interval Sequence m = 1 f ( N ) = − 1 + √ 1 + 4 x Theorem Reliability Implementation 2 ∞ m ( m + 2 ) K 1 = 1 m = 1 f ( N ) = N + 1 √ 2 − 1 = 1 2 1 2 1 1 + 1 + 1 + 1 + 1 + · · · √ f ( 2 N ) → 2 − 1 √ f ( 2 N + 1 ) → 2 15 / 34

Interval Sequence Theorem Special functions Building blocks sequence of element sets { E n } n ∈ N Equivalence and notation ⊲ Interval sequence of value sets { V n } n ∈ N Sequence Theorem Reliability a n Implementation a n ∈ E n ⇒ ⊆ V n − 1 , n � 1 1 + V n f ( N ) ∈ ¯ V N , f ( x ) ∈ ¯ V 0 w N ∈ V N ⇒ f N ( x ; w N ) ∈ V 0 Example: 2 1 1 2 1 1 2 1 + 1 − 1 + 1 + 1 − 1 + 1 + · · · E 3 n = { − 1 } , E 3 n + 1 = { 2 } , E 3 n + 2 = { 1 } V 3 n = { 1 / 2 } , V 3 n + 1 = { 3 } , V 3 n + 2 = { − 2 / 3 } 16 / 34

Special functions Building blocks Equivalence and ◮ element sets notation ⊲ Interval Sequence Theorem a m ( x ) ∈ E m = [ α ( l ) m , α ( u ) m ] Reliability Implementation α ( l ) m α ( u ) m � 0 ◮ value sets V k = [ L k , R k ] f ( k ) , f ( k ) L k = min R k = max a m ∈ E m a m ∈ E m m � k + 1 m � k + 1 17 / 34

Special functions Building blocks Equivalence and notation ⊲ Interval Sequence Theorem Reliability M k := max {| L k / ( 1 + L k ) | , | R k / ( 1 + R k ) |} Implementation w N ∈ V N = [ L N , R N ] N − 1 � � f ( x ) − f N ( x ; w N ) � � R N − L N � � � M k � � f ( x ) 1 + L N � k = 1 18 / 34

Reliability Special functions Building blocks Equivalence and notation Interval Sequence Theorem ⊲ Reliability � � f ( x ) − f N ( x ; w N ) Implementation � � Relative truncation error = � � ǫ T � � f ( x ) � � � f N ( x ; w N ) − F N ( x ; w N ) � � Relative rounding error = � � ǫ R � � f ( x ) � ǫ := ǫ T + ǫ R f ( x ) ∈ [ F N ( x ; w N )( 1 − 2 ǫ ) , F N ( x ; w N )( 1 + 2 ǫ )] 19 / 34

Special functions Building blocks Equivalence and notation Example: Interval Sequence Theorem ⊲ Reliability Implementation ( m − 1 ) x a f ( a , x ) = aγ ( a , x ) e x ∞ ( a + m − 1 − x )( a + m − 2 − x ) a − x K = x a 1 + 1 m = 2 Evaluate f ( 9 / 2, 1 ) with ǫ T + ǫ R � 10 − d + 1 , d = 73, 74, . . . , 80 20 / 34

Special functions ◮ Continued fraction library output Building blocks Equivalence and notation 73 1.214009591773512617777498734645198390079596056622283491877162409691879700 Interval Sequence Theorem 74 1.2140095917735126177774987346451983900795960566222834918771624096918797000 ⊲ Reliability 75 1.21400959177351261777749873464519839007959605662228349187716240969187969998 76 1.214009591773512617777498734645198390079596056622283491877162409691879699983 Implementation 77 1.2140095917735126177774987346451983900795960566222834918771624096918796999829 78 1.21400959177351261777749873464519839007959605662228349187716240969187969998292 79 1.214009591773512617777498734645198390079596056622283491877162409691879699982919 80 1.2140095917735126177774987346451983900795960566222834918771624096918796999829190 ◮ Maple output 73 1.214009591773512617777498734645198390079596056622283491877162409691879774 74 1.2140095917735126177774987346451983900795960566222834918771624096918797015 75 1.21400959177351261777749873464519839007959605662228349187716240969187969966 76 1.214009591773512617777498734645198390079596056622283491877162409691879700001 77 1.2140095917735126177774987346451983900795960566222834918771624096918796999764 78 1.21400959177351261777749873464519839007959605662228349187716240969187969998223 79 1.214009591773512617777498734645198390079596056622283491877162409691879699982930 80 1.2140095917735126177774987346451983900795960566222834918771624096918796999829239 21 / 34

Implementation Special functions Building blocks ⊲ Implementation Book Demo ◮ Book with tables and graphs ◮ C ++ library (all functions ) ◮ Maple library (all fractions ) ◮ Web version to explore both 22 / 34

Formulas Special functions Building blocks Implementation ⊲ Book Demo 23 / 34

Tables Special functions Building blocks Implementation ⊲ Book Demo 24 / 34

Special functions Building blocks Implementation ⊲ Book Demo 25 / 34

Figures Special functions Building blocks Implementation ⊲ Book Demo 26 / 34

Special functions Building blocks Implementation ⊲ Book Demo 27 / 34