[PPT] - the Correctness of math.h Implementations Wonyeol Lee 1,2 Rahul PowerPoint Presentation

SLIDE 1

/146

On Automatically Proving the Correctness of math.h Implementations

Wonyeol Lee1,2 Rahul Sharma3 Alex Aiken1

1 Stanford University 2 KAIST 3 Microsoft Research

1

POPL 2018

SLIDE 2

/146

Industry standard implementations of math.h claim: “precision loss is less than 1 ulp”

Our Goal

𝑚𝑝𝑕 𝑦

mathematical specification 2

SLIDE 3

/146

Industry standard implementations of math.h claim: “precision loss is less than 1 ulp”

Our Goal

math.h implementation

<log>

... mulpd %xmm2, %xmm6 addsd %xmm0, %xmm5 addsd %xmm7, %xmm3 addsd %xmm5, %xmm1 ...

𝑚𝑝𝑕 𝑦

mathematical specification 3

SLIDE 4

/146

Industry standard implementations of math.h claim: “precision loss is less than 1 ulp”

Our Goal

math.h implementation

<log>

... mulpd %xmm2, %xmm6 addsd %xmm0, %xmm5 addsd %xmm7, %xmm3 addsd %xmm5, %xmm1 ...

𝑚𝑝𝑕 𝑦

mathematical specification 4

SLIDE 5

/146

Industry standard implementations of math.h claim: “precision loss is less than 1 ulp”

Our Goal

math.h implementation

<log>

... mulpd %xmm2, %xmm6 addsd %xmm0, %xmm5 addsd %xmm7, %xmm3 addsd %xmm5, %xmm1 ...

𝑚𝑝𝑕 𝑦

mathematical specification 5 infinite bits fixed bits

SLIDE 6

/146

Industry standard implementations of math.h claim: “precision loss is less than 1 ulp”

Our Goal

≠

math.h implementation

<log>

... mulpd %xmm2, %xmm6 addsd %xmm0, %xmm5 addsd %xmm7, %xmm3 addsd %xmm5, %xmm1 ...

𝑚𝑝𝑕 𝑦

mathematical specification 6 infinite bits fixed bits

SLIDE 7

/146

Industry standard implementations of math.h claim: “precision loss is less than 1 ulp”

Our Goal

≠

math.h implementation

<log>

... mulpd %xmm2, %xmm6 addsd %xmm0, %xmm5 addsd %xmm7, %xmm3 addsd %xmm5, %xmm1 ...

𝑚𝑝𝑕 𝑦

mathematical specification 7

∃ precision loss

infinite bits fixed bits

SLIDE 8

/146

Industry standard implementations of math.h claim: “precision loss is less than 1 ulp”

Our Goal

≠

math.h implementation

<log>

... mulpd %xmm2, %xmm6 addsd %xmm0, %xmm5 addsd %xmm7, %xmm3 addsd %xmm5, %xmm1 ...

𝑚𝑝𝑕 𝑦

mathematical specification 8

∃ precision loss

infinite bits fixed bits

SLIDE 9

/146

Industry standard implementations of math.h claim: “precision loss is less than 1 ulp”

Our Goal

9

… −1 1 …

SLIDE 10

/146

Industry standard implementations of math.h claim: “precision loss is less than 1 ulp”

Our Goal

10

𝑚𝑝𝑕 𝑦 … −1 1 …

SLIDE 11

/146

Industry standard implementations of math.h claim: “precision loss is less than 1 ulp”

Our Goal

11

𝑚𝑝𝑕 𝑦 … −1 1 … log has precision loss of < 1 ulp:

SLIDE 12

/146

Industry standard implementations of math.h claim: “precision loss is less than 1 ulp”

Our Goal

12

𝑚𝑝𝑕 𝑦 log(𝑦) … −1 1 … log has precision loss of < 1 ulp:

SLIDE 13

/146

Industry standard implementations of math.h claim: “precision loss is less than 1 ulp”

Our Goal

13

𝑚𝑝𝑕 𝑦 log(𝑦) … −1 1 … log has precision loss of < 1 ulp:

SLIDE 14

/146

Industry standard implementations of math.h claim: “precision loss is less than 1 ulp”

Our Goal

14

𝑚𝑝𝑕 𝑦 log(𝑦) … −1 1 … log has precision loss of < 1 ulp:

Goal: Prove this claim automatically!

SLIDE 15

/146

Example:
Rounding errors
Floating-point arithmetic ≠ Real arithmetic

1 ⊕ 2100 ⊖ 2100 = 1 ≠ 0 = 1 ⊕ 2100 ⊖ 2100

Floating-point implementations often have precision loss

Floating-Point Numbers/Operations

15

1 01111111111 1100⋯00 (2) = −1 1 ∙ 21023−1023 ∙ 1.110 ⋯ 00 2

significand

SLIDE 16

/146

Example:
Rounding errors
Floating-point arithmetic ≠ Real arithmetic

1 ⊕ 2100 ⊖ 2100 = 1 ≠ 0 = 1 ⊕ 2100 ⊖ 2100

Floating-point implementations often have precision loss

Floating-Point Numbers/Operations

16

1 01111111111 1100⋯00 (2) = −1 1 ∙ 21023−1023 ∙ 1.110 ⋯ 00 2

significand

SLIDE 17

/146

Example:
Rounding errors
Floating-point arithmetic ≠ Real arithmetic

1 ⊕ 2100 ⊖ 2100 = 1 ≠ 0 = 1 ⊕ 2100 ⊖ 2100

Floating-point implementations often have precision loss

Floating-Point Numbers/Operations

17

1 01111111111 1100⋯00 (2) = −1 1 ∙ 21023−1023 ∙ 1.110 ⋯ 00 2

significand

floating-point

peration

SLIDE 18

/146

Example:
Rounding errors
Floating-point arithmetic ≠ Real arithmetic

1 ⊕ 2100 ⊖ 2100 = 1 ≠ 0 = 1 ⊕ 2100 ⊖ 2100

Floating-point implementations often have precision loss

Floating-Point Numbers/Operations

18

1 01111111111 1100⋯00 (2) = −1 1 ∙ 21023−1023 ∙ 1.110 ⋯ 00 2

significand

floating-point

peration

SLIDE 19

/146

Ulp Error

Typically used to measure accuracy of numeric libraries
Definition:
𝑏 = exact value, 𝑐 = approximate value

ErrUlp 𝑏, 𝑐 = 𝑏 − 𝑐 ulp 𝑏

19

SLIDE 20

/146

Ulp Error

Typically used to measure accuracy of numeric libraries
Definition:
𝑏 = exact value, 𝑐 = approximate value

20

𝑏 𝑐

SLIDE 21

/146

Ulp Error

Typically used to measure accuracy of numeric libraries
Definition:
𝑏 = exact value, 𝑐 = approximate value
ulp 𝑏 = gap between two doubles that surround 𝑏

21

𝑏 𝑐 ulp 𝑏

SLIDE 22

/146

Ulp Error

Typically used to measure accuracy of numeric libraries
Definition:
𝑏 = exact value, 𝑐 = approximate value
ulp 𝑏 = gap between two doubles that surround 𝑏

ErrUlp 𝑏, 𝑐 = 𝑏 − 𝑐 ulp 𝑏

22

𝑏 𝑐 ulp 𝑏

SLIDE 23

/146

Ulp Error

Typically used to measure accuracy of numeric libraries
Definition:
𝑏 = exact value, 𝑐 = approximate value
ulp 𝑏 = gap between two doubles that surround 𝑏

ErrUlp 𝑏, 𝑐 = 𝑏 − 𝑐 ulp 𝑏

23

𝑏 𝑐 ulp 𝑏 2.7 ulps

SLIDE 24

/146

Problem Statement

Problem: Find a small Θ > 0 automatically such that

ErrUlp(𝑔 𝑦 , 𝑄 𝑦 ) ≤ Θ for all 𝑦 ∈ 𝑌

i.e., prove a bound on the maximum precision loss

24 math.h implementation 𝑄

<log> ... mulpd %xmm2, %xmm6 addsd %xmm0, %xmm5 addsd %xmm7, %xmm3 addsd %xmm5, %xmm1 ...

𝑚𝑝𝑕 𝑦

mathematical specification 𝑔

SLIDE 25

/146

Problem Statement

Problem: Find a small Θ > 0 automatically such that

ErrUlp(𝑔 𝑦 , 𝑄 𝑦 ) ≤ Θ for all 𝑦 ∈ 𝑌

i.e., prove a bound on the maximum precision loss

25 input interval 𝑌

(0, 21024)

math.h implementation 𝑄

<log> ... mulpd %xmm2, %xmm6 addsd %xmm0, %xmm5 addsd %xmm7, %xmm3 addsd %xmm5, %xmm1 ...

𝑚𝑝𝑕 𝑦

mathematical specification 𝑔

SLIDE 26

/146

Problem Statement

Problem: Find a small Θ > 0 automatically such that

ErrUlp(𝑔 𝑦 , 𝑄 𝑦 ) ≤ Θ for all 𝑦 ∈ 𝑌

i.e., prove a bound on the maximum precision loss

26 input interval 𝑌

(0, 21024)

math.h implementation 𝑄

<log> ... mulpd %xmm2, %xmm6 addsd %xmm0, %xmm5 addsd %xmm7, %xmm3 addsd %xmm5, %xmm1 ...

𝑚𝑝𝑕 𝑦

mathematical specification 𝑔

SLIDE 27

/146

Problem Statement

Problem: Find a small Θ > 0 automatically such that

ErrUlp(𝑔 𝑦 , 𝑄 𝑦 ) ≤ Θ for all 𝑦 ∈ 𝑌

i.e., prove a bound on the maximum precision loss

27 input interval 𝑌

(0, 21024)

math.h implementation 𝑄

<log> ... mulpd %xmm2, %xmm6 addsd %xmm0, %xmm5 addsd %xmm7, %xmm3 addsd %xmm5, %xmm1 ...

𝑚𝑝𝑕 𝑦

mathematical specification 𝑔

SLIDE 28

/146

Problem Statement

Problem: Find a small Θ > 0 automatically such that

ErrUlp(𝑔 𝑦 , 𝑄 𝑦 ) ≤ Θ for all 𝑦 ∈ 𝑌

i.e., prove a bound on the maximum precision loss

28 input interval 𝑌

(0, 21024)

math.h implementation 𝑄

<log> ... mulpd %xmm2, %xmm6 addsd %xmm0, %xmm5 addsd %xmm7, %xmm3 addsd %xmm5, %xmm1 ...

𝑚𝑝𝑕 𝑦

mathematical specification 𝑔

Goal: want to find Θ < 1

SLIDE 29

/146

Previous Work

Machine-checkable proofs
Harrison used HOL Light to prove [FMCAD’00]:

“Intel’s sin has precision loss of < 0.574 ulps”

“Construction of these proofs often requires

considerable persistence.” [FMSD’00]

Automatic techniques
Astree [PLDI’03], Fluctuat [FMICS’09], Gappa [TOMS’10],

MathSAT [FMCAD’12], Rosa [POPL’14], FPTaylor [FM’15], Lee/Sharma/Aiken [PLDI’16], ⋯

None of them can prove < 1 ulp error bound automatically

29

SLIDE 30

/146

Previous Work

Machine-checkable proofs
Harrison used HOL Light to prove [FMCAD’00]:

“Intel’s sin has precision loss of < 0.574 ulps”

“Construction of these proofs often requires

considerable persistence.” [FMSD’00]

Automatic techniques
Astree [PLDI’03], Fluctuat [FMICS’09], Gappa [TOMS’10],

MathSAT [FMCAD’12], Rosa [POPL’14], FPTaylor [FM’15], Lee/Sharma/Aiken [PLDI’16], ⋯

None of them can prove < 1 ulp error bound automatically

30

SLIDE 31

/146

Previous Work

Machine-checkable proofs
Harrison used HOL Light to prove [FMCAD’00]:

“Intel’s sin has precision loss of < 0.574 ulps”

“Construction of these proofs often requires

considerable persistence.” [FMSD’00]

Automatic techniques
Astree [PLDI’03], Fluctuat [FMICS’09], Gappa [TOMS’10],

MathSAT [FMCAD’12], Rosa [POPL’14], FPTaylor [FM’15], Lee/Sharma/Aiken [PLDI’16], ⋯

None of them can prove < 1 ulp error bound automatically

31

SLIDE 32

/146

Previous Work

Machine-checkable proofs
Harrison used HOL Light to prove [FMCAD’00]:

“Intel’s sin has precision loss of < 0.574 ulps”

“Construction of these proofs often requires

considerable persistence.” [FMSD’00]

Automatic techniques
Astree [PLDI’03], Fluctuat [FMICS’09], Gappa [TOMS’10],

MathSAT [FMCAD’12], Rosa [POPL’14], FPTaylor [FM’15], Lee/Sharma/Aiken [PLDI’16], ⋯

None of them can prove < 1 ulp error bound automatically

32

SLIDE 33

/146

1 + 𝜗 Property

A standard way to model rounding errors
Theorem For any double 𝑏 and 𝑐, and ∗ ∈ {+, −,×,/},

𝑏 ⊛ 𝑐 = 𝑏 ∗ 𝑐 1 + 𝜀 for some 𝜀 < 𝜗

where 𝜗 = 2−53

That is,

33

SLIDE 34

/146

1 + 𝜗 Property

A standard way to model rounding errors
Theorem For any double 𝑏 and 𝑐, and ∗ ∈ {+, −,×,/},

𝑏 ⊛ 𝑐 = 𝑏 ∗ 𝑐 1 + 𝜀 for some 𝜀 < 𝜗

where 𝜗 = 2−53

34

SLIDE 35

/146

1 + 𝜗 Property

A standard way to model rounding errors
Theorem For any double 𝑏 and 𝑐, and ∗ ∈ {+, −,×,/},

𝑏 ⊛ 𝑐 = 𝑏 ∗ 𝑐 1 + 𝜀 for some 𝜀 < 𝜗

where 𝜗 = 2−53

That is,

35

𝑏 ⊛ 𝑐

1

SLIDE 36

/146

1 + 𝜗 Property

A standard way to model rounding errors
Theorem For any double 𝑏 and 𝑐, and ∗ ∈ {+, −,×,/},

𝑏 ⊛ 𝑐 = 𝑏 ∗ 𝑐 1 + 𝜀 for some 𝜀 < 𝜗

where 𝜗 = 2−53

That is,

36

𝑏 ⊛ 𝑐 𝑏 ∗ 𝑐

1

SLIDE 37

/146

1 + 𝜗 Property

A standard way to model rounding errors
Theorem For any double 𝑏 and 𝑐, and ∗ ∈ {+, −,×,/},

𝑏 ⊛ 𝑐 = 𝑏 ∗ 𝑐 1 + 𝜀 for some 𝜀 < 𝜗

where 𝜗 = 2−53

That is,

37

𝑏 ⊛ 𝑐 𝑏 ∗ 𝑐

1

𝑏 ∗ 𝑐 1 + 𝜀 ∶ 𝜀 < 𝜗

SLIDE 38

/146

1 + 𝜗 Property

A standard way to model rounding errors
Theorem For any double 𝑏 and 𝑐, and ∗ ∈ {+, −,×,/},

𝑏 ⊛ 𝑐 = 𝑏 ∗ 𝑐 1 + 𝜀 for some 𝜀 < 𝜗

where 𝜗 = 2−53

That is,

38

𝑏 ⊛ 𝑐 𝑏 ∗ 𝑐

1

𝑏 ∗ 𝑐 1 + 𝜀 ∶ 𝜀 < 𝜗

SLIDE 39

/146

Example:

𝑔 𝑦 = 1 + 𝑦 over 𝑌 = [−0.1, 0.1] 𝑓 = 1 ⊕ 0.5 ⊗ 𝑦 ⟵ implementation of 𝑔 𝑦

1. Construct abstraction 𝐵𝜀 of 𝑓:

𝐵𝜀 𝑦 = 1 + 0.5 × 𝑦 1 + 𝜀1 1 + 𝜀2

𝐵𝜀 is a sound abstraction of 𝑓 (or 𝑓 ⊑ 𝐵𝜀):

𝑓 𝑦 ∈ {𝐵𝜀 𝑦 ∶ 𝜀1 , 𝜀2 < 𝜗} for all 𝑦 ∈ 𝑌

Standard Error Analysis

39

SLIDE 40

/146

Example:

𝑔 𝑦 = 1 + 𝑦 over 𝑌 = [−0.1, 0.1] 𝑓 = 1 ⊕ 0.5 ⊗ 𝑦 ⟵ implementation of 𝑔 𝑦

1. Construct abstraction 𝐵𝜀 of 𝑓:

𝐵𝜀 𝑦 = 1 + 0.5 × 𝑦 1 + 𝜀1 1 + 𝜀2

𝐵𝜀 is a sound abstraction of 𝑓 (or 𝑓 ⊑ 𝐵𝜀):

𝑓 𝑦 ∈ {𝐵𝜀 𝑦 ∶ 𝜀1 , 𝜀2 < 𝜗} for all 𝑦 ∈ 𝑌

Standard Error Analysis

40

SLIDE 41

/146

Example:

𝑔 𝑦 = 1 + 𝑦 over 𝑌 = [−0.1, 0.1] 𝑓 = 1 ⊕ 0.5 ⊗ 𝑦 ⟵ implementation of 𝑔 𝑦

1. Construct abstraction 𝐵𝜀 of 𝑓:

𝐵𝜀 𝑦 = 1 + 0.5 × 𝑦 1 + 𝜀1 1 + 𝜀2

𝐵𝜀 is a sound abstraction of 𝑓 (or 𝑓 ⊑ 𝐵𝜀):

𝑓 𝑦 ∈ {𝐵𝜀 𝑦 ∶ 𝜀1 , 𝜀2 < 𝜗} for all 𝑦 ∈ 𝑌

Standard Error Analysis

41

SLIDE 42

/146

Example:

𝑔 𝑦 = 1 + 𝑦 over 𝑌 = [−0.1, 0.1] 𝑓 = 1 ⊕ 0.5 ⊗ 𝑦 ⟵ implementation of 𝑔 𝑦

1. Construct abstraction 𝐵𝜀 of 𝑓:

𝐵𝜀 𝑦 = 1 + 0.5 × 𝑦 1 + 𝜀1 1 + 𝜀2

𝐵𝜀 is a sound abstraction of 𝑓 (or 𝑓 ⊑ 𝐵𝜀):

𝑓 𝑦 ∈ {𝐵𝜀 𝑦 ∶ 𝜀1 , 𝜀2 < 𝜗} for all 𝑦 ∈ 𝑌

Standard Error Analysis

42

⊕ ⊗ =

SLIDE 43

/146

Example:

𝑔 𝑦 = 1 + 𝑦 over 𝑌 = [−0.1, 0.1] 𝑓 = 1 ⊕ 0.5 ⊗ 𝑦 ⟵ implementation of 𝑔 𝑦

1. Construct abstraction 𝐵𝜀 of 𝑓:

𝐵𝜀 𝑦 = 1 + 0.5 × 𝑦 1 + 𝜀1 1 + 𝜀2

𝐵𝜀 is a sound abstraction of 𝑓 (or 𝑓 ⊑ 𝐵𝜀):

𝑓 𝑦 ∈ {𝐵𝜀 𝑦 ∶ 𝜀1 , 𝜀2 < 𝜗} for all 𝑦 ∈ 𝑌

Standard Error Analysis

43

=

SLIDE 44

/146

Example:

𝑔 𝑦 = 1 + 𝑦 over 𝑌 = [−0.1, 0.1] 𝑓 = 1 ⊕ 0.5 ⊗ 𝑦 ⟵ implementation of 𝑔 𝑦

1. Construct abstraction 𝐵𝜀 of 𝑓:

𝐵𝜀 𝑦 = 1 + 0.5 × 𝑦 1 + 𝜀1 1 + 𝜀2

𝐵𝜀 is a sound abstraction of 𝑓 (or 𝑓 ⊑ 𝐵𝜀):

𝑓 𝑦 ∈ {𝐵𝜀 𝑦 ∶ 𝜀1 , 𝜀2 < 𝜗} for all 𝑦 ∈ 𝑌

Standard Error Analysis

44

SLIDE 45

/146

Example:

𝑔 𝑦 = 1 + 𝑦 over 𝑌 = [−0.1, 0.1] 𝑓 = 1 ⊕ 0.5 ⊗ 𝑦 ⟵ implementation of 𝑔 𝑦

1. Construct abstraction 𝐵𝜀 of 𝑓:

𝐵𝜀 𝑦 = 1 + 0.5 × 𝑦 1 + 𝜀1 1 + 𝜀2

𝐵𝜀 is a sound abstraction of 𝑓 (or 𝑓 ⊑ 𝐵𝜀):

𝑓 𝑦 ∈ {𝐵𝜀 𝑦 ∶ 𝜀1 , 𝜀2 < 𝜗} for all 𝑦 ∈ 𝑌

Standard Error Analysis

45

SLIDE 46

/146

Example:

𝑔 𝑦 = 1 + 𝑦 over 𝑌 = [−0.1, 0.1] 𝑓 = 1 ⊕ 0.5 ⊗ 𝑦 ⟵ implementation of 𝑔 𝑦

2. Compute a bound on relative error of 𝑓:

Θ𝑠𝑓𝑚 = max

𝑦∈ −0.1,0.1 , 𝜀1,𝜀2∈ −𝜗,𝜗 𝑔(𝑦)−𝐵𝜀 𝑦 𝑔(𝑦)

= max

𝑦∈ −0.1,0.1 , 𝜀1,𝜀2∈ −𝜗,𝜗 1+𝑦− 1+ 0.5 × 𝑦 1+𝜀1 1+𝜀2 1+𝑦

3. Compute a bound on ulp error of 𝑓:

Θ𝑣𝑚𝑞 = Θ𝑠𝑓𝑚/𝜗

Standard Error Analysis

46

SLIDE 47

/146

Example:

𝑔 𝑦 = 1 + 𝑦 over 𝑌 = [−0.1, 0.1] 𝑓 = 1 ⊕ 0.5 ⊗ 𝑦 ⟵ implementation of 𝑔 𝑦

2. Compute a bound on relative error of 𝑓:

Θ𝑠𝑓𝑚 = max

𝑦∈ −0.1,0.1 , 𝜀1,𝜀2∈ −𝜗,𝜗 𝑔(𝑦)−𝐵𝜀 𝑦 𝑔(𝑦)

= max

𝑦∈ −0.1,0.1 , 𝜀1,𝜀2∈ −𝜗,𝜗 1+𝑦− 1+ 0.5 × 𝑦 1+𝜀1 1+𝜀2 1+𝑦

3. Compute a bound on ulp error of 𝑓:

Θ𝑣𝑚𝑞 = Θ𝑠𝑓𝑚/𝜗

Standard Error Analysis

47

SLIDE 48

/146

Example:

𝑔 𝑦 = 1 + 𝑦 over 𝑌 = [−0.1, 0.1] 𝑓 = 1 ⊕ 0.5 ⊗ 𝑦 ⟵ implementation of 𝑔 𝑦

2. Compute a bound on relative error of 𝑓:

Θ𝑠𝑓𝑚 = max

𝑦∈ −0.1,0.1 , 𝜀1,𝜀2∈ −𝜗,𝜗 𝑔(𝑦)−𝐵𝜀 𝑦 𝑔(𝑦)

= max

𝑦∈ −0.1,0.1 , 𝜀1,𝜀2∈ −𝜗,𝜗 1+𝑦− 1+ 0.5 × 𝑦 1+𝜀1 1+𝜀2 1+𝑦

3. Compute a bound on ulp error of 𝑓:

Θ𝑣𝑚𝑞 = Θ𝑠𝑓𝑚/𝜗

Standard Error Analysis

48

Theorem ErrUlp 𝑏, 𝑐 ≤ ErrRel 𝑏, 𝑐 /𝜗

SLIDE 49

/146

Example:

𝑔 𝑦 = 1 + 𝑦 over 𝑌 = [−0.1, 0.1] 𝑓 = 1 ⊕ 0.5 ⊗ 𝑦 ⟵ implementation of 𝑔 𝑦

2. Compute a bound on relative error of 𝑓:

Θ𝑠𝑓𝑚 = max

𝑦∈ −0.1,0.1 , 𝜀1,𝜀2∈ −𝜗,𝜗 𝑔(𝑦)−𝐵𝜀 𝑦 𝑔(𝑦)

= max

𝑦∈ −0.1,0.1 , 𝜀1,𝜀2∈ −𝜗,𝜗 1+𝑦− 1+ 0.5 × 𝑦 1+𝜀1 1+𝜀2 1+𝑦

3. Compute a bound on ulp error of 𝑓:

Θ𝑣𝑚𝑞 = Θ𝑠𝑓𝑚/𝜗

Standard Error Analysis

49

Theorem ErrUlp 𝑏, 𝑐 ≤ ErrRel 𝑏, 𝑐 /𝜗

SLIDE 50

/146

Standard Error Analysis: Limitation

50

SLIDE 51

/146

Standard Error Analysis: Limitation

51

error bounds from [PLDI’16] input value

Intel’s log

1014 ⋮

ulp error

based on standard error analysis

SLIDE 52

/146

Standard Error Analysis: Limitation

52

error bounds from [PLDI’16] input value

Intel’s log

1014 ⋮

ulp error 1 ulp

based on standard error analysis

3 ulps

SLIDE 53

/146

Standard Error Analysis: Limitation

53

error bounds from [PLDI’16] input value

Intel’s log

1014 ⋮

ulp error

4095 4096 , 1

1 ulp

based on standard error analysis

3 ulps

SLIDE 54

/146

Why the Loose Error Bound?

Floating-point operations are sometimes exact

0.125 ⊗ 𝑦 = 0.125 × 𝑦 if 0.125𝑦 ≥ 2−1022 𝑦 ⊖ 1 0 = 𝑦 − 1 if 0.5 ≤ 𝑦 ≤ 2

Standard error analysis constructs imprecise abstractions

0.125 ⊗ 𝑦 ⊑ 0.125 × 𝑦 if 0.125𝑦 ≥ 2−1022 0.125 ⊗ 𝑦 ⊑ 0.125 × 𝑦 (1 + 𝜀) 𝑦 ⊖ 1 0 ⊑ 𝑦 − 1 if 0.5 ≤ 𝑦 ≤ 2 𝑦 ⊖ 1 0 ⊑ 𝑦 − 1 (1 + 𝜀′)

54

SLIDE 55

/146

Why the Loose Error Bound?

Floating-point operations are sometimes exact

0.125 ⊗ 𝑦 = 0.125 × 𝑦 if 0.125𝑦 ≥ 2−1022 𝑦 ⊖ 1 0 = 𝑦 − 1 if 0.5 ≤ 𝑦 ≤ 2

Standard error analysis constructs imprecise abstractions

0.125 ⊗ 𝑦 ⊑ 0.125 × 𝑦 if 0.125𝑦 ≥ 2−1022 0.125 ⊗ 𝑦 ⊑ 0.125 × 𝑦 (1 + 𝜀) 𝑦 ⊖ 1 0 ⊑ 𝑦 − 1 if 0.5 ≤ 𝑦 ≤ 2 𝑦 ⊖ 1 0 ⊑ 𝑦 − 1 (1 + 𝜀′)

55

SLIDE 56

/146

Why the Loose Error Bound?

Floating-point operations are sometimes exact

0.125 ⊗ 𝑦 = 0.125 × 𝑦 if 0.125𝑦 ≥ 2−1022 𝑦 ⊖ 1 0 = 𝑦 − 1 if 0.5 ≤ 𝑦 ≤ 2

Standard error analysis constructs imprecise abstractions

0.125 ⊗ 𝑦 ⊑ 0.125 × 𝑦 if 0.125𝑦 ≥ 2−1022 0.125 ⊗ 𝑦 ⊑ 0.125 × 𝑦 (1 + 𝜀) 𝑦 ⊖ 1 0 ⊑ 𝑦 − 1 if 0.5 ≤ 𝑦 ≤ 2 𝑦 ⊖ 1 0 ⊑ 𝑦 − 1 (1 + 𝜀′)

56

SLIDE 57

/146

Why the Loose Error Bound?

Floating-point operations are sometimes exact

0.125 ⊗ 𝑦 = 0.125 × 𝑦 if 0.125𝑦 ≥ 2−1022 𝑦 ⊖ 1 0 = 𝑦 − 1 if 0.5 ≤ 𝑦 ≤ 2

Standard error analysis constructs imprecise abstractions

0.125 ⊗ 𝑦 ⊑ 0.125 × 𝑦 if 0.125𝑦 ≥ 2−1022 0.125 ⊗ 𝑦 ⊑ 0.125 × 𝑦 (1 + 𝜀) 𝑦 ⊖ 1 0 ⊑ 𝑦 − 1 if 0.5 ≤ 𝑦 ≤ 2 𝑦 ⊖ 1 0 ⊑ 𝑦 − 1 (1 + 𝜀′)

57

SLIDE 58

/146

Why the Loose Error Bound?

Floating-point operations are sometimes exact

0.125 ⊗ 𝑦 = 0.125 × 𝑦 if 0.125𝑦 ≥ 2−1022 𝑦 ⊖ 1 0 = 𝑦 − 1 if 0.5 ≤ 𝑦 ≤ 2

Standard error analysis constructs imprecise abstractions

0.125 ⊗ 𝑦 ⊑ 0.125 × 𝑦 if 0.125𝑦 ≥ 2−1022 0.125 ⊗ 𝑦 ⊑ 0.125 × 𝑦 (1 + 𝜀) 𝑦 ⊖ 1 0 ⊑ 𝑦 − 1 if 0.5 ≤ 𝑦 ≤ 2 𝑦 ⊖ 1 0 ⊑ 𝑦 − 1 (1 + 𝜀′)

58

SLIDE 59

/146

Why the Loose Error Bound?

Floating-point operations are sometimes exact

0.125 ⊗ 𝑦 = 0.125 × 𝑦 if 0.125𝑦 ≥ 2−1022 𝑦 ⊖ 1 0 = 𝑦 − 1 if 0.5 ≤ 𝑦 ≤ 2

Standard error analysis constructs imprecise abstractions

0.125 ⊗ 𝑦 ⊑ 0.125 × 𝑦 if 0.125𝑦 ≥ 2−1022 0.125 ⊗ 𝑦 ⊑ 0.125 × 𝑦 (1 + 𝜀) 𝑦 ⊖ 1 0 ⊑ 𝑦 − 1 if 0.5 ≤ 𝑦 ≤ 2 𝑦 ⊖ 1 0 ⊑ 𝑦 − 1 (1 + 𝜀′)

59

SLIDE 60

/146

Why the Loose Error Bound?

Floating-point operations are sometimes exact

0.125 ⊗ 𝑦 = 0.125 × 𝑦 if 0.125𝑦 ≥ 2−1022 𝑦 ⊖ 1 0 = 𝑦 − 1 if 0.5 ≤ 𝑦 ≤ 2

Standard error analysis constructs imprecise abstractions

0.125 ⊗ 𝑦 ⊑ 0.125 × 𝑦 if 0.125𝑦 ≥ 2−1022 0.125 ⊗ 𝑦 ⊑ 0.125 × 𝑦 (1 + 𝜀) 𝑦 ⊖ 1 0 ⊑ 𝑦 − 1 if 0.5 ≤ 𝑦 ≤ 2 𝑦 ⊖ 1 0 ⊑ 𝑦 − 1 (1 + 𝜀′)

60

SLIDE 61

/146

Analysis of log

For 𝑦 ∈

4095 4096 , 1 , log computes

𝑠 𝑦 = 2 ⊗ 𝑦 ⊖ 255 128 ⊗ 1 2 ⊕ 255 128 ⊗ 1 2 ⊖ 1 ≈ 𝑦 − 1

and returns

𝑠 𝑦 − 1 2 𝑠 𝑦 2 + ⋯ + 1 7 𝑠 𝑦 7

Roughly speaking,

(precision loss of log) ≥ (precision loss of 𝑠 𝑦 )

61

SLIDE 62

/146

Analysis of log

For 𝑦 ∈

4095 4096 , 1 , log computes

𝑠 𝑦 = 2 ⊗ 𝑦 ⊖ 255 128 ⊗ 1 2 ⊕ 255 128 ⊗ 1 2 ⊖ 1 ≈ 𝑦 − 1

and returns

𝑠 𝑦 − 1 2 𝑠 𝑦 2 + ⋯ + 1 7 𝑠 𝑦 7

Roughly speaking,

(precision loss of log) ≥ (precision loss of 𝑠 𝑦 )

62

SLIDE 63

/146

Analysis of log

For 𝑦 ∈

4095 4096 , 1 , log computes

𝑠 𝑦 = 2 ⊗ 𝑦 ⊖ 255 128 ⊗ 1 2 ⊕ 255 128 ⊗ 1 2 ⊖ 1 ≈ 𝑦 − 1

and returns

𝑠 𝑦 − 1 2 𝑠 𝑦 2 + ⋯ + 1 7 𝑠 𝑦 7

Roughly speaking,

(precision loss of log) ≥ (precision loss of 𝑠 𝑦 )

63

SLIDE 64

/146

Analysis of log

For 𝑦 ∈

4095 4096 , 1 , log computes

𝑠 𝑦 = 2 ⊗ 𝑦 ⊖ 255 128 ⊗ 1 2 ⊕ 255 128 ⊗ 1 2 ⊖ 1 ≈ 𝑦 − 1

and returns

𝑠 𝑦 − 1 2 𝑠 𝑦 2 + ⋯ + 1 7 𝑠 𝑦 7

Roughly speaking,

(precision loss of log) ≥ (precision loss of 𝑠 𝑦 )

64

SLIDE 65

/146

Analysis of log

For 𝑦 ∈

4095 4096 , 1 , log computes

𝑠 𝑦 = 2 ⊗ 𝑦 ⊖ 255 128 ⊗ 1 2 ⊕ 255 128 ⊗ 1 2 ⊖ 1 ≈ 𝑦 − 1

and returns

𝑠 𝑦 − 1 2 𝑠 𝑦 2 + ⋯ + 1 7 𝑠 𝑦 7

Roughly speaking,

(precision loss of log) ≥ (precision loss of 𝑠 𝑦 )

65

SLIDE 66

/146

Standard Analysis of log

𝑠 𝑦 = 2 ⊗ 𝑦 ⊖ 255 128 ⊗ 1 2 ⊕ 255 128 ⊗ 1 2 ⊖ 1 4095 4096 ≤ 𝑦 < 1

Abstraction of 𝑠 𝑦 :

𝐵𝜀 𝑦 = 2 × 𝑦 1 + 𝜀0 − 255 128 1 + 𝜀1 × 1 2 1 + 𝜀2 + ⋯ = 𝑦 − 1 + 𝑦 − 255 256 𝜀1 + ⋯ 00000000000000

Bound on relative error of 𝑠 𝑦 :

max

. 4095 4096≤𝑦<1, 𝜀𝑗 <𝜗

𝐵𝜀 𝑦 − (𝑦 − 1) 𝑦 − 1 ≥ max

. 4095 4096≤𝑦<1

𝑦 − 255 256 𝑦 − 1 𝜗

66

SLIDE 67

/146

Standard Analysis of log

𝑠 𝑦 = 2 ⊗ 𝑦 ⊖ 255 128 ⊗ 1 2 ⊕ 255 128 ⊗ 1 2 ⊖ 1 4095 4096 ≤ 𝑦 < 1

Abstraction of 𝑠 𝑦 :

𝐵𝜀 𝑦 = 2 × 𝑦 1 + 𝜀0 − 255 128 1 + 𝜀1 × 1 2 1 + 𝜀2 + ⋯ = 𝑦 − 1 + 𝑦 − 255 256 𝜀1 + ⋯ 00000000000000

Bound on relative error of 𝑠 𝑦 :

max

. 4095 4096≤𝑦<1, 𝜀𝑗 <𝜗

𝐵𝜀 𝑦 − (𝑦 − 1) 𝑦 − 1 ≥ max

. 4095 4096≤𝑦<1

𝑦 − 255 256 𝑦 − 1 𝜗

67

SLIDE 68

/146

Standard Analysis of log

𝑠 𝑦 = 2 ⊗ 𝑦 ⊖ 255 128 ⊗ 1 2 ⊕ 255 128 ⊗ 1 2 ⊖ 1 4095 4096 ≤ 𝑦 < 1

Abstraction of 𝑠 𝑦 :

𝐵𝜀 𝑦 = 2 × 𝑦 1 + 𝜀0 − 255 128 1 + 𝜀1 × 1 2 1 + 𝜀2 + ⋯ = 𝑦 − 1 + 𝑦 − 255 256 𝜀1 + ⋯ 00000000000000

Bound on relative error of 𝑠 𝑦 :

max

. 4095 4096≤𝑦<1, 𝜀𝑗 <𝜗

𝐵𝜀 𝑦 − (𝑦 − 1) 𝑦 − 1 ≥ max

. 4095 4096≤𝑦<1

𝑦 − 255 256 𝑦 − 1 𝜗

68

SLIDE 69

/146

Standard Analysis of log

𝑠 𝑦 = 2 ⊗ 𝑦 ⊖ 255 128 ⊗ 1 2 ⊕ 255 128 ⊗ 1 2 ⊖ 1 4095 4096 ≤ 𝑦 < 1

Abstraction of 𝑠 𝑦 :

𝐵𝜀 𝑦 = 2 × 𝑦 1 + 𝜀0 − 255 128 1 + 𝜀1 × 1 2 1 + 𝜀2 + ⋯ = 𝑦 − 1 + 𝑦 − 255 256 𝜀1 + ⋯ 00000000000000

Bound on relative error of 𝑠 𝑦 :

max

. 4095 4096≤𝑦<1, 𝜀𝑗 <𝜗

𝐵𝜀 𝑦 − (𝑦 − 1) 𝑦 − 1 ≥ max

. 4095 4096≤𝑦<1

𝑦 − 255 256 𝑦 − 1 𝜗

69

SLIDE 70

/146

Standard Analysis of log

𝑠 𝑦 = 2 ⊗ 𝑦 ⊖ 255 128 ⊗ 1 2 ⊕ 255 128 ⊗ 1 2 ⊖ 1 4095 4096 ≤ 𝑦 < 1

Abstraction of 𝑠 𝑦 :

𝐵𝜀 𝑦 = 2 × 𝑦 1 + 𝜀0 − 255 128 1 + 𝜀1 × 1 2 1 + 𝜀2 + ⋯ = 𝑦 − 1 + 𝑦 − 255 256 𝜀1 + ⋯ 00000000000000

Bound on relative error of 𝑠 𝑦 :

max

. 4095 4096≤𝑦<1, 𝜀𝑗 <𝜗

𝐵𝜀 𝑦 − (𝑦 − 1) 𝑦 − 1 ≥ max

. 4095 4096≤𝑦<1

𝑦 − 255 256 𝑦 − 1 𝜗

70

SLIDE 71

/146

Standard Analysis of log

𝑠 𝑦 = 2 ⊗ 𝑦 ⊖ 255 128 ⊗ 1 2 ⊕ 255 128 ⊗ 1 2 ⊖ 1 4095 4096 ≤ 𝑦 < 1

Abstraction of 𝑠 𝑦 :

𝐵𝜀 𝑦 = 2 × 𝑦 1 + 𝜀0 − 255 128 1 + 𝜀1 × 1 2 1 + 𝜀2 + ⋯ = 𝑦 − 1 + 𝑦 − 255 256 𝜀1 + ⋯ 00000000000000

Bound on relative error of 𝑠 𝑦 :

max

. 4095 4096≤𝑦<1, 𝜀𝑗 <𝜗

𝐵𝜀 𝑦 − (𝑦 − 1) 𝑦 − 1 ≥ max

. 4095 4096≤𝑦<1

𝑦 − 255 256 𝑦 − 1 𝜗

71

SLIDE 72

/146

Standard Analysis of log

𝑠 𝑦 = 2 ⊗ 𝑦 ⊖ 255 128 ⊗ 1 2 ⊕ 255 128 ⊗ 1 2 ⊖ 1 4095 4096 ≤ 𝑦 < 1

Abstraction of 𝑠 𝑦 :

𝐵𝜀 𝑦 = 2 × 𝑦 1 + 𝜀0 − 255 128 1 + 𝜀1 × 1 2 1 + 𝜀2 + ⋯ = 𝑦 − 1 + 𝑦 − 255 256 𝜀1 + ⋯ 00000000000000

Bound on relative error of 𝑠 𝑦 :

max

. 4095 4096≤𝑦<1, 𝜀𝑗 <𝜗

𝐵𝜀 𝑦 − (𝑦 − 1) 𝑦 − 1 ≥ max

. 4095 4096≤𝑦<1

𝑦 − 255 256 𝑦 − 1 𝜗

72

unbounded

SLIDE 73

/146

Standard Analysis of log

𝑠 𝑦 = 2 ⊗ 𝑦 ⊖ 255 128 ⊗ 1 2 ⊕ 255 128 ⊗ 1 2 ⊖ 1 4095 4096 ≤ 𝑦 < 1

Abstraction of 𝑠 𝑦 :

𝐵𝜀 𝑦 = 2 × 𝑦 1 + 𝜀0 − 255 128 1 + 𝜀1 × 1 2 1 + 𝜀2 + ⋯ = 𝑦 − 1 + 𝑦 − 255 256 𝜀1 + ⋯ 00000000000000

Bound on relative error of 𝑠 𝑦 :

max

. 4095 4096≤𝑦<1, 𝜀𝑗 <𝜗

𝐵𝜀 𝑦 − (𝑦 − 1) 𝑦 − 1 ≥ max

. 4095 4096≤𝑦<1

𝑦 − 255 256 𝑦 − 1 𝜗

73

1014 ulps for log

near 𝑦 = 1 [PLDI’16]

unbounded

SLIDE 74

/146

𝑠 𝑦 = 2 ⊗ 𝑦 ⊖ 255 128 ⊗ 1 2 ⊕ 255 128 ⊗ 1 2 ⊖ 1 4095 4096 ≤ 𝑦 < 1 max

. 4095 4096≤𝑦<1, 𝜀′ <𝜗

𝐵′𝜀 𝑦 − (𝑦 − 1) 𝑦 − 1 = max

. 4095 4096≤𝑦<1

𝑦 − 1 𝑦 − 1 𝜗 0000000000000 = 𝜗

Precise Analysis of log

74

SLIDE 75

/146

𝑠 𝑦 = 2 ⊗ 𝑦 ⊖ 255 128 ⊗ 1 2 ⊕ 255 128 ⊗ 1 2 ⊖ 1 4095 4096 ≤ 𝑦 < 1

Precise Analysis of log

75

exact

SLIDE 76

/146

𝑠 𝑦 = 2 ⊗ 𝑦 ⊖ 255 128 ⊗ 1 2 ⊕ 255 128 ⊗ 1 2 ⊖ 1 4095 4096 ≤ 𝑦 < 1

Precise Analysis of log

76

exact

SLIDE 77

/146

𝑠 𝑦 = 2 ⊗ 𝑦 ⊖ 255 128 ⊗ 1 2 ⊕ 255 128 ⊗ 1 2 ⊖ 1 4095 4096 ≤ 𝑦 < 1

Precise Analysis of log

77

exact

SLIDE 78

/146

𝑠 𝑦 = 2 ⊗ 𝑦 ⊖ 255 128 ⊗ 1 2 ⊕ 255 128 ⊗ 1 2 ⊖ 1 4095 4096 ≤ 𝑦 < 1

Precise Analysis of log

78

exact

SLIDE 79

/146

𝑠 𝑦 = 2 ⊗ 𝑦 ⊖ 255 128 ⊗ 1 2 ⊕ 255 128 ⊗ 1 2 ⊖ 1 4095 4096 ≤ 𝑦 < 1

More precise abstraction of 𝑠 𝑦 :

𝐵′𝜀 𝑦 = 𝑦 − 1 + 𝑦 − 1 𝜀′

Precise Analysis of log

79

exact

SLIDE 80

/146

𝑠 𝑦 = 2 ⊗ 𝑦 ⊖ 255 128 ⊗ 1 2 ⊕ 255 128 ⊗ 1 2 ⊖ 1 4095 4096 ≤ 𝑦 < 1

More precise abstraction of 𝑠 𝑦 :

𝐵′𝜀 𝑦 = 𝑦 − 1 + 𝑦 − 1 𝜀′

Tighter bound on relative error of 𝑠 𝑦 :

max

. 4095 4096≤𝑦<1, 𝜀′ <𝜗

𝐵′𝜀 𝑦 − (𝑦 − 1) 𝑦 − 1 = max

. 4095 4096≤𝑦<1

𝑦 − 1 𝑦 − 1 𝜗 0000000000000 = 𝜗

Precise Analysis of log

80

exact

SLIDE 81

/146

𝑠 𝑦 = 2 ⊗ 𝑦 ⊖ 255 128 ⊗ 1 2 ⊕ 255 128 ⊗ 1 2 ⊖ 1 4095 4096 ≤ 𝑦 < 1

More precise abstraction of 𝑠 𝑦 :

𝐵′𝜀 𝑦 = 𝑦 − 1 + 𝑦 − 1 𝜀′

Tighter bound on relative error of 𝑠 𝑦 :

max

. 4095 4096≤𝑦<1, 𝜀′ <𝜗

𝐵′𝜀 𝑦 − (𝑦 − 1) 𝑦 − 1 = max

. 4095 4096≤𝑦<1

𝑦 − 1 𝑦 − 1 𝜗 0000000000000 = 𝜗

Precise Analysis of log

81

exact

We prove a bound of

0.583 ulps for log

SLIDE 82

/146

𝑠 𝑦 = 2 ⊗ 𝑦 ⊖ 255 128 ⊗ 1 2 ⊕ 255 128 ⊗ 1 2 ⊖ 1 4095 4096 ≤ 𝑦 < 1

More precise abstraction of 𝑠 𝑦 :

𝐵′𝜀 𝑦 = 𝑦 − 1 + 𝑦 − 1 𝜀′

Tighter bound on relative error of 𝑠 𝑦 :

max

. 4095 4096≤𝑦<1, 𝜀′ <𝜗

𝐵′𝜀 𝑦 − (𝑦 − 1) 𝑦 − 1 = max

. 4095 4096≤𝑦<1

𝑦 − 1 𝑦 − 1 𝜗 0000000000000 = 𝜗

Precise Analysis of log

82

To prove the 1 ulp error bound,

Need to construct more precise abstractions

Need to use floating-point exactness results

exact

We prove a bound of

0.583 ulps for log

SLIDE 83

/146

𝑠 𝑦 = 2 ⊗ 𝑦 ⊖ 255 128 ⊗ 1 2 ⊕ 255 128 ⊗ 1 2 ⊖ 1 4095 4096 ≤ 𝑦 < 1

More precise abstraction of 𝑠 𝑦 :

𝐵′𝜀 𝑦 = 𝑦 − 1 + 𝑦 − 1 𝜀′

Tighter bound on relative error of 𝑠 𝑦 :

max

. 4095 4096≤𝑦<1, 𝜀′ <𝜗

𝐵′𝜀 𝑦 − (𝑦 − 1) 𝑦 − 1 = max

. 4095 4096≤𝑦<1

𝑦 − 1 𝑦 − 1 𝜗 0000000000000 = 𝜗

Precise Analysis of log

83

To prove the 1 ulp error bound,

Need to construct more precise abstractions
Need to use floating-point exactness results

exact

We prove a bound of

0.583 ulps for log

SLIDE 84

/146

Sterbenz’s Theorem

Theorem [Sterbenz, 1973]

1 2 𝑏 ≤ 𝑐 ≤ 2𝑏

⟹ 𝑏 ⊖ 𝑐 = 𝑏 − 𝑐

Example:

84

SLIDE 85

/146

Sterbenz’s Theorem

Theorem [Sterbenz, 1973]

1 2 𝑏 ≤ 𝑐 ≤ 2𝑏

⟹ 𝑏 ⊖ 𝑐 = 𝑏 − 𝑐

Example: log for 𝑦 ∈

4095 4096 , 1 85

SLIDE 86

/146

Sterbenz’s Theorem

Theorem [Sterbenz, 1973]

1 2 𝑏 ≤ 𝑐 ≤ 2𝑏

⟹ 𝑏 ⊖ 𝑐 = 𝑏 − 𝑐

Example: log for 𝑦 ∈

4095 4096 , 1

Example: sin for 𝑦 ∈ 2𝜌 −

𝜌 64 , 2𝜌 + 𝜌 64

1. Compute 𝑧 = 𝑦 ⊖ 2𝜌
2. Return 𝑧 −

1 3! 𝑧3 + ⋯ + 1 9! 𝑧9

≈ 𝑡𝑗𝑜 𝑧

86

SLIDE 87

/146

Sterbenz’s Theorem

Theorem [Sterbenz, 1973]

1 2 𝑏 ≤ 𝑐 ≤ 2𝑏

⟹ 𝑏 ⊖ 𝑐 = 𝑏 − 𝑐

Example: log for 𝑦 ∈

4095 4096 , 1

Example: sin for 𝑦 ∈ 2𝜌 −

𝜌 64 , 2𝜌 + 𝜌 64

1. Compute 𝑧 = 𝑦 ⊖ 2𝜌
2. Return 𝑧 −

1 3! 𝑧3 + ⋯ + 1 9! 𝑧9

≈ 𝑡𝑗𝑜 𝑧

87

SLIDE 88

/146

Sterbenz’s Theorem

Theorem [Sterbenz, 1973]

1 2 𝑏 ≤ 𝑐 ≤ 2𝑏

⟹ 𝑏 ⊖ 𝑐 = 𝑏 − 𝑐

Example: log for 𝑦 ∈

4095 4096 , 1

Example: sin for 𝑦 ∈ 2𝜌 −

𝜌 64 , 2𝜌 + 𝜌 64

1. Compute 𝑧 = 𝑦 ⊖ 2𝜌
2. Return 𝑧 −

1 3! 𝑧3 + ⋯ + 1 9! 𝑧9

≈ 𝑡𝑗𝑜 𝑧

88

exact

SLIDE 89

/146

Sterbenz’s Theorem

Theorem [Sterbenz, 1973]

1 2 𝑏 ≤ 𝑐 ≤ 2𝑏

⟹ 𝑏 ⊖ 𝑐 = 𝑏 − 𝑐

Example: log for 𝑦 ∈

4095 4096 , 1

Example: sin for 𝑦 ∈ 2𝜌 −

𝜌 64 , 2𝜌 + 𝜌 64

1. Compute 𝑧 = 𝑦 ⊖ 2𝜌
2. Return 𝑧 −

1 3! 𝑧3 + ⋯ + 1 9! 𝑧9

≈ 𝑡𝑗𝑜 𝑧

89

exact

precision loss

f sin is small

SLIDE 90

/146

Sterbenz’s Theorem

How to apply the theorem automatically?

𝑓 ⊖ 𝑓′ ⊑ ? 𝑓 ⊖ 𝑓′

Construct: 𝑓 ⊑ 𝐵𝜀, 𝑓′ ⊑ 𝐵′

𝜀

Check:

min

𝑦, 𝜀

𝐵′𝜀 𝑦 − 1

2 𝐵𝜀 𝑦

≥ 0 and

max

𝑦, 𝜀

𝐵′𝜀 𝑦 − 2𝐵𝜀 𝑦 ≤ 0

⟹

1 2 𝐵𝜀 𝑦 ≤ 𝐵′ 𝜀 𝑦 ≤ 2𝐵𝜀 𝑦

for all 𝑦 ∈ 𝑌, all 𝜀𝑗 < 𝜗

Precondition of theorem:

⟹

1 2 𝑓 𝑦 ≤ 𝑓′ 𝑦 ≤ 2𝑓 𝑦

for all 𝑦 ∈ 𝑌

90

SLIDE 91

/146

Sterbenz’s Theorem

How to apply the theorem automatically?

𝑓 ⊖ 𝑓′ ⊑ ? 𝑓 ⊖ 𝑓′

Construct: 𝑓 ⊑ 𝐵𝜀, 𝑓′ ⊑ 𝐵′

𝜀

Check:

min

𝑦, 𝜀

𝐵′𝜀 𝑦 − 1

2 𝐵𝜀 𝑦

≥ 0 and

max

𝑦, 𝜀

𝐵′𝜀 𝑦 − 2𝐵𝜀 𝑦 ≤ 0

⟹

1 2 𝐵𝜀 𝑦 ≤ 𝐵′ 𝜀 𝑦 ≤ 2𝐵𝜀 𝑦

for all 𝑦 ∈ 𝑌, all 𝜀𝑗 < 𝜗

Precondition of theorem:

⟹

1 2 𝑓 𝑦 ≤ 𝑓′ 𝑦 ≤ 2𝑓 𝑦

for all 𝑦 ∈ 𝑌

91

SLIDE 92

/146

Sterbenz’s Theorem

How to apply the theorem automatically?

𝑓 ⊖ 𝑓′ ⊑ ? 𝑓 ⊖ 𝑓′

Construct: 𝑓 ⊑ 𝐵𝜀, 𝑓′ ⊑ 𝐵′

𝜀

Check:

min

𝑦, 𝜀

𝐵′𝜀 𝑦 − 1

2 𝐵𝜀 𝑦

≥ 0 and

max

𝑦, 𝜀

𝐵′𝜀 𝑦 − 2𝐵𝜀 𝑦 ≤ 0

⟹

1 2 𝐵𝜀 𝑦 ≤ 𝐵′ 𝜀 𝑦 ≤ 2𝐵𝜀 𝑦

for all 𝑦 ∈ 𝑌, all 𝜀𝑗 < 𝜗

Precondition of theorem:

⟹

1 2 𝑓 𝑦 ≤ 𝑓′ 𝑦 ≤ 2𝑓 𝑦

for all 𝑦 ∈ 𝑌

92

SLIDE 93

/146

Sterbenz’s Theorem

How to apply the theorem automatically?

𝑓 ⊖ 𝑓′ ⊑ ? 𝑓 ⊖ 𝑓′

Construct: 𝑓 ⊑ 𝐵𝜀, 𝑓′ ⊑ 𝐵′

𝜀

Check:

min

𝑦, 𝜀

𝐵′𝜀 𝑦 − 1

2 𝐵𝜀 𝑦

≥ 0 and

max

𝑦, 𝜀

𝐵′𝜀 𝑦 − 2𝐵𝜀 𝑦 ≤ 0

⟹

1 2 𝐵𝜀 𝑦 ≤ 𝐵′ 𝜀 𝑦 ≤ 2𝐵𝜀 𝑦

for all 𝑦 ∈ 𝑌, all 𝜀𝑗 < 𝜗

Precondition of theorem:

⟹

1 2 𝑓 𝑦 ≤ 𝑓′ 𝑦 ≤ 2𝑓 𝑦

for all 𝑦 ∈ 𝑌

93

SLIDE 94

/146

Sterbenz’s Theorem

How to apply the theorem automatically?

𝑓 ⊖ 𝑓′ ⊑ ? 𝑓 ⊖ 𝑓′

Construct: 𝑓 ⊑ 𝐵𝜀, 𝑓′ ⊑ 𝐵′

𝜀

Check:

min

𝑦, 𝜀

𝐵′𝜀 𝑦 − 1

2 𝐵𝜀 𝑦

≥ 0 and

max

𝑦, 𝜀

𝐵′𝜀 𝑦 − 2𝐵𝜀 𝑦 ≤ 0

⟹

1 2 𝐵𝜀 𝑦 ≤ 𝐵′ 𝜀 𝑦 ≤ 2𝐵𝜀 𝑦

for all 𝑦 ∈ 𝑌, all 𝜀𝑗 < 𝜗

Precondition of theorem:

⟹

1 2 𝑓 𝑦 ≤ 𝑓′ 𝑦 ≤ 2𝑓 𝑦

for all 𝑦 ∈ 𝑌

94

SLIDE 95

/146

Sterbenz’s Theorem

How to apply the theorem automatically?

𝑓 ⊖ 𝑓′ ⊑ ? 𝑓 ⊖ 𝑓′

Construct: 𝑓 ⊑ 𝐵𝜀, 𝑓′ ⊑ 𝐵′

𝜀

Check:

min

𝑦, 𝜀

𝐵′𝜀 𝑦 − 1

2 𝐵𝜀 𝑦

≥ 0 and

max

𝑦, 𝜀

𝐵′𝜀 𝑦 − 2𝐵𝜀 𝑦 ≤ 0

⟹

1 2 𝐵𝜀 𝑦 ≤ 𝐵′ 𝜀 𝑦 ≤ 2𝐵𝜀 𝑦

for all 𝑦 ∈ 𝑌, all 𝜀𝑗 < 𝜗

Precondition of theorem:

⟹

1 2 𝑓 𝑦 ≤ 𝑓′ 𝑦 ≤ 2𝑓 𝑦

for all 𝑦 ∈ 𝑌

95

SLIDE 96

/146

Sterbenz’s Theorem

How to apply the theorem automatically?

𝑓 ⊖ 𝑓′ ⊑ ? 𝑓 ⊖ 𝑓′

Construct: 𝑓 ⊑ 𝐵𝜀, 𝑓′ ⊑ 𝐵′

𝜀

Check:

min

𝑦, 𝜀

𝐵′𝜀 𝑦 − 1

2 𝐵𝜀 𝑦

≥ 0 and

max

𝑦, 𝜀

𝐵′𝜀 𝑦 − 2𝐵𝜀 𝑦 ≤ 0

⟹

1 2 𝐵𝜀 𝑦 ≤ 𝐵′ 𝜀 𝑦 ≤ 2𝐵𝜀 𝑦

for all 𝑦 ∈ 𝑌, all 𝜀𝑗 < 𝜗

Precondition of theorem:

⟹

1 2 𝑓 𝑦 ≤ 𝑓′ 𝑦 ≤ 2𝑓 𝑦

for all 𝑦 ∈ 𝑌

96

true

SLIDE 97

/146

Sterbenz’s Theorem

How to apply the theorem automatically?

𝑓 ⊖ 𝑓′ ⊑ ? 𝑓 ⊖ 𝑓′

Construct: 𝑓 ⊑ 𝐵𝜀, 𝑓′ ⊑ 𝐵′

𝜀

Check:

min

𝑦, 𝜀

𝐵′𝜀 𝑦 − 1

2 𝐵𝜀 𝑦

≥ 0 and

max

𝑦, 𝜀

𝐵′𝜀 𝑦 − 2𝐵𝜀 𝑦 ≤ 0

⟹

1 2 𝐵𝜀 𝑦 ≤ 𝐵′ 𝜀 𝑦 ≤ 2𝐵𝜀 𝑦

for all 𝑦 ∈ 𝑌, all 𝜀𝑗 < 𝜗

Precondition of theorem:

⟹

1 2 𝑓 𝑦 ≤ 𝑓′ 𝑦 ≤ 2𝑓 𝑦

for all 𝑦 ∈ 𝑌

97

true

SLIDE 98

/146

Sterbenz’s Theorem

How to apply the theorem automatically?

𝑓 ⊖ 𝑓′ ⊑ ? 𝑓 ⊖ 𝑓′

Construct: 𝑓 ⊑ 𝐵𝜀, 𝑓′ ⊑ 𝐵′

𝜀

Check:

min

𝑦, 𝜀

𝐵′𝜀 𝑦 − 1

2 𝐵𝜀 𝑦

≥ 0 and

max

𝑦, 𝜀

𝐵′𝜀 𝑦 − 2𝐵𝜀 𝑦 ≤ 0

⟹

1 2 𝐵𝜀 𝑦 ≤ 𝐵′ 𝜀 𝑦 ≤ 2𝐵𝜀 𝑦

for all 𝑦 ∈ 𝑌, all 𝜀𝑗 < 𝜗

Precondition of theorem:

⟹

1 2 𝑓 𝑦 ≤ 𝑓′ 𝑦 ≤ 2𝑓 𝑦

for all 𝑦 ∈ 𝑌

98

true

SLIDE 99

/146

Sterbenz’s Theorem

How to apply the theorem automatically?

𝑓 ⊖ 𝑓′ ⊑ ? 𝑓 ⊖ 𝑓′

Construct: 𝑓 ⊑ 𝐵𝜀, 𝑓′ ⊑ 𝐵′

𝜀

Check:

min

𝑦, 𝜀

𝐵′𝜀 𝑦 − 1

2 𝐵𝜀 𝑦

≥ 0 and

max

𝑦, 𝜀

𝐵′𝜀 𝑦 − 2𝐵𝜀 𝑦 ≤ 0

⟹

1 2 𝐵𝜀 𝑦 ≤ 𝐵′ 𝜀 𝑦 ≤ 2𝐵𝜀 𝑦

for all 𝑦 ∈ 𝑌, all 𝜀𝑗 < 𝜗

Precondition of theorem:

⟹

1 2 𝑓 𝑦 ≤ 𝑓′ 𝑦 ≤ 2𝑓 𝑦

for all 𝑦 ∈ 𝑌

99

true

can apply the theorem

SLIDE 100

/146

Sterbenz’s Theorem

How to apply the theorem automatically?

𝑓 ⊖ 𝑓′ ⊑ ? 𝑓 ⊖ 𝑓′

Construct: 𝑓 ⊑ 𝐵𝜀, 𝑓′ ⊑ 𝐵′

𝜀

Check:

min

𝑦, 𝜀

𝐵′𝜀 𝑦 − 1

2 𝐵𝜀 𝑦

≥ 0 and

max

𝑦, 𝜀

𝐵′𝜀 𝑦 − 2𝐵𝜀 𝑦 ≤ 0

⟹

1 2 𝐵𝜀 𝑦 ≤ 𝐵′ 𝜀 𝑦 ≤ 2𝐵𝜀 𝑦

for all 𝑦 ∈ 𝑌, all 𝜀𝑗 < 𝜗

Precondition of theorem:

⟹

1 2 𝑓 𝑦 ≤ 𝑓′ 𝑦 ≤ 2𝑓 𝑦

for all 𝑦 ∈ 𝑌

100

𝐵𝜀 𝑦 − 𝐵′

𝜀 𝑦

true

can apply the theorem

SLIDE 101

/146

Sterbenz’s Theorem

How to apply the theorem automatically?

𝑓 ⊖ 𝑓′ ⊑ ? 𝑓 ⊖ 𝑓′

Construct: 𝑓 ⊑ 𝐵𝜀, 𝑓′ ⊑ 𝐵′

𝜀

Check:

min

𝑦, 𝜀

𝐵′𝜀 𝑦 − 1

2 𝐵𝜀 𝑦

≥ 0 and

max

𝑦, 𝜀

𝐵′𝜀 𝑦 − 2𝐵𝜀 𝑦 ≤ 0

⟹

1 2 𝐵𝜀 𝑦 ≤ 𝐵′ 𝜀 𝑦 ≤ 2𝐵𝜀 𝑦

for all 𝑦 ∈ 𝑌, all 𝜀𝑗 < 𝜗

Precondition of theorem:

⟹

1 2 𝑓 𝑦 ≤ 𝑓′ 𝑦 ≤ 2𝑓 𝑦

for all 𝑦 ∈ 𝑌

101

𝐵𝜀 𝑦 − 𝐵′

𝜀 𝑦

true

can apply the theorem

SLIDE 102

/146

Dekker’s Theorem

Theorem [Dekker, 1971]

𝑏 ≥ 𝑐 and 𝑠 = 𝑏 ⊕ 𝑐 ⊖ 𝑏 ⊖ 𝑐 ⟹ 𝑠 = 𝑏 ⊕ 𝑐 − 𝑏 + 𝑐

102

SLIDE 103

/146

Dekker’s Theorem

Theorem [Dekker, 1971]

𝑏 ≥ 𝑐 and 𝑠 = 𝑏 ⊕ 𝑐 ⊖ 𝑏 ⊖ 𝑐 ⟹ 𝑠 = 𝑏 ⊕ 𝑐 − 𝑏 + 𝑐

103

SLIDE 104

/146

Dekker’s Theorem

Theorem [Dekker, 1971]

𝑏 ≥ 𝑐 and 𝑠 = 𝑏 ⊕ 𝑐 ⊖ 𝑏 ⊖ 𝑐 ⟹ 𝑠 = 𝑏 ⊕ 𝑐 − 𝑏 + 𝑐

Example: 𝑏 + 𝑐 + 𝑑, where 𝑏 ≥ 𝑐 ≥ 𝑑 > 0
Naïve ways: (𝑏 ⊕ 𝑐) ⊕ 𝑑, 𝑏 ⊕ 𝑐 ⊕ 𝑑 , ⋯
More accurate way:
1. Compute 𝑠 = 𝑏 ⊕ 𝑐 ⊖ 𝑏 ⊖ 𝑐
2. Return

𝑏 ⊕ 𝑐 ⊕ 𝑑 ⊖ 𝑠

104

SLIDE 105

/146

Dekker’s Theorem

Theorem [Dekker, 1971]

𝑏 ≥ 𝑐 and 𝑠 = 𝑏 ⊕ 𝑐 ⊖ 𝑏 ⊖ 𝑐 ⟹ 𝑠 = 𝑏 ⊕ 𝑐 − 𝑏 + 𝑐

Example: 𝑏 + 𝑐 + 𝑑, where 𝑏 ≥ 𝑐 ≥ 𝑑 > 0
Naïve ways: (𝑏 ⊕ 𝑐) ⊕ 𝑑, 𝑏 ⊕ 𝑐 ⊕ 𝑑 , ⋯

105

SLIDE 106

/146

Dekker’s Theorem

Theorem [Dekker, 1971]

𝑏 ≥ 𝑐 and 𝑠 = 𝑏 ⊕ 𝑐 ⊖ 𝑏 ⊖ 𝑐 ⟹ 𝑠 = 𝑏 ⊕ 𝑐 − 𝑏 + 𝑐

Example: 𝑏 + 𝑐 + 𝑑, where 𝑏 ≥ 𝑐 ≥ 𝑑 > 0
Naïve ways: (𝑏 ⊕ 𝑐) ⊕ 𝑑, 𝑏 ⊕ 𝑐 ⊕ 𝑑 , ⋯
More accurate way:
1. Compute 𝑠 = 𝑏 ⊕ 𝑐 ⊖ 𝑏 ⊖ 𝑐
2. Return

𝑏 ⊕ 𝑐 ⊕ 𝑑 ⊖ 𝑠

106

SLIDE 107

/146

Dekker’s Theorem

Theorem [Dekker, 1971]

𝑏 ≥ 𝑐 and 𝑠 = 𝑏 ⊕ 𝑐 ⊖ 𝑏 ⊖ 𝑐 ⟹ 𝑠 = 𝑏 ⊕ 𝑐 − 𝑏 + 𝑐

Example: 𝑏 + 𝑐 + 𝑑, where 𝑏 ≥ 𝑐 ≥ 𝑑 > 0
Naïve ways: (𝑏 ⊕ 𝑐) ⊕ 𝑑, 𝑏 ⊕ 𝑐 ⊕ 𝑑 , ⋯
More accurate way:
1. Compute 𝑠 = 𝑏 ⊕ 𝑐 ⊖ 𝑏 ⊖ 𝑐
2. Return

𝑏 ⊕ 𝑐 ⊕ 𝑑 ⊖ 𝑠

107

SLIDE 108

/146

Dekker’s Theorem

Theorem [Dekker, 1971]

𝑏 ≥ 𝑐 and 𝑠 = 𝑏 ⊕ 𝑐 ⊖ 𝑏 ⊖ 𝑐 ⟹ 𝑠 = 𝑏 ⊕ 𝑐 − 𝑏 + 𝑐

Example: 𝑏 + 𝑐 + 𝑑, where 𝑏 ≥ 𝑐 ≥ 𝑑 > 0
Naïve ways: (𝑏 ⊕ 𝑐) ⊕ 𝑑, 𝑏 ⊕ 𝑐 ⊕ 𝑑 , ⋯
More accurate way:
1. Compute 𝑠 = 𝑏 ⊕ 𝑐 ⊖ 𝑏 ⊖ 𝑐
2. Return

𝑏 ⊕ 𝑐 ⊕ 𝑑 ⊖ 𝑠

108

𝑏 + 𝑐 + 𝑠 =

SLIDE 109

/146

Dekker’s Theorem

Theorem [Dekker, 1971]

𝑏 ≥ 𝑐 and 𝑠 = 𝑏 ⊕ 𝑐 ⊖ 𝑏 ⊖ 𝑐 ⟹ 𝑠 = 𝑏 ⊕ 𝑐 − 𝑏 + 𝑐

Example: 𝑏 + 𝑐 + 𝑑, where 𝑏 ≥ 𝑐 ≥ 𝑑 > 0
Naïve ways: (𝑏 ⊕ 𝑐) ⊕ 𝑑, 𝑏 ⊕ 𝑐 ⊕ 𝑑 , ⋯
More accurate way:
1. Compute 𝑠 = 𝑏 ⊕ 𝑐 ⊖ 𝑏 ⊖ 𝑐
2. Return

𝑏 ⊕ 𝑐 ⊕ 𝑑 ⊖ 𝑠

109

𝑏 + 𝑐 + 𝑠 =

SLIDE 110

/146

Dekker’s Theorem

How to apply the theorem automatically?

𝑓 ⊕ 𝑓′ ⊖ 𝑓 ⊖ 𝑓′ ⊑ ? 𝑓 ⊕ 𝑓′ ⊖ 𝑓 ⊖ 𝑓′

Construct: 𝑓 ⊑ 𝐵𝜀, 𝑓′ ⊑ 𝐵′

𝜀, 𝑓 ⊕ 𝑓′ ⊑ 𝐵𝜀 𝑦 + 𝐵′ 𝜀 𝑦

1 + 𝜀

Check:

min

𝑦, 𝜀

𝐵𝜀 𝑦 ≥ max

𝑦, 𝜀

𝐵′𝜀 𝑦

⟹

𝐵𝜀 𝑦 ≥ 𝐵′𝜀 𝑦

for all 𝑦 ∈ 𝑌, all 𝜀𝑗 < 𝜗

Precondition of theorem:

⟹

𝑓 𝑦 ≥ 𝑓′ 𝑦

for all 𝑦 ∈ 𝑌

110

SLIDE 111

/146

Dekker’s Theorem

How to apply the theorem automatically?

𝑓 ⊕ 𝑓′ ⊖ 𝑓 ⊖ 𝑓′ ⊑ ? 𝑓 ⊕ 𝑓′ ⊖ 𝑓 ⊖ 𝑓′

Construct: 𝑓 ⊑ 𝐵𝜀, 𝑓′ ⊑ 𝐵′

𝜀, 𝑓 ⊕ 𝑓′ ⊑ 𝐵𝜀 𝑦 + 𝐵′ 𝜀 𝑦

1 + 𝜀

Check:

min

𝑦, 𝜀

𝐵𝜀 𝑦 ≥ max

𝑦, 𝜀

𝐵′𝜀 𝑦

⟹

𝐵𝜀 𝑦 ≥ 𝐵′𝜀 𝑦

for all 𝑦 ∈ 𝑌, all 𝜀𝑗 < 𝜗

Precondition of theorem:

⟹

𝑓 𝑦 ≥ 𝑓′ 𝑦

for all 𝑦 ∈ 𝑌

111

SLIDE 112

/146

Dekker’s Theorem

How to apply the theorem automatically?

𝑓 ⊕ 𝑓′ ⊖ 𝑓 ⊖ 𝑓′ ⊑ ? 𝑓 ⊕ 𝑓′ ⊖ 𝑓 ⊖ 𝑓′

Construct: 𝑓 ⊑ 𝐵𝜀, 𝑓′ ⊑ 𝐵′

𝜀, 𝑓 ⊕ 𝑓′ ⊑ 𝐵𝜀 𝑦 + 𝐵′ 𝜀 𝑦

1 + 𝜀

Check:

min

𝑦, 𝜀

𝐵𝜀 𝑦 ≥ max

𝑦, 𝜀

𝐵′𝜀 𝑦

⟹

𝐵𝜀 𝑦 ≥ 𝐵′𝜀 𝑦

for all 𝑦 ∈ 𝑌, all 𝜀𝑗 < 𝜗

Precondition of theorem:

⟹

𝑓 𝑦 ≥ 𝑓′ 𝑦

for all 𝑦 ∈ 𝑌

112

SLIDE 113

/146

Dekker’s Theorem

How to apply the theorem automatically?

𝑓 ⊕ 𝑓′ ⊖ 𝑓 ⊖ 𝑓′ ⊑ ? 𝑓 ⊕ 𝑓′ ⊖ 𝑓 ⊖ 𝑓′

Construct: 𝑓 ⊑ 𝐵𝜀, 𝑓′ ⊑ 𝐵′

𝜀, 𝑓 ⊕ 𝑓′ ⊑ 𝐵𝜀 𝑦 + 𝐵′ 𝜀 𝑦

1 + 𝜀

Check:

min

𝑦, 𝜀

𝐵𝜀 𝑦 ≥ max

𝑦, 𝜀

𝐵′𝜀 𝑦

⟹

𝐵𝜀 𝑦 ≥ 𝐵′𝜀 𝑦

for all 𝑦 ∈ 𝑌, all 𝜀𝑗 < 𝜗

Precondition of theorem:

⟹

𝑓 𝑦 ≥ 𝑓′ 𝑦

for all 𝑦 ∈ 𝑌

113

SLIDE 114

/146

Dekker’s Theorem

How to apply the theorem automatically?

𝑓 ⊕ 𝑓′ ⊖ 𝑓 ⊖ 𝑓′ ⊑ ? 𝑓 ⊕ 𝑓′ ⊖ 𝑓 ⊖ 𝑓′

Construct: 𝑓 ⊑ 𝐵𝜀, 𝑓′ ⊑ 𝐵′

𝜀, 𝑓 ⊕ 𝑓′ ⊑ 𝐵𝜀 𝑦 + 𝐵′ 𝜀 𝑦

1 + 𝜀

Check:

min

𝑦, 𝜀

𝐵𝜀 𝑦 ≥ max

𝑦, 𝜀

𝐵′𝜀 𝑦

⟹

𝐵𝜀 𝑦 ≥ 𝐵′𝜀 𝑦

for all 𝑦 ∈ 𝑌, all 𝜀𝑗 < 𝜗

Precondition of theorem:

⟹

𝑓 𝑦 ≥ 𝑓′ 𝑦

for all 𝑦 ∈ 𝑌

114

SLIDE 115

/146

Dekker’s Theorem

How to apply the theorem automatically?

𝑓 ⊕ 𝑓′ ⊖ 𝑓 ⊖ 𝑓′ ⊑ ? 𝑓 ⊕ 𝑓′ ⊖ 𝑓 ⊖ 𝑓′

Construct: 𝑓 ⊑ 𝐵𝜀, 𝑓′ ⊑ 𝐵′

𝜀, 𝑓 ⊕ 𝑓′ ⊑ 𝐵𝜀 𝑦 + 𝐵′ 𝜀 𝑦

1 + 𝜀

Check:

min

𝑦, 𝜀

𝐵𝜀 𝑦 ≥ max

𝑦, 𝜀

𝐵′𝜀 𝑦

⟹

𝐵𝜀 𝑦 ≥ 𝐵′𝜀 𝑦

for all 𝑦 ∈ 𝑌, all 𝜀𝑗 < 𝜗

Precondition of theorem:

⟹

𝑓 𝑦 ≥ 𝑓′ 𝑦

for all 𝑦 ∈ 𝑌

115

true

SLIDE 116

/146

Dekker’s Theorem

How to apply the theorem automatically?

𝑓 ⊕ 𝑓′ ⊖ 𝑓 ⊖ 𝑓′ ⊑ ? 𝑓 ⊕ 𝑓′ ⊖ 𝑓 ⊖ 𝑓′

Construct: 𝑓 ⊑ 𝐵𝜀, 𝑓′ ⊑ 𝐵′

𝜀, 𝑓 ⊕ 𝑓′ ⊑ 𝐵𝜀 𝑦 + 𝐵′ 𝜀 𝑦

1 + 𝜀

Check:

min

𝑦, 𝜀

𝐵𝜀 𝑦 ≥ max

𝑦, 𝜀

𝐵′𝜀 𝑦

⟹

𝐵𝜀 𝑦 ≥ 𝐵′𝜀 𝑦

for all 𝑦 ∈ 𝑌, all 𝜀𝑗 < 𝜗

Precondition of theorem:

⟹

𝑓 𝑦 ≥ 𝑓′ 𝑦

for all 𝑦 ∈ 𝑌

116

true

SLIDE 117

/146

Dekker’s Theorem

How to apply the theorem automatically?

𝑓 ⊕ 𝑓′ ⊖ 𝑓 ⊖ 𝑓′ ⊑ ? 𝑓 ⊕ 𝑓′ ⊖ 𝑓 ⊖ 𝑓′

Construct: 𝑓 ⊑ 𝐵𝜀, 𝑓′ ⊑ 𝐵′

𝜀, 𝑓 ⊕ 𝑓′ ⊑ 𝐵𝜀 𝑦 + 𝐵′ 𝜀 𝑦

1 + 𝜀

Check:

min

𝑦, 𝜀

𝐵𝜀 𝑦 ≥ max

𝑦, 𝜀

𝐵′𝜀 𝑦

⟹

𝐵𝜀 𝑦 ≥ 𝐵′𝜀 𝑦

for all 𝑦 ∈ 𝑌, all 𝜀𝑗 < 𝜗

Precondition of theorem:

⟹

𝑓 𝑦 ≥ 𝑓′ 𝑦

for all 𝑦 ∈ 𝑌

117

true

SLIDE 118

/146

Dekker’s Theorem

How to apply the theorem automatically?

𝑓 ⊕ 𝑓′ ⊖ 𝑓 ⊖ 𝑓′ ⊑ ? 𝑓 ⊕ 𝑓′ ⊖ 𝑓 ⊖ 𝑓′

Construct: 𝑓 ⊑ 𝐵𝜀, 𝑓′ ⊑ 𝐵′

𝜀, 𝑓 ⊕ 𝑓′ ⊑ 𝐵𝜀 𝑦 + 𝐵′ 𝜀 𝑦

1 + 𝜀

Check:

min

𝑦, 𝜀

𝐵𝜀 𝑦 ≥ max

𝑦, 𝜀

𝐵′𝜀 𝑦

⟹

𝐵𝜀 𝑦 ≥ 𝐵′𝜀 𝑦

for all 𝑦 ∈ 𝑌, all 𝜀𝑗 < 𝜗

Precondition of theorem:

⟹

𝑓 𝑦 ≥ 𝑓′ 𝑦

for all 𝑦 ∈ 𝑌

118

true

can apply the theorem

SLIDE 119

/146

Dekker’s Theorem

How to apply the theorem automatically?

𝑓 ⊕ 𝑓′ ⊖ 𝑓 ⊖ 𝑓′ ⊑ ? 𝑓 ⊕ 𝑓′ ⊖ 𝑓 ⊖ 𝑓′

Construct: 𝑓 ⊑ 𝐵𝜀, 𝑓′ ⊑ 𝐵′

𝜀, 𝑓 ⊕ 𝑓′ ⊑ 𝐵𝜀 𝑦 + 𝐵′ 𝜀 𝑦

1 + 𝜀

Check:

min

𝑦, 𝜀

𝐵𝜀 𝑦 ≥ max

𝑦, 𝜀

𝐵′𝜀 𝑦

⟹

𝐵𝜀 𝑦 ≥ 𝐵′𝜀 𝑦

for all 𝑦 ∈ 𝑌, all 𝜀𝑗 < 𝜗

Precondition of theorem:

⟹

𝑓 𝑦 ≥ 𝑓′ 𝑦

for all 𝑦 ∈ 𝑌

119

𝐵𝜀 𝑦 + 𝐵′

𝜀 𝑦

𝜀

true

can apply the theorem

SLIDE 120

/146

Dekker’s Theorem

How to apply the theorem automatically?

𝑓 ⊕ 𝑓′ ⊖ 𝑓 ⊖ 𝑓′ ⊑ ? 𝑓 ⊕ 𝑓′ ⊖ 𝑓 ⊖ 𝑓′

Construct: 𝑓 ⊑ 𝐵𝜀, 𝑓′ ⊑ 𝐵′

𝜀, 𝑓 ⊕ 𝑓′ ⊑ 𝐵𝜀 𝑦 + 𝐵′ 𝜀 𝑦

1 + 𝜀

Check:

min

𝑦, 𝜀

𝐵𝜀 𝑦 ≥ max

𝑦, 𝜀

𝐵′𝜀 𝑦

⟹

𝐵𝜀 𝑦 ≥ 𝐵′𝜀 𝑦

for all 𝑦 ∈ 𝑌, all 𝜀𝑗 < 𝜗

Precondition of theorem:

⟹

𝑓 𝑦 ≥ 𝑓′ 𝑦

for all 𝑦 ∈ 𝑌

120

𝐵𝜀 𝑦 + 𝐵′

𝜀 𝑦

𝜀

true

can apply the theorem

SLIDE 121

/146

More in Our Paper

More complicated exactness results
Their explicit statements
How to apply them automatically

121

SLIDE 122

/146

More in Our Paper

More complicated exactness results
Their explicit statements
How to apply them automatically
We check preconditions of exactness results
Abstractions
In practice, we cannot use naïve abstractions
because solving max

𝑦, 𝜀

⋯ 𝐵𝜀 𝑦 ⋯ is difficult

We apply our “compress” operation to abstractions
to reduce # of 𝜀 variables without losing much precision

122

SLIDE 123

/146

More in Our Paper

More complicated exactness results
Their explicit statements
How to apply them automatically
We check preconditions of exactness results
Abstractions
In practice, we cannot use naïve abstractions
because solving max

𝑦, 𝜀

⋯ 𝐵𝜀 𝑦 ⋯ is difficult

We apply our “compress” operation to abstractions
to reduce # of 𝜀 variables without losing much precision

123

Solve

ptimization problems

SLIDE 124

/146

More in Our Paper

More complicated exactness results
Their explicit statements
How to apply them automatically
We check preconditions of exactness results
Abstractions
In practice, we cannot use naïve abstractions
because solving max

𝑦, 𝜀

⋯ 𝐵𝜀 𝑦 ⋯ is difficult

We apply our “compress” operation to abstractions
to reduce # of 𝜀 variables without losing much precision

124

Solve

ptimization problems

SLIDE 125

/146

More in Our Paper

More complicated exactness results
Their explicit statements
How to apply them automatically
We check preconditions of exactness results
Abstractions
In practice, we cannot use naïve abstractions
because solving max

𝑦, 𝜀

⋯ 𝐵𝜀 𝑦 ⋯ is difficult

We apply our “compress” operation to abstractions
to reduce # of 𝜀 variables without losing much precision

125

Solve

ptimization problems

SLIDE 126

/146

More in Our Paper

More complicated exactness results
Their explicit statements
How to apply them automatically
We check preconditions of exactness results
Abstractions
In practice, we cannot use naïve abstractions
because solving max

𝑦, 𝜀

⋯ 𝐵𝜀 𝑦 ⋯ is difficult

We apply our “compress” operation to abstractions
to reduce # of 𝜀 variables without losing much precision

126

Solve

ptimization problems

SLIDE 127

/146

Case Studies

Benchmarks
sin, log, tan: Intel’s implementation of math.h
These are written directly in x86 assembly
Use bit-level / integer arithmetic operations (bit-shift, int add, ⋯)
Mixed in with floating-point operations
Eliminate them using the technique from [PLDI’16]

→ Results in multiple independent subproblems (on subintervals)

Apply our technique to the resulting floating-point expressions
Use Mathematica to solve optimization problems analytically

127

SLIDE 128

/146

Case Studies

Benchmarks
sin, log, tan: Intel’s implementation of math.h
These are written directly in x86 assembly
Use bit-level / integer arithmetic operations (bit-shift, int add, ⋯)
Mixed in with floating-point operations
Eliminate them using the technique from [PLDI’16]

→ Results in multiple independent subproblems (on subintervals)

Apply our technique to the resulting floating-point expressions
Use Mathematica to solve optimization problems analytically

128

SLIDE 129

/146

Case Studies

Benchmarks
sin, log, tan: Intel’s implementation of math.h
These are written directly in x86 assembly
Use bit-level / integer arithmetic operations (bit-shift, int add, ⋯)
Mixed in with floating-point operations
Eliminate them using the technique from [PLDI’16]

→ Results in multiple independent subproblems (on subintervals)

Apply our technique to the resulting floating-point expressions
Use Mathematica to solve optimization problems analytically

129

SLIDE 130

/146

Case Studies

Benchmarks
sin, log, tan: Intel’s implementation of math.h
These are written directly in x86 assembly
Use bit-level / integer arithmetic operations (bit-shift, int add, ⋯)
Mixed in with floating-point operations
Eliminate them using the technique from [PLDI’16]

→ Results in multiple independent subproblems (on subintervals)

Apply our technique to the resulting floating-point expressions
Use Mathematica to solve optimization problems analytically

130

SLIDE 131

/146

Results: sin on −𝜌, 𝜌

131

error bounds from [PLDI’16]

ur error bounds

1 ulp actual ulp errors ulp error (log scale) input value

SLIDE 132

/146

Results: sin on −𝜌, 𝜌

132

error bounds from [PLDI’16]

ur error bounds

1 ulp actual ulp errors ulp error (log scale) input value

SLIDE 133

/146

Results: sin on −𝜌, 𝜌

133

error bounds from [PLDI’16]

ur error bounds

1 ulp actual ulp errors ulp error (log scale) input value 0.530 ulps

SLIDE 134

/146

Results: sin on −𝜌, 𝜌

134

error bounds from [PLDI’16]

ur error bounds

1 ulp actual ulp errors ulp error (log scale) input value 0.530 ulps

SLIDE 135

/146

Results: log on 2−1022, 21024

135

error bounds from [PLDI’16]

ur error bounds

1 ulp actual ulp errors ulp error (log scale) input value (log scale)

SLIDE 136

/146

Results: log on 2−1022, 21024

136

error bounds from [PLDI’16]

ur error bounds

1 ulp actual ulp errors ulp error (log scale) input value (log scale)

SLIDE 137

/146

Results: log on 2−1022, 21024

137

error bounds from [PLDI’16]

ur error bounds

1 ulp actual ulp errors ulp error (log scale) input value (log scale) 0.583 ulps

SLIDE 138

/146

Results: log on 2−1022, 21024

138

error bounds from [PLDI’16]

ur error bounds

1 ulp actual ulp errors ulp error (log scale) input value (log scale)

461 hours on 16 cores
Highly parallelizable:
4 million subintervals
< 16 sec on 1 core per subinterval

0.583 ulps

SLIDE 139

/146

Results: log on 2−1022, 21024

139

error bounds from [PLDI’16]

ur error bounds

1 ulp actual ulp errors ulp error (log scale) input value (log scale)

461 hours on 16 cores
Highly parallelizable:
4 million subintervals
< 16 sec on 1 core per subinterval

0.583 ulps

SLIDE 140

/146

Results: tan on

13 128 , 17𝜌 64 , 17𝜌 64 , 𝜌 2

140

error bounds from [PLDI’16]

ur error bounds

1 ulp actual ulp errors ulp error (log scale) input value

SLIDE 141

/146

Results: tan on

13 128 , 17𝜌 64 , 17𝜌 64 , 𝜌 2

141

error bounds from [PLDI’16]

ur error bounds

1 ulp actual ulp errors ulp error (log scale) input value 0.595 ulps

SLIDE 142

/146

Results: tan on

13 128 , 17𝜌 64 , 17𝜌 64 , 𝜌 2

142

error bounds from [PLDI’16]

ur error bounds

1 ulp actual ulp errors ulp error (log scale) input value 13 ulps 0.595 ulps

SLIDE 143

/146

Results: tan on

13 128 , 17𝜌 64 , 17𝜌 64 , 𝜌 2

143

error bounds from [PLDI’16]

ur error bounds

1 ulp actual ulp errors ulp error (log scale) input value 13 ulps 0.595 ulps

SLIDE 144

/146

Results: tan on

13 128 , 17𝜌 64 , 17𝜌 64 , 𝜌 2

144

error bounds from [PLDI’16]

ur error bounds

1 ulp actual ulp errors ulp error (log scale) input value 13 ulps 0.595 ulps Our abstraction must be linear in 𝜀 variables → Results in imprecise abstractions

SLIDE 145

/146

Summary

A novel static analysis for verifying math.h implementations

automatically

A reduction of this verification task to optimization problems
Our technique verified some of Intel’s math.h implementations

automatically

145

SLIDE 146

/146

Summary

A novel static analysis for verifying math.h implementations

automatically

A reduction of this verification task to optimization problems
Our technique verified some of Intel’s math.h implementations

automatically

146