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

▶
reliable multiprecision implementation of a class of
SMART_READER_LITE
LIVE PREVIEW

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

Oct 28, 2022 199 likes •540 views

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 Sr-Trndelag University College Norwegian University

slide-1
SLIDE 1

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

slide-2
SLIDE 2

⊲ Special functions

Type I Type II Type III

Building blocks Implementation

Special functions

◮ Incomplete gamma and related ◮ Error function and related ◮ Exponential integrals ◮ Hypergeometric 2F1(a, n; c; x), n ∈ Z ◮ Confluent hypergeometric 1F1(n; b; x), n ∈ Z ◮ Legendre functions ◮ Bessel functions of integer and fractional order ◮ Coulomb wave functions

2 / 34

slide-3
SLIDE 3

Special functions

⊲ Type I

Type II Type III

Building blocks Implementation

Special functions

◮ Incomplete gamma and related ◮ Error function and related ◮ Exponential integrals ◮ Hypergeometric 2F1(a, n; c; x), n ∈ Z ◮ Confluent hypergeometric 1F1(n; b; x), n ∈ Z ◮ Legendre functions ◮ Bessel functions of integer and fractional order ◮ Coulomb wave functions

3 / 34

slide-4
SLIDE 4

Special functions

⊲ Type I

Type II Type III

Building blocks Implementation

Type I

monotone behaviour of am(x), bm(x) f(x) = b0(x) +

∞

K

m=1

am(x) bm(x) = b0(x) + a1(x) b1(x) + a2(x) b2(x) + a3(x) b3(x) + · · ·

4 / 34

slide-5
SLIDE 5

Special functions

⊲ Type I

Type II Type III

Building blocks Implementation

Example: γ(a, x) = xae−x a − x +

∞

K

m=2

(m − 1)x a + (m − 1) − x x = 0, a > 0 a1(x) = xae−x b1(x) = a − x am(x) = (m − 1)x bm(x) = a + (m − 1) − x am(x)/(bm(x)bm−1(x)) ց 0

5 / 34

slide-6
SLIDE 6

Special functions

Type I

⊲ Type II

Type III

Building blocks Implementation

Special functions

◮ Incomplete gamma and related ◮ Error function and related ◮ Exponential integrals ◮ Hypergeometric 2F1(a, n; c; x), n ∈ Z ◮ Confluent hypergeometric 1F1(n; b; x), n ∈ Z ◮ Legendre functions ◮ Bessel functions of integer and fractional order ◮ Coulomb wave functions

6 / 34

slide-7
SLIDE 7

Special functions

Type I

⊲ Type II

Type III

Building blocks Implementation

Type II

bound a(k)

m (x), b(k) m (x), for m nk, 1 k l

f(x) =

l

  • k=1

fk(x) fk(x) = b0(x) +

∞

K

m=1

a(k)

m (x)

b(k)

m (x)

7 / 34

slide-8
SLIDE 8

Special functions

Type I

⊲ Type II

Type III

Building blocks Implementation

Example: fk(x) =

2F1(a, b; c; x) 2F1(a, b + 1; c + 1; x) = 1 + ∞

K

m=1

amx 1 x < 1 a2m+1 = −(a + m)(c − b + m) (c + 2m)(c + 2m + 1), m 0 a2m = −(b + m)(c − a + m) (c + 2m − 1)(c + 2m), m 1

8 / 34

slide-9
SLIDE 9

Special functions

Type I Type II

⊲ Type III

Building blocks Implementation

Special functions

◮ Incomplete gamma and related ◮ Error function and related ◮ Exponential integrals ◮ Hypergeometric 2F1(a, n; c; x), n ∈ Z ◮ Confluent hypergeometric 1F1(n; b; x), n ∈ Z ◮ Legendre functions ◮ Bessel functions of integer and fractional order ◮ Coulomb wave functions

9 / 34

slide-10
SLIDE 10

Special functions

Type I Type II

⊲ Type III

Building blocks Implementation

Type III

each fk(x) either type I or type II f(x) = f0(x)

l

  • k=1

fk(x) fk(x) = b0(x) +

∞

K

m=1

a(k)

m (x)

b(k)

m (x)

f0(x) obtained separately

10 / 34

slide-11
SLIDE 11

Special functions

Type I Type II

⊲ Type III

Building blocks Implementation

Example: J5/2(x) = J1/2(x)J3/2(x) J1/2(x) J5/2(x) J3/2(x) J1/2(x) =

  • 2

πx sin x Jν+1(x) Jν(x) = x/(2ν + 2) 1 +

∞

K

m=2

−x2/(4(ν + m − 1)(ν + m)) 1 x 0, v 0 am(x) ր 0

11 / 34

slide-12
SLIDE 12

Special functions ⊲ Building blocks

Equivalence and notation Interval Sequence Theorem Reliability

Implementation

Building blocks

◮ equivalence transformation ◮ continued fraction approximant ◮ continued fraction tail ◮ truncation error bound ◮ roundoff error bound

12 / 34

slide-13
SLIDE 13

Special functions Building blocks

⊲ Equivalence and

notation Interval Sequence Theorem Reliability

Implementation

Equivalence and notation

bm(x) = 1 non-terminating fraction: am(x) = 0, bm(x) = 0 a1(x) b1(x) + a2(x) b2(x) + a3(x) b3(x) + · · · = a1(x)/b1(x) 1 + a2(x)/b1(x)b2(x) 1 + a3(x)/b2(x)b3(x) 1 + · · ·

13 / 34

slide-14
SLIDE 14

Special functions Building blocks

⊲ Equivalence and

notation Interval Sequence Theorem Reliability

Implementation

Equivalence and notation

approximant fN(x; wN) =

N−1

K

m=1

am(x) 1 + aN(x) 1 + wN lim

N→∞ fN(x; wN) = f(x)

tail f(N) =

∞

K

m=N+1

am(x) 1 lim

N→∞ f(N) ?

= 0

14 / 34

slide-15
SLIDE 15

Special functions Building blocks

⊲ Equivalence and

notation Interval Sequence Theorem Reliability

Implementation

Example:

∞

K

m=1

x 1 = −1 + √1 + 4x 2 f(N) = −1 + √1 + 4x 2 1 =

∞

K

m=1

m(m + 2) 1 f(N) = N + 1 √ 2 − 1 = 1 1 + 2 1 + 1 1 + 2 1 + 1 1 + · · · f(2N) → √ 2 − 1 f(2N+1) → √ 2

15 / 34

slide-16
SLIDE 16

Special functions Building blocks

Equivalence and notation

⊲ Interval

Sequence Theorem Reliability

Implementation

Interval Sequence Theorem

sequence of element sets {En}n∈N sequence of value sets {Vn}n∈N an ∈ En ⇒ an 1 + Vn ⊆ Vn−1, n 1 f(N) ∈ ¯ VN, f(x) ∈ ¯ V0 wN ∈ VN ⇒ fN(x; wN) ∈ V0 Example: 2 1 + 1 1 − 1 1 + 2 1 + 1 1 − 1 1 + 2 1 + · · · E3n = {−1}, E3n+1 = {2}, E3n+2 = {1} V3n = {1/2}, V3n+1 = {3}, V3n+2 = {−2/3}

16 / 34

slide-17
SLIDE 17

Special functions Building blocks

Equivalence and notation

⊲ Interval

Sequence Theorem Reliability

Implementation

◮ element sets

am(x) ∈ Em = [α(l)

m , α(u) m ]

α(l)

m α(u) m 0 ◮ value sets

Vk = [Lk, Rk] Lk = min

am∈Em mk+1

f(k), Rk = max

am∈Em mk+1

f(k)

17 / 34

slide-18
SLIDE 18

Special functions Building blocks

Equivalence and notation

⊲ Interval

Sequence Theorem Reliability

Implementation

Mk := max {|Lk/(1 + Lk)|, |Rk/(1 + Rk)|} wN ∈ VN = [LN, RN]

  • f(x) − fN(x; wN)

f(x)

  • RN − LN

1 + LN

N−1

  • k=1

Mk

18 / 34

slide-19
SLIDE 19

Special functions Building blocks

Equivalence and notation Interval Sequence Theorem

⊲ Reliability

Implementation

Reliability

Relative truncation error =

  • f(x) − fN(x; wN)

f(x)

  • ǫT

Relative rounding error =

  • fN(x; wN) − FN(x; wN)

f(x)

  • ǫR

ǫ := ǫT + ǫR f(x) ∈ [FN(x; wN)(1 − 2ǫ), FN(x; wN)(1 + 2ǫ)]

19 / 34

slide-20
SLIDE 20

Special functions Building blocks

Equivalence and notation Interval Sequence Theorem

⊲ Reliability

Implementation

Example: f(a, x) = aγ(a, x)ex xa =

a a−x

1 +

∞

K

m=2 (m−1)x (a+m−1−x)(a+m−2−x)

1 Evaluate f(9/2, 1) with ǫT + ǫR 10−d+1, d = 73, 74, . . . , 80

20 / 34

slide-21
SLIDE 21

Special functions Building blocks

Equivalence and notation Interval Sequence Theorem

⊲ Reliability

Implementation

◮ Continued fraction library output

73 1.214009591773512617777498734645198390079596056622283491877162409691879700 74 1.2140095917735126177774987346451983900795960566222834918771624096918797000 75 1.21400959177351261777749873464519839007959605662228349187716240969187969998 76 1.214009591773512617777498734645198390079596056622283491877162409691879699983 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

slide-22
SLIDE 22

Special functions Building blocks ⊲ Implementation

Book Demo

Implementation

◮ Book with tables and graphs ◮ C++ library (all functions) ◮ Maple library (all fractions) ◮ Web version to explore both

22 / 34

slide-23
SLIDE 23

Special functions Building blocks Implementation

⊲ Book

Demo

Formulas

23 / 34

slide-24
SLIDE 24

Special functions Building blocks Implementation

⊲ Book

Demo

Tables

24 / 34

slide-25
SLIDE 25

Special functions Building blocks Implementation

⊲ Book

Demo

25 / 34

slide-26
SLIDE 26

Special functions Building blocks Implementation

⊲ Book

Demo

Figures

26 / 34

slide-27
SLIDE 27

Special functions Building blocks Implementation

⊲ Book

Demo

27 / 34

slide-28
SLIDE 28

Special functions Building blocks Implementation

Book

⊲ Demo

Demo

28 / 34

slide-29
SLIDE 29

Special functions Building blocks Implementation

Book

⊲ Demo

29 / 34

slide-30
SLIDE 30

Special functions Building blocks Implementation

Book

⊲ Demo

30 / 34

slide-31
SLIDE 31

Special functions Building blocks Implementation

Book

⊲ Demo

31 / 34

slide-32
SLIDE 32

Special functions Building blocks Implementation

Book

⊲ Demo

32 / 34

slide-33
SLIDE 33

Special functions Building blocks Implementation

Book

⊲ Demo

33 / 34

slide-34
SLIDE 34

Special functions Building blocks Implementation

Book

⊲ Demo

34 / 34