Formal verification of numerical programs: from C annotated programs - - PowerPoint PPT Presentation

▶
formal verification of numerical programs from c
SMART_READER_LITE
LIVE PREVIEW

Formal verification of numerical programs: from C annotated programs - - PowerPoint PPT Presentation

Sep 08, 2023 342 likes •1.72k views

Third International Workshop on Numerical Software Verification Formal verification of numerical programs: from C annotated programs to Coq proofs Sylvie Boldo INRIA Saclay - Ile-de-France July 15th, 2010 Thanks to the organizers !

slide-1
SLIDE 1

Third International Workshop on Numerical Software Verification

Formal verification of numerical programs: from C annotated programs to Coq proofs

Sylvie Boldo INRIA Saclay - ˆ Ile-de-France July 15th, 2010

slide-2
SLIDE 2

Thanks to

the organizers !

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 2 / 41

slide-3
SLIDE 3

Thanks to

the organizers ! all collaborators of these works

◮ F. Cl´

ement

◮ J.-C. Filliˆ

atre

◮ G. Melquiond ◮ T. Nguyen Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 2 / 41

slide-4
SLIDE 4

Motivations

Numerical Software Verification

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 3 / 41

slide-5
SLIDE 5

Motivations

Numerical Software Verification ⇒ software with floating-point computations

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 3 / 41

slide-6
SLIDE 6

Floating-point number

This is only a string of bits. 11100011010010011110000111000000

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 4 / 41

slide-7
SLIDE 7

Floating-point number

This is only a string of bits. 11100011010010011110000111000000 We interpret it depending on the respective values of s (sign), e (exponent) and f (fraction). 1 11000110 10010011110000111000000 1 11000110 10010011110000111000000 s e f

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 4 / 41

slide-8
SLIDE 8

Floating-point number

We associate a real value : 1 11000110 10010011110000111000000 s e f ↓ ↓ ↓ (−1)s× 2e−B × 1 • f (−1)1× 2198−127 ×1.100100111100001110000002 −254 × 206727 ≈ −3.724 × 1021

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 5 / 41

slide-9
SLIDE 9

Floating-point number

We associate a real value : 1 11000110 10010011110000111000000 s e f ↓ ↓ ↓ (−1)s× 2e−B × 1 • f (−1)1× 2198−127 ×1.100100111100001110000002 −254 × 206727 ≈ −3.724 × 1021 except for the special values of e : ±0, ±∞, NaN, subnormals.

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 5 / 41

slide-10
SLIDE 10

Floating-point number repartition

R

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 6 / 41

slide-11
SLIDE 11

Floating-point number repartition

R subnormals

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 6 / 41

slide-12
SLIDE 12

Floating-point number repartition

R subnormals binade (common exponent)

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 6 / 41

slide-13
SLIDE 13

Floating-point operations

Thanks to the IEEE-754 standard, the computed results of +, −, ×, /, √ should be the same as if they were first computed with infinite precision and then rounded. ⇒ computations with 3 more bits (see J. Coonen)

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 7 / 41

slide-14
SLIDE 14

Floating-point operations

Thanks to the IEEE-754 standard, the computed results of +, −, ×, /, √ should be the same as if they were first computed with infinite precision and then rounded. ⇒ computations with 3 more bits (see J. Coonen) ⇒ mathematical properties such that : when a real value fits exactly in a floating-point number in a given format, then it is exactly computed.

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 7 / 41

slide-15
SLIDE 15

Motivations

Numerical Software Verification

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 8 / 41

slide-16
SLIDE 16

Motivations

Numerical Software Verification Critical C code ֒ → formal proof ⇒ high guarantee

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 8 / 41

slide-17
SLIDE 17

Related work

static analyzers

◮ Astr´

ee

◮ Fluctuat

specification languages

◮ JML

formal proofs about floating-point arithmetic

◮ trigonometric functions (HOL Light) ◮ verification of the FPU (ACL2) Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 9 / 41

slide-18
SLIDE 18

Motivations

C program

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 10 / 41

slide-19
SLIDE 19

Motivations

Annotated C program Human

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 10 / 41

slide-20
SLIDE 20

Motivations

Annotated C program Human Mathematical theorems Frama-C Jessie

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 10 / 41

slide-21
SLIDE 21

Motivations

Annotated C program Human Mathematical theorems Frama-C Jessie Coq ← Human

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 10 / 41

slide-22
SLIDE 22

Motivations

Annotated C program Human Mathematical theorems Frama-C Jessie Coq ← Human

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 10 / 41

slide-23
SLIDE 23

Motivations

Annotated C program Human Mathematical theorems Frama-C Jessie Coq ← Human

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 10 / 41

slide-24
SLIDE 24

Motivations

Annotated C program Human Mathematical theorems Frama-C Jessie Coq ← Human

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 10 / 41

slide-25
SLIDE 25

Motivations

Annotated C program Human Mathematical theorems Frama-C Jessie Coq ← Human

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 10 / 41

slide-26
SLIDE 26

Motivations

Annotated C program Human Mathematical theorems Frama-C Jessie Coq ← Human

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 10 / 41

slide-27
SLIDE 27

Motivations

Annotated C program Human Mathematical theorems Frama-C Jessie Coq ← Human

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 10 / 41

slide-28
SLIDE 28

Motivations

Annotated C program Human Mathematical theorems Frama-C Jessie Coq ← Human

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 10 / 41

slide-29
SLIDE 29

Motivations

Annotated C program Human Mathematical theorems Frama-C Jessie Proved theorems Coq ← Human

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 10 / 41

slide-30
SLIDE 30

Motivations

Annotated C program Human Mathematical theorems Frama-C Jessie Proved theorems Coq ← Human Program is correct w.r.t. its specification

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 10 / 41

slide-31
SLIDE 31

Plan

1

Motivations

2

Tools Formal proof Frama-C/Jessie/Why

3

Examples

4

Conclusions

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 11 / 41

slide-32
SLIDE 32

Formal proof

Certified formal proof The proof is checked in its deep details until the computer agrees with it. We often use formal proof checkers, meaning programs that only check a proof (they may also generate easy demonstrations). Therefore the checker is a very short program (de Bruijn criteria : the correctness of the system as a whole depends on the correctness of a very small ”kernel”).

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 12 / 41

slide-33
SLIDE 33

The Coq proof assistant (http://coq.inria.fr)

Based on the Curry-Howard isomorphism. (equivalence between proofs and λ-terms) Few automations. Comprehensive libraries, including on Z and R. Coq kernel mechanically checks each step of each proof. The method is to apply successively tactics (theorem application, rewriting, simplifications. . .) to transform or reduce the goal down to the hypotheses. The proof is handled starting from the conclusion.

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 13 / 41

slide-34
SLIDE 34

Coq formalization (by L. Th´ ery)

Float = pair of signed integers (mantissa, exponent)

(n, e) ∈ Z2

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 14 / 41

slide-35
SLIDE 35

Coq formalization (by L. Th´ ery)

Float = pair of signed integers (mantissa, exponent) associated to a real value.

(n, e) ∈ Z2 ֒ → n × βe ∈ R

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 14 / 41

slide-36
SLIDE 36

Coq formalization (by L. Th´ ery)

Float = pair of signed integers (mantissa, exponent) associated to a real value.

(n, e) ∈ Z2 ֒ → n × βe ∈ R

1.000102 E 4 → (1000102, −1)2 ֒ → 17 IEEE-754 significant of 754R real value ⇒ normal floats, subnormal floats, overflow. Many floats may represent the same real value, but we can exhibit a canon- ical representation.

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 14 / 41

slide-37
SLIDE 37

Example using Coq 8.2

Theorem Rle Fexp eq Zle : forall x y :float, (x <= y)%R -> Fexp x = Fexp y -> (Fnum x <= Fnum y)%Z. intros x y H’ H’0. apply le IZR. apply (Rle monotony contra exp radix) with (z := Fexp x); auto with real arith. pattern (Fexp x) at 2 in |- *; rewrite H’0; auto. Qed. With keywords, stating of the theorem, tactics and names of used theorems.

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 15 / 41

slide-38
SLIDE 38

Example using Coq 8.2

Theorem Rle Fexp eq Zle : forall x y :float, (x <= y)%R -> Fexp x = Fexp y -> (Fnum x <= Fnum y)%Z. intros x y H’ H’0. apply le IZR. apply (Rle monotony contra exp radix) with (z := Fexp x); auto with real arith. pattern (Fexp x) at 2 in |- *; rewrite H’0; auto. Qed. With keywords, stating of the theorem, tactics and names of used theorems.

Theorem (Rle Fexp eq Zle)

If two floats x = (nx, ex) and y = (ny, ey) verifies x ≤ y, and ex = ey, then nx ≤ ny.

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 15 / 41

slide-39
SLIDE 39

Plan

1

Motivations

2

Tools Formal proof Frama-C/Jessie/Why

3

Examples

4

Conclusions

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 16 / 41

slide-40
SLIDE 40

Frama-C/Jessie/Why

Frama-C is a framework dedicated to the analysis of the source code

  • f software written in C.

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 17 / 41

slide-41
SLIDE 41

Frama-C/Jessie/Why

Frama-C is a framework dedicated to the analysis of the source code

  • f software written in C.

Available plugins :

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 17 / 41

slide-42
SLIDE 42

Frama-C/Jessie/Why

Frama-C is a framework dedicated to the analysis of the source code

  • f software written in C.

Available plugins :

◮ value analysis Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 17 / 41

slide-43
SLIDE 43

Frama-C/Jessie/Why

Frama-C is a framework dedicated to the analysis of the source code

  • f software written in C.

Available plugins :

◮ value analysis ◮ Jessie, the deductive verification plug-in

(based on weakest precondition computation techniques)

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 17 / 41

slide-44
SLIDE 44

Frama-C/Jessie/Why

Frama-C is a framework dedicated to the analysis of the source code

  • f software written in C.

Available plugins :

◮ value analysis ◮ Jessie, the deductive verification plug-in

(based on weakest precondition computation techniques)

◮ . . . Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 17 / 41

slide-45
SLIDE 45

Frama-C/Jessie/Why

Frama-C is a framework dedicated to the analysis of the source code

  • f software written in C.

Available plugins :

◮ value analysis ◮ Jessie, the deductive verification plug-in

(based on weakest precondition computation techniques)

◮ . . .

Free softwares in CAML available at http://frama-c.com/ and http://why.lri.fr/.

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 17 / 41

slide-46
SLIDE 46

Frama-C/Jessie/Why

ACSL-annotated C program

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 18 / 41

slide-47
SLIDE 47

Frama-C/Jessie/Why

ACSL-annotated C program Frama-C/Jessie plug-in WHY verification condition generator Verification conditions

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 18 / 41

slide-48
SLIDE 48

Frama-C/Jessie/Why

ACSL-annotated C program Frama-C/Jessie plug-in WHY verification condition generator Verification conditions Automatic provers (Alt-Ergo,Gappa,CVC3,etc.) Interactive provers (Coq,PVS,etc.)

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 18 / 41

slide-49
SLIDE 49

ACSL

ANSI/ISO C Specification Language

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 19 / 41

slide-50
SLIDE 50

ACSL

ANSI/ISO C Specification Language behavioral specification language for C programs

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 19 / 41

slide-51
SLIDE 51

ACSL

ANSI/ISO C Specification Language behavioral specification language for C programs pre-conditions and post-conditions to functions (and which variables are modified).

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 19 / 41

slide-52
SLIDE 52

ACSL

ANSI/ISO C Specification Language behavioral specification language for C programs pre-conditions and post-conditions to functions (and which variables are modified). variants and invariants of the loops.

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 19 / 41

slide-53
SLIDE 53

ACSL

ANSI/ISO C Specification Language behavioral specification language for C programs pre-conditions and post-conditions to functions (and which variables are modified). variants and invariants of the loops. assertions

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 19 / 41

slide-54
SLIDE 54

ACSL

ANSI/ISO C Specification Language behavioral specification language for C programs pre-conditions and post-conditions to functions (and which variables are modified). variants and invariants of the loops. assertions In annotations, all computations are exact.

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 19 / 41

slide-55
SLIDE 55

ACSL and floating-point numbers

A floating-point number is a triple : the floating-point number, really computed by the program, x → xf floating-point part

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 20 / 41

slide-56
SLIDE 56

ACSL and floating-point numbers

A floating-point number is a triple : the floating-point number, really computed by the program, x → xf floating-point part the value that would have been obtained with exact computations, x → xe exact part

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 20 / 41

slide-57
SLIDE 57

ACSL and floating-point numbers

A floating-point number is a triple : the floating-point number, really computed by the program, x → xf floating-point part the value that would have been obtained with exact computations, x → xe exact part the value that we ideally wanted to compute x → xm model part

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 20 / 41

slide-58
SLIDE 58

ACSL and floating-point numbers

A floating-point number is a triple : the floating-point number, really computed by the program, x → xf floating-point part 1+x+x*x/2 the value that would have been obtained with exact computations, x → xe exact part 1 + x + x2

2

the value that we ideally wanted to compute x → xm model part exp(x)

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 20 / 41

slide-59
SLIDE 59

ACSL and floating-point numbers

A floating-point number is a triple : the floating-point number, really computed by the program, x → xf floating-point part 1+x+x*x/2 the value that would have been obtained with exact computations, x → xe exact part 1 + x + x2

2

the value that we ideally wanted to compute x → xm model part exp(x) ⇒ easy to split into method error and rounding error

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 20 / 41

slide-60
SLIDE 60

ACSL and floating-point numbers

A floating-point number is a triple : the floating-point number, really computed by the program, x → xf floating-point part 1+x+x*x/2 the value that would have been obtained with exact computations, x → xe exact part 1 + x + x2

2

the value that we ideally wanted to compute x → xm model part exp(x) ⇒ easy to split into method error and rounding error For a float f, we have macros such as \rounding error(f) and \exact(f), while f (as a real) is its floating-point value.

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 20 / 41

slide-61
SLIDE 61

Pragmas

Several pragmas corresponding to different formalization for floating-point numbers. defensive (default pragma) : IEEE roundings occur. We prove that no exceptional behavior may happen (Overflow, NaN creation. . .)

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 21 / 41

slide-62
SLIDE 62

Pragmas

Several pragmas corresponding to different formalization for floating-point numbers. defensive (default pragma) : IEEE roundings occur. We prove that no exceptional behavior may happen (Overflow, NaN creation. . .) math : all computations are exact.

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 21 / 41

slide-63
SLIDE 63

Pragmas

Several pragmas corresponding to different formalization for floating-point numbers. defensive (default pragma) : IEEE roundings occur. We prove that no exceptional behavior may happen (Overflow, NaN creation. . .) math : all computations are exact. full : IEEE roundings occur. Exceptional behaviors may happen.

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 21 / 41

slide-64
SLIDE 64

Pragmas

Several pragmas corresponding to different formalization for floating-point numbers. defensive (default pragma) : IEEE roundings occur. We prove that no exceptional behavior may happen (Overflow, NaN creation. . .) math : all computations are exact. full : IEEE roundings occur. Exceptional behaviors may happen. multi-rounding : we may have any hardware and compiler (80-bit extended registers, FMA)

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 21 / 41

slide-65
SLIDE 65

Plan

1

Motivations

2

Tools Formal proof Frama-C/Jessie/Why

3

Examples

4

Conclusions

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 22 / 41

slide-66
SLIDE 66

Examples

All examples use Frama-C Boron and Why 2.26.

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 23 / 41

slide-67
SLIDE 67

Examples

All examples use Frama-C Boron and Why 2.26. All proof obligations are proved using Coq. (except 2 inequalities in the last example).

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 23 / 41

slide-68
SLIDE 68

Examples

All examples use Frama-C Boron and Why 2.26. All proof obligations are proved using Coq. (except 2 inequalities in the last example). Code & proofs available on http://www.lri.fr/~sboldo/research.html.

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 23 / 41

slide-69
SLIDE 69

Sterbenz

Theorem (Sterbenz)

If x and y are FP numbers in a given precision such that y 2 ≤ x ≤ 2y, then x − y fits in a FP number in the same precision and is therefore computed without error.

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 24 / 41

slide-70
SLIDE 70

Sterbenz – program

/*@ requires y/2. <= x <= 2.*y; @ ensures \result == x-y; @*/ f l o a t Sterbenz ( f l o a t x , f l o a t y ) { return x−y ; }

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 25 / 41

slide-71
SLIDE 71

Sterbenz – program

/*@ requires y/2. <= x <= 2.*y; @ ensures \result == x-y; @*/ f l o a t Sterbenz ( f l o a t x , f l o a t y ) { return x−y ; } Exact subtraction

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 25 / 41

slide-72
SLIDE 72

Sterbenz – program

/*@ requires y/2. <= x <= 2.*y; @ ensures \result == x-y; @*/ f l o a t Sterbenz ( f l o a t x , f l o a t y ) { return x−y ; } 1 PO : exact subtraction

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 25 / 41

slide-73
SLIDE 73

Sterbenz – program

/*@ requires y/2. <= x <= 2.*y; @ ensures \result == x-y; @*/ f l o a t Sterbenz ( f l o a t x , f l o a t y ) { return x−y ; } 1 PO : exact subtraction 1 PO : no overflow

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 25 / 41

slide-74
SLIDE 74

Veltkamp/Dekker

Theorem (Veltkamp/Dekker)

Provided no Overflow and no Underflow occur, there is an algorithm computing the exact error of the multiplication using only FP operations.

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 26 / 41

slide-75
SLIDE 75

Veltkamp/Dekker

Theorem (Veltkamp/Dekker)

Provided no Overflow and no Underflow occur, there is an algorithm computing the exact error of the multiplication using only FP operations. Idea : split your floats in 2, multiply all the parts, add them in the correct order.

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 26 / 41

slide-76
SLIDE 76

Veltkamp : how to split a floating-point number

Let C = 227 + 1 for double precision numbers.

27 bits x 227 × x + p = ◦(x × C) q = ◦(x − p) x1 = p + q

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 27 / 41

slide-77
SLIDE 77

Dekker : how to get the error of the multiplication

x1 × y1 r1 = ◦(x × y) x1 × y2 t2 = t1 + x1 × y2 x2 × y1 x2 × y2 r2 = t3 + x2 × y2 t1 = x1 × y1 − r t3 = t2 + x2 × y1

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 28 / 41

slide-78
SLIDE 78

Veltkamp/Dekker – program

/*@ requires xy == \ round_double (\ NearestEven ,x*y) && @ \abs(x) <= 0x1.p995 && @ \abs(y) <= 0x1.p995 && @ \abs(x*y) <= 0x1.p1021; @ ensures ((x*y == 0 || 0x1.p -969 <= \abs(x*y)) @ ==> x*y == xy+\ result ); @*/ double Dekker ( double x , double y , double xy ) { double C, px , qx , hx , py , qy , hy , tx , ty , r2 ; i n t i ; [ . . . ] /*@ assert C == \pow (2. ,27) + 1. */ px=x∗C ; qx=x−px ; hx=px+qx ; tx=x−hx ; py=y∗C ; qy=y−py ; hy=py+qy ; ty=y−hy ; r2= −xy+hx∗hy ; r2+=hx∗ ty ; r2+=hy∗ tx ; r2+=tx ∗ ty ; return r2 ; }

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 29 / 41

slide-79
SLIDE 79

Veltkamp/Dekker – program

Split x and y

/*@ requires xy == \ round_double (\ NearestEven ,x*y) && @ \abs(x) <= 0x1.p995 && @ \abs(y) <= 0x1.p995 && @ \abs(x*y) <= 0x1.p1021; @ ensures ((x*y == 0 || 0x1.p -969 <= \abs(x*y)) @ ==> x*y == xy+\ result ); @*/ double Dekker ( double x , double y , double xy ) { double C, px , qx , hx , py , qy , hy , tx , ty , r2 ; i n t i ; [ . . . ] /*@ assert C == \pow (2. ,27) + 1. */ px=x∗C ; qx=x−px ; hx=px+qx ; tx=x−hx ; py=y∗C ; qy=y−py ; hy=py+qy ; ty=y−hy ; r2= −xy+hx∗hy ; r2+=hx∗ ty ; r2+=hy∗ tx ; r2+=tx ∗ ty ; return r2 ; }

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 29 / 41

slide-80
SLIDE 80

Veltkamp/Dekker – program

Multiply all halves and add all the results

/*@ requires xy == \ round_double (\ NearestEven ,x*y) && @ \abs(x) <= 0x1.p995 && @ \abs(y) <= 0x1.p995 && @ \abs(x*y) <= 0x1.p1021; @ ensures ((x*y == 0 || 0x1.p -969 <= \abs(x*y)) @ ==> x*y == xy+\ result ); @*/ double Dekker ( double x , double y , double xy ) { double C, px , qx , hx , py , qy , hy , tx , ty , r2 ; i n t i ; [ . . . ] /*@ assert C == \pow (2. ,27) + 1. */ px=x∗C ; qx=x−px ; hx=px+qx ; tx=x−hx ; py=y∗C ; qy=y−py ; hy=py+qy ; ty=y−hy ; r2= −xy+hx∗hy ; r2+=hx∗ ty ; r2+=hy∗ tx ; r2+=tx ∗ ty ; return r2 ; }

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 29 / 41

slide-81
SLIDE 81

Veltkamp/Dekker – program

xy = ◦(xy)

/*@ requires xy == \ round_double (\ NearestEven ,x*y) && @ \abs(x) <= 0x1.p995 && @ \abs(y) <= 0x1.p995 && @ \abs(x*y) <= 0x1.p1021; @ ensures ((x*y == 0 || 0x1.p -969 <= \abs(x*y)) @ ==> x*y == xy+\ result ); @*/ double Dekker ( double x , double y , double xy ) { double C, px , qx , hx , py , qy , hy , tx , ty , r2 ; i n t i ; [ . . . ] /*@ assert C == \pow (2. ,27) + 1. */ px=x∗C ; qx=x−px ; hx=px+qx ; tx=x−hx ; py=y∗C ; qy=y−py ; hy=py+qy ; ty=y−hy ; r2= −xy+hx∗hy ; r2+=hx∗ ty ; r2+=hy∗ tx ; r2+=tx ∗ ty ; return r2 ; }

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 29 / 41

slide-82
SLIDE 82

Veltkamp/Dekker – program

Overflow

/*@ requires xy == \ round_double (\ NearestEven ,x*y) && @ \abs(x) <= 0x1.p995 && @ \abs(y) <= 0x1.p995 && @ \abs(x*y) <= 0x1.p1021; @ ensures ((x*y == 0 || 0x1.p -969 <= \abs(x*y)) @ ==> x*y == xy+\ result ); @*/ double Dekker ( double x , double y , double xy ) { double C, px , qx , hx , py , qy , hy , tx , ty , r2 ; i n t i ; [ . . . ] /*@ assert C == \pow (2. ,27) + 1. */ px=x∗C ; qx=x−px ; hx=px+qx ; tx=x−hx ; py=y∗C ; qy=y−py ; hy=py+qy ; ty=y−hy ; r2= −xy+hx∗hy ; r2+=hx∗ ty ; r2+=hy∗ tx ; r2+=tx ∗ ty ; return r2 ; }

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 29 / 41

slide-83
SLIDE 83

Veltkamp/Dekker – program

If no Underflow

/*@ requires xy == \ round_double (\ NearestEven ,x*y) && @ \abs(x) <= 0x1.p995 && @ \abs(y) <= 0x1.p995 && @ \abs(x*y) <= 0x1.p1021; @ ensures ((x*y == 0 || 0x1.p -969 <= \abs(x*y)) @ ==> x*y == xy+\ result ); @*/ double Dekker ( double x , double y , double xy ) { double C, px , qx , hx , py , qy , hy , tx , ty , r2 ; i n t i ; [ . . . ] /*@ assert C == \pow (2. ,27) + 1. */ px=x∗C ; qx=x−px ; hx=px+qx ; tx=x−hx ; py=y∗C ; qy=y−py ; hy=py+qy ; ty=y−hy ; r2= −xy+hx∗hy ; r2+=hx∗ ty ; r2+=hy∗ tx ; r2+=tx ∗ ty ; return r2 ; }

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 29 / 41

slide-84
SLIDE 84

Veltkamp/Dekker – program

Exact error of ⊗

/*@ requires xy == \ round_double (\ NearestEven ,x*y) && @ \abs(x) <= 0x1.p995 && @ \abs(y) <= 0x1.p995 && @ \abs(x*y) <= 0x1.p1021; @ ensures ((x*y == 0 || 0x1.p -969 <= \abs(x*y)) @ ==> x*y == xy+\ result ); @*/ double Dekker ( double x , double y , double xy ) { double C, px , qx , hx , py , qy , hy , tx , ty , r2 ; i n t i ; [ . . . ] /*@ assert C == \pow (2. ,27) + 1. */ px=x∗C ; qx=x−px ; hx=px+qx ; tx=x−hx ; py=y∗C ; qy=y−py ; hy=py+qy ; ty=y−hy ; r2= −xy+hx∗hy ; r2+=hx∗ ty ; r2+=hy∗ tx ; r2+=tx ∗ ty ; return r2 ; }

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 29 / 41

slide-85
SLIDE 85

Accurate discriminant

It is pretty hard to compute b2 − ac accurately.

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 30 / 41

slide-86
SLIDE 86

Accurate discriminant

It is pretty hard to compute b2 − ac accurately.

Theorem (Kahan)

Provided no Overflow and no Underflow occur, there is an algorithm computing the b2 − a ∗ c within 2 ulps.

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 30 / 41

slide-87
SLIDE 87

Accurate discriminant – program

/*@ requires @ (b==0. || 0x1.p -916 <= \abs(b*b)) && @ (a*c==0. || 0x1.p -916 <= \abs(a*c)) && @ \abs(b) <= 0x1.p510 && @ \abs(a) <= 0x1.p995 && \abs(c) <= 0x1.p995 && @ \abs(a*c) <= 0x1.p1021; @ ensures \result ==0. @ || \abs (\ result -(b*b-a*c)) <= 2.* ulp (\ result ); @ */ double d i s c r i m i n a n t ( double a , double b , double c ) { double p , q , d , dp , dq ; p=b∗b ; q=a∗c ; i f ( p+q <= 3∗ f a b s (p−q )) d=p−q ; e l s e { dp=Dekker (b , b , p ) ; dq=Dekker ( a , c , q ) ; d=(p−q)+(dp−dq ) ; } return d ; }

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 31 / 41

slide-88
SLIDE 88

Accurate discriminant – program

/*@ requires @ (b==0. || 0x1.p -916 <= \abs(b*b)) && @ (a*c==0. || 0x1.p -916 <= \abs(a*c)) && @ \abs(b) <= 0x1.p510 && @ \abs(a) <= 0x1.p995 && \abs(c) <= 0x1.p995 && @ \abs(a*c) <= 0x1.p1021; @ ensures \result ==0. @ || \abs (\ result -(b*b-a*c)) <= 2.* ulp (\ result ); @ */ double d i s c r i m i n a n t ( double a , double b , double c ) { double p , q , d , dp , dq ; p=b∗b ; q=a∗c ; i f ( p+q <= 3∗ f a b s (p−q )) d=p−q ; e l s e { dp=Dekker (b , b , p ) ; dq=Dekker ( a , c , q ) ; d=(p−q)+(dp−dq ) ; } return d ; }

Test whether ac ≈ b2

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 31 / 41

slide-89
SLIDE 89

Accurate discriminant – program

If ac ≈ b2, compute naively

/*@ requires @ (b==0. || 0x1.p -916 <= \abs(b*b)) && @ (a*c==0. || 0x1.p -916 <= \abs(a*c)) && @ \abs(b) <= 0x1.p510 && @ \abs(a) <= 0x1.p995 && \abs(c) <= 0x1.p995 && @ \abs(a*c) <= 0x1.p1021; @ ensures \result ==0. @ || \abs (\ result -(b*b-a*c)) <= 2.* ulp (\ result ); @ */ double d i s c r i m i n a n t ( double a , double b , double c ) { double p , q , d , dp , dq ; p=b∗b ; q=a∗c ; i f ( p+q <= 3∗ f a b s (p−q )) d=p−q ; e l s e { dp=Dekker (b , b , p ) ; dq=Dekker ( a , c , q ) ; d=(p−q)+(dp−dq ) ; } return d ; }

Test whether ac ≈ b2

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 31 / 41

slide-90
SLIDE 90

Accurate discriminant – program

If ac ≈ b2, compute accurately using errors of the multiplications

/*@ requires @ (b==0. || 0x1.p -916 <= \abs(b*b)) && @ (a*c==0. || 0x1.p -916 <= \abs(a*c)) && @ \abs(b) <= 0x1.p510 && @ \abs(a) <= 0x1.p995 && \abs(c) <= 0x1.p995 && @ \abs(a*c) <= 0x1.p1021; @ ensures \result ==0. @ || \abs (\ result -(b*b-a*c)) <= 2.* ulp (\ result ); @ */ double d i s c r i m i n a n t ( double a , double b , double c ) { double p , q , d , dp , dq ; p=b∗b ; q=a∗c ; i f ( p+q <= 3∗ f a b s (p−q )) d=p−q ; e l s e { dp=Dekker (b , b , p ) ; dq=Dekker ( a , c , q ) ; d=(p−q)+(dp−dq ) ; } return d ; }

Test whether ac ≈ b2

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 31 / 41

slide-91
SLIDE 91

Accurate discriminant – program

Underflow

/*@ requires @ (b==0. || 0x1.p -916 <= \abs(b*b)) && @ (a*c==0. || 0x1.p -916 <= \abs(a*c)) && @ \abs(b) <= 0x1.p510 && @ \abs(a) <= 0x1.p995 && \abs(c) <= 0x1.p995 && @ \abs(a*c) <= 0x1.p1021; @ ensures \result ==0. @ || \abs (\ result -(b*b-a*c)) <= 2.* ulp (\ result ); @ */ double d i s c r i m i n a n t ( double a , double b , double c ) { double p , q , d , dp , dq ; p=b∗b ; q=a∗c ; i f ( p+q <= 3∗ f a b s (p−q )) d=p−q ; e l s e { dp=Dekker (b , b , p ) ; dq=Dekker ( a , c , q ) ; d=(p−q)+(dp−dq ) ; } return d ; }

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 31 / 41

slide-92
SLIDE 92

Accurate discriminant – program

Overflow

/*@ requires @ (b==0. || 0x1.p -916 <= \abs(b*b)) && @ (a*c==0. || 0x1.p -916 <= \abs(a*c)) && @ \abs(b) <= 0x1.p510 && @ \abs(a) <= 0x1.p995 && \abs(c) <= 0x1.p995 && @ \abs(a*c) <= 0x1.p1021; @ ensures \result ==0. @ || \abs (\ result -(b*b-a*c)) <= 2.* ulp (\ result ); @ */ double d i s c r i m i n a n t ( double a , double b , double c ) { double p , q , d , dp , dq ; p=b∗b ; q=a∗c ; i f ( p+q <= 3∗ f a b s (p−q )) d=p−q ; e l s e { dp=Dekker (b , b , p ) ; dq=Dekker ( a , c , q ) ; d=(p−q)+(dp−dq ) ; } return d ; }

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 31 / 41

slide-93
SLIDE 93

Accurate discriminant – program

2 ulps

/*@ requires @ (b==0. || 0x1.p -916 <= \abs(b*b)) && @ (a*c==0. || 0x1.p -916 <= \abs(a*c)) && @ \abs(b) <= 0x1.p510 && @ \abs(a) <= 0x1.p995 && \abs(c) <= 0x1.p995 && @ \abs(a*c) <= 0x1.p1021; @ ensures \result ==0. @ || \abs (\ result -(b*b-a*c)) <= 2.* ulp (\ result ); @ */ double d i s c r i m i n a n t ( double a , double b , double c ) { double p , q , d , dp , dq ; p=b∗b ; q=a∗c ; i f ( p+q <= 3∗ f a b s (p−q )) d=p−q ; e l s e { dp=Dekker (b , b , p ) ; dq=Dekker ( a , c , q ) ; d=(p−q)+(dp−dq ) ; } return d ; }

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 31 / 41

slide-94
SLIDE 94

Accurate discriminant – program

/*@ requires @ (b==0. || 0x1.p -916 <= \abs(b*b)) && @ (a*c==0. || 0x1.p -916 <= \abs(a*c)) && @ \abs(b) <= 0x1.p510 && @ \abs(a) <= 0x1.p995 && \abs(c) <= 0x1.p995 && @ \abs(a*c) <= 0x1.p1021; @ ensures \result ==0. @ || \abs (\ result -(b*b-a*c)) <= 2.* ulp (\ result ); @ */ double d i s c r i m i n a n t ( double a , double b , double c ) { double p , q , d , dp , dq ; p=b∗b ; q=a∗c ; i f ( p+q <= 3∗ f a b s (p−q )) d=p−q ; e l s e { dp=Dekker (b , b , p ) ; dq=Dekker ( a , c , q ) ; d=(p−q)+(dp−dq ) ; } return d ; }

Function calls ⇒ pre-conditions to prove ⇒ post-conditions guaranteed

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 31 / 41

slide-95
SLIDE 95

Accurate discriminant – program

/*@ requires @ (b==0. || 0x1.p -916 <= \abs(b*b)) && @ (a*c==0. || 0x1.p -916 <= \abs(a*c)) && @ \abs(b) <= 0x1.p510 && @ \abs(a) <= 0x1.p995 && \abs(c) <= 0x1.p995 && @ \abs(a*c) <= 0x1.p1021; @ ensures \result ==0. @ || \abs (\ result -(b*b-a*c)) <= 2.* ulp (\ result ); @ */ double d i s c r i m i n a n t ( double a , double b , double c ) { double p , q , d , dp , dq ; p=b∗b ; q=a∗c ; i f ( p+q <= 3∗ f a b s (p−q )) d=p−q ; e l s e { dp=Dekker (b , b , p ) ; dq=Dekker ( a , c , q ) ; d=(p−q)+(dp−dq ) ; } return d ; }

In initial proof, test assumed correct ⇒ Additional proof when test is incorrect

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 31 / 41

slide-96
SLIDE 96

Wave equation resolution scheme

∂2u(x, t) ∂t2 − c2 ∂2u(x, t) ∂x2 = s(x, t)

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 32 / 41

slide-97
SLIDE 97

Wave equation resolution scheme

∂2u(x, t) ∂t2 − c2 ∂2u(x, t) ∂x2 = s(x, t) ֒ →

x t xj+1 xj−1 xj tk tk−1 tk−2

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 32 / 41

slide-98
SLIDE 98

Wave equation resolution scheme

∂2u(x, t) ∂t2 − c2 ∂2u(x, t) ∂x2 = s(x, t) ֒ →

x t xj+1 xj−1 xj tk tk−1 tk−2

uk

j − 2uk−1 j

+ uk−2

j

∆t2 − c2 uk−1

j+1 − 2uk−1 j

+ uk−1

j−1

∆x2 = sk−1

j

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 32 / 41

slide-99
SLIDE 99

Wave equation resolution scheme – program

double ∗∗ forward prop ( i n t ni , i n t nk , double dx , double dt , double v , double xs , double l ) { double ∗∗p ; i n t i , k ; double a1 , a , dp ; a1 = dt /dx∗v ; a = a1∗a1 ; [ . . . ] // i n i t i a l i z a t i o n s

  • f

p [ . . . ] [ 0 ] and p [ . . . ] [ 1 ] /* propagation = time loop */ /*@ loop invariant 1 <= k <= nk && analytic_error (p,ni ,ni ,k,a); @ loop variant nk -k; */ f o r ( k=1; k<nk ; k++) { p [ 0 ] [ k+1] = 0 . ; /* time iteration = space loop */ /*@ loop invariant 1 <= i <= ni && analytic_error (p,ni ,i-1,k+1,a) @ loop variant ni -i; */ f o r ( i =1; i <n i ; i ++) { dp = p [ i +1][ k ] − 2.∗ p [ i ] [ k ] + p [ i −1][ k ] ; p [ i ] [ k+1] = 2.∗ p [ i ] [ k ] − p [ i ] [ k−1] + a∗dp ; } p [ n i ] [ k+1] = 0 . ; } return p ; }

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 33 / 41

slide-100
SLIDE 100

Wave equation resolution scheme – program

Space loop Time loop

double ∗∗ forward prop ( i n t ni , i n t nk , double dx , double dt , double v , double xs , double l ) { double ∗∗p ; i n t i , k ; double a1 , a , dp ; a1 = dt /dx∗v ; a = a1∗a1 ; [ . . . ] // i n i t i a l i z a t i o n s

  • f

p [ . . . ] [ 0 ] and p [ . . . ] [ 1 ] /* propagation = time loop */ /*@ loop invariant 1 <= k <= nk && analytic_error (p,ni ,ni ,k,a); @ loop variant nk -k; */ f o r ( k=1; k<nk ; k++) { p [ 0 ] [ k+1] = 0 . ; /* time iteration = space loop */ /*@ loop invariant 1 <= i <= ni && analytic_error (p,ni ,i-1,k+1,a) @ loop variant ni -i; */ f o r ( i =1; i <n i ; i ++) { dp = p [ i +1][ k ] − 2.∗ p [ i ] [ k ] + p [ i −1][ k ] ; p [ i ] [ k+1] = 2.∗ p [ i ] [ k ] − p [ i ] [ k−1] + a∗dp ; } p [ n i ] [ k+1] = 0 . ; } return p ; }

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 33 / 41

slide-101
SLIDE 101

Wave equation resolution scheme – program

Space loop Time loop Loop invariant Loop invariant

double ∗∗ forward prop ( i n t ni , i n t nk , double dx , double dt , double v , double xs , double l ) { double ∗∗p ; i n t i , k ; double a1 , a , dp ; a1 = dt /dx∗v ; a = a1∗a1 ; [ . . . ] // i n i t i a l i z a t i o n s

  • f

p [ . . . ] [ 0 ] and p [ . . . ] [ 1 ] /* propagation = time loop */ /*@ loop invariant 1 <= k <= nk && analytic_error (p,ni ,ni ,k,a); @ loop variant nk -k; */ f o r ( k=1; k<nk ; k++) { p [ 0 ] [ k+1] = 0 . ; /* time iteration = space loop */ /*@ loop invariant 1 <= i <= ni && analytic_error (p,ni ,i-1,k+1,a) @ loop variant ni -i; */ f o r ( i =1; i <n i ; i ++) { dp = p [ i +1][ k ] − 2.∗ p [ i ] [ k ] + p [ i −1][ k ] ; p [ i ] [ k+1] = 2.∗ p [ i ] [ k ] − p [ i ] [ k−1] + a∗dp ; } p [ n i ] [ k+1] = 0 . ; } return p ; }

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 33 / 41

slide-102
SLIDE 102

Wave equation resolution scheme – program

Main computations

double ∗∗ forward prop ( i n t ni , i n t nk , double dx , double dt , double v , double xs , double l ) { double ∗∗p ; i n t i , k ; double a1 , a , dp ; a1 = dt /dx∗v ; a = a1∗a1 ; [ . . . ] // i n i t i a l i z a t i o n s

  • f

p [ . . . ] [ 0 ] and p [ . . . ] [ 1 ] /* propagation = time loop */ /*@ loop invariant 1 <= k <= nk && analytic_error (p,ni ,ni ,k,a); @ loop variant nk -k; */ f o r ( k=1; k<nk ; k++) { p [ 0 ] [ k+1] = 0 . ; /* time iteration = space loop */ /*@ loop invariant 1 <= i <= ni && analytic_error (p,ni ,i-1,k+1,a) @ loop variant ni -i; */ f o r ( i =1; i <n i ; i ++) { dp = p [ i +1][ k ] − 2.∗ p [ i ] [ k ] + p [ i −1][ k ] ; p [ i ] [ k+1] = 2.∗ p [ i ] [ k ] − p [ i ] [ k−1] + a∗dp ; } p [ n i ] [ k+1] = 0 . ; } return p ; }

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 33 / 41

slide-103
SLIDE 103

Wave equation resolution scheme – program

Main computations

double ∗∗ forward prop ( i n t ni , i n t nk , double dx , double dt , double v , double xs , double l ) { double ∗∗p ; i n t i , k ; double a1 , a , dp ; a1 = dt /dx∗v ; a = a1∗a1 ; [ . . . ] // i n i t i a l i z a t i o n s

  • f

p [ . . . ] [ 0 ] and p [ . . . ] [ 1 ] /* propagation = time loop */ /*@ loop invariant 1 <= k <= nk && analytic_error (p,ni ,ni ,k,a); @ loop variant nk -k; */ f o r ( k=1; k<nk ; k++) { p [ 0 ] [ k+1] = 0 . ; /* time iteration = space loop */ /*@ loop invariant 1 <= i <= ni && analytic_error (p,ni ,i-1,k+1,a) @ loop variant ni -i; */ f o r ( i =1; i <n i ; i ++) { dp = p [ i +1][ k ] − 2.∗ p [ i ] [ k ] + p [ i −1][ k ] ; p [ i ] [ k+1] = 2.∗ p [ i ] [ k ] − p [ i ] [ k−1] + a∗dp ; } p [ n i ] [ k+1] = 0 . ; } return p ; }

Accumulation of rounding errors

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 33 / 41

slide-104
SLIDE 104

Wave equation resolution scheme – rounding error

Interval arithmetic ⇒ pk

i has error 2k2−53.

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 34 / 41

slide-105
SLIDE 105

Wave equation resolution scheme – rounding error

Interval arithmetic ⇒ pk

i has error 2k2−53.

We define εk

i as the signed rounding error made at step (i, k).

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 34 / 41

slide-106
SLIDE 106

Wave equation resolution scheme – rounding error

Interval arithmetic ⇒ pk

i has error 2k2−53.

We define εk

i as the signed rounding error made at step (i, k).

The predicate analytic error(x,t) is defined in Coq as : For all steps (i, k) that are under (x, t), |εk

i | ≤ 78 × 2−52

pk

i − exact(pk i ) = k

  • l=0

l

  • j=−l

αl

j εk−l i+j ,

with known αl

j

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 34 / 41

slide-107
SLIDE 107

Wave equation resolution scheme – rounding error

Interval arithmetic ⇒ pk

i has error 2k2−53.

We define εk

i as the signed rounding error made at step (i, k).

The predicate analytic error(x,t) is defined in Coq as : For all steps (i, k) that are under (x, t), |εk

i | ≤ 78 × 2−52

pk

i − exact(pk i ) = k

  • l=0

l

  • j=−l

αl

j εk−l i+j ,

with known αl

j

  • pk

i − exact

  • pk

i

  • ≤ 85 × 2−53 × (k + 1) × (k + 2)

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 34 / 41

slide-108
SLIDE 108

Wave equation resolution scheme – proof

33 proof obligations for the behavior (assertions, loop invariants, post-conditions. . .)

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 35 / 41

slide-109
SLIDE 109

Wave equation resolution scheme – proof

33 proof obligations for the behavior (assertions, loop invariants, post-conditions. . .) 84 proof obligations for the safety (loop variants, Overflow, pointer dereferencing. . .)

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 35 / 41

slide-110
SLIDE 110

Wave equation resolution scheme – proof

33 proof obligations for the behavior (assertions, loop invariants, post-conditions. . .) 84 proof obligations for the safety (loop variants, Overflow, pointer dereferencing. . .) 2 admits corresponding to the boundedness of the exact(pk

i )

(by scheme properties)

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 35 / 41

slide-111
SLIDE 111

Wave equation resolution scheme – proof

33 proof obligations for the behavior (assertions, loop invariants, post-conditions. . .) 84 proof obligations for the safety (loop variants, Overflow, pointer dereferencing. . .) 2 admits corresponding to the boundedness of the exact(pk

i )

(by scheme properties) 26000 lines of Coq (including less than 3700 lines of proof) (Note that the method error proof was presented at ITP on July 11th)

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 35 / 41

slide-112
SLIDE 112

Plan

1

Motivations

2

Tools Formal proof Frama-C/Jessie/Why

3

Examples

4

Conclusions

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 36 / 41

slide-113
SLIDE 113

Conclusion : advantages

Very high guarantee

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 37 / 41

slide-114
SLIDE 114

Conclusion : advantages

Very high guarantee not only rounding errors :

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 37 / 41

slide-115
SLIDE 115

Conclusion : advantages

Very high guarantee not only rounding errors :

◮ all other errors such as pointer dereferencing or division by zero Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 37 / 41

slide-116
SLIDE 116

Conclusion : advantages

Very high guarantee not only rounding errors :

◮ all other errors such as pointer dereferencing or division by zero ◮ link with mathematical properties Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 37 / 41

slide-117
SLIDE 117

Conclusion : advantages

Very high guarantee not only rounding errors :

◮ all other errors such as pointer dereferencing or division by zero ◮ link with mathematical properties ◮ any property can be checked Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 37 / 41

slide-118
SLIDE 118

Conclusion : advantages

Very high guarantee not only rounding errors :

◮ all other errors such as pointer dereferencing or division by zero ◮ link with mathematical properties ◮ any property can be checked

expressive annotation language (as expressive as Coq) ⇒ exactly the specification you want

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 37 / 41

slide-119
SLIDE 119

Conclusion : limits (1/2)

long and tedious

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 38 / 41

slide-120
SLIDE 120

Conclusion : limits (1/2)

long and tedious ⇒ automations !

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 38 / 41

slide-121
SLIDE 121

Conclusion : limits (1/2)

long and tedious ⇒ automations ! for example Gappa (G. Melquiond)

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 38 / 41

slide-122
SLIDE 122

Conclusion : limits (1/2)

long and tedious ⇒ automations ! for example Gappa (G. Melquiond) Verification conditions Automatic provers (Alt-Ergo,Gappa,CVC3,etc.) Interactive provers (Coq,PVS,etc.) ⇒ Use automatic provers to prove part of the verification conditions

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 38 / 41

slide-123
SLIDE 123

Conclusion : limits (1/2)

long and tedious ⇒ automations ! for example Gappa (G. Melquiond) Verification conditions Automatic provers (Alt-Ergo,Gappa,CVC3,etc.) Interactive provers (Coq,PVS,etc.) ⇒ Use automatic provers to prove part of the verification conditions ⇒ Use Gappa inside Coq to ease proofs

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 38 / 41

slide-124
SLIDE 124

Conclusion : limits (2/2)

We assume all double operations are direct 64-bit roundings.

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 39 / 41

slide-125
SLIDE 125

Conclusion : limits (2/2)

We assume all double operations are direct 64-bit roundings. On recent processors, we have x86 extended registers (80-bit long) and FMA (◦(ax + b) with one single rounding).

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 39 / 41

slide-126
SLIDE 126

Conclusion : limits (2/2)

We assume all double operations are direct 64-bit roundings. On recent processors, we have x86 extended registers (80-bit long) and FMA (◦(ax + b) with one single rounding). How does we know how the program was compiled and what will be the result ?

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 39 / 41

slide-127
SLIDE 127

Conclusion : limits (2/2)

We assume all double operations are direct 64-bit roundings. On recent processors, we have x86 extended registers (80-bit long) and FMA (◦(ax + b) with one single rounding). How does we know how the program was compiled and what will be the result ? Solution 1 : cover all cases. The result of an operation is a real near the correct result (it covers, 64-bit, 80-bit, double roundings and all uses of FMA) pragma multi-rounding Less precise, but always correct !

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 39 / 41

slide-128
SLIDE 128

Conclusion : limits (2/2)

We assume all double operations are direct 64-bit roundings. On recent processors, we have x86 extended registers (80-bit long) and FMA (◦(ax + b) with one single rounding). How does we know how the program was compiled and what will be the result ? Solution 1 : cover all cases. The result of an operation is a real near the correct result (it covers, 64-bit, 80-bit, double roundings and all uses of FMA) pragma multi-rounding Less precise, but always correct ! Solution 2 : look into the assembly. . .

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 39 / 41

slide-129
SLIDE 129

Perspectives

How to find correct specifications ?

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 40 / 41

slide-130
SLIDE 130

Perspectives

How to find correct specifications ? ⇒ use other tools. . .

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 40 / 41

slide-131
SLIDE 131

Perspectives

How to find correct specifications ? ⇒ use other tools. . . What about the Coq library ?

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 40 / 41

slide-132
SLIDE 132

Perspectives

How to find correct specifications ? ⇒ use other tools. . . What about the Coq library ? PFF (Boldo, Th´ ery,

  • Rideau. . .)

Gappa (Melquiond)

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 40 / 41

slide-133
SLIDE 133

Perspectives

How to find correct specifications ? ⇒ use other tools. . . What about the Coq library ? PFF (Boldo, Th´ ery,

  • Rideau. . .)

Gappa (Melquiond) 1 theo

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 40 / 41

slide-134
SLIDE 134

Perspectives

How to find correct specifications ? ⇒ use other tools. . . What about the Coq library ? PFF (Boldo, Th´ ery,

  • Rideau. . .)

Gappa (Melquiond) 1 theo Flocq (Boldo,Melquiond)

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 40 / 41

slide-135
SLIDE 135

Perspectives

How to find correct specifications ? ⇒ use other tools. . . What about the Coq library ? PFF (Boldo, Th´ ery,

  • Rideau. . .)

Gappa (Melquiond) 1 theo Flocq (Boldo,Melquiond) PFF theos Gappa theos

Sylvie Boldo (INRIA) Formal verification of numerical programs July 15th, 2010 40 / 41

slide-136
SLIDE 136

Thank you for your attention

Tools :

◮ http://frama-c.com/ ◮ http://why.lri.fr/ ◮ http://coq.inria.fr/

Code & proofs :

◮ http://www.lri.fr/~sboldo/research.html.

Formal proofs about scientific computations :

◮ http://fost.saclay.inria.fr/