Floating p oin t v erication in HOL Ligh t: the exp onen - PDF document

Floating p oin t v eri�cation in HOL Ligh t: the exp onen tial function 1 Floating p oin t v eri�cation in HOL Ligh t: the exp onen tial function John Harrison Univ ersit y of Cam bridge � In tro duction � Floating p oin t correctness � Our implemen tation language � The algorithm � Outline of the HOL pro of � General conclusions John Harrison Univ ersit y of Cam bridge, Decem b er 1997

Floating p oin t v eri�cation in HOL Ligh t: the exp onen tial function 2 In tro duction � Floating p oin t algorithms are fairly small, but often complicated mathematically . � There ha v e b een errors in commercial systems, e.g. the P en tium FDIV bug in 1994. � In the case of transcenden tal functions it's di�cult ev en to sa y what correctness me ans . � V eri�cation using mo del c hec k ers is di�cult b ecause of the need for mathematical apparatus. � It can ev en b e di�cult using theorem pro v ers since not man y of them ha v e go o d theories of real n um b ers etc. John Harrison Univ ersit y of Cam bridge, Decem b er 1997

Floating p oin t v eri�cation in HOL Ligh t: the exp onen tial function 3 Floating p oin t correctness W e w an t to sp ecify the correctness according to the follo wing diagram: exp exp ( v ( a )) - v ( a ) v ( E X P ( a )) 6 6 v v - a E X P ( a ) E X P W e measure the di�erence b et w een v ( E X P ( a )) and exp ( v ( a )) in `units in the last place' of E X P ( a ). John Harrison Univ ersit y of Cam bridge, Decem b er 1997

Floating p oin t v eri�cation in HOL Ligh t: the exp onen tial function 4 Our implemen tation language This includes the follo wing constructs: c ommand = variable := expr ession | c ommand ; c ommand | if expr ession then c ommand else c ommand | if expr ession then c ommand | while expr ession do c ommand | do c ommand while expr ession | skip | f expr ession g W e de�ne a simple relational seman tics in HOL, and deriv e w eak est preconditions and total correctness rules. W e then pro v e total correctness via V C generation. The idea is that this language can b e formally link ed to C, V erilog, Handel, . . . John Harrison Univ ersit y of Cam bridge, Decem b er 1997

Floating p oin t v eri�cation in HOL Ligh t: the exp onen tial function 5 The algorithm The algorithm w e v erify is tak en from a pap er b y T ang in A CM T r ansactions on Mathematic al Softwar e , 1989. Similar tec hniques are widely used for �oating p oin t libraries, and, probably , for hardw are implemen tations. The algorithm relies on a table of precomputed constan ts. T ang's pap er giv es actual v alues as hex represen tations of IEEE n um b ers. The algorithm w orks in three phases: � P erform range reduction � Use p olynomial appro ximation � Reconstruct answ er using tables The correctness pro of re�ects this. John Harrison Univ ersit y of Cam bridge, Decem b er 1997

Floating p oin t v eri�cation in HOL Ligh t: the exp onen tial function 6 Co de for the algorithm if Isnan(X) then E := X else if X == Plus_infinity then E := Plus_infinity else if X == Minus_infinity then E := Plus_zero else if abs(X) > THRESHOLD_1 then if X > Plus_zero then E := Plus_infinity else E := Plus_zero else if abs(X) < THRESHOLD_2 then E := Plus_one + X else (N := INTRND(X * Inv_L); N2 := N % Int_32; N1 := N - N2; if abs(N) >= Int_2e9 then R1 := (X - Tofloat(N1) * L1) - Tofloat(N2) * L1 else R1 := X - Tofloat(N) * L1; R2 := Tofloat(--N) * L2; M := N1 / Int_32; J := N2; R := R1 + R2; Q := R * R * (A1 + R * A2); P := R1 + (R2 + Q); S := S_Lead(J) + S_Trail(J); E1 := S_Lead(J) + (S_Trail(J) + S * P); E := Scalb(E1,M) ) John Harrison Univ ersit y of Cam bridge, Decem b er 1997

Floating p oin t v eri�cation in HOL Ligh t: the exp onen tial function 7 Structure of the HOL pro of Real numbers / \ / \ / \ Programming / \ language IEEE spec Real analysis | / | | | / | | | / | | | / | Squarefree decomp & | / | Sturm's theorem | / | / | / | / Algorithm | / \ | / \ | / \ | / \ FP lemmas / \ | / \ | / \ | / \ | / \ | / Verification John Harrison Univ ersit y of Cam bridge, Decem b er 1997

Floating p oin t v eri�cation in HOL Ligh t: the exp onen tial function 8 Floating p oin t lemmas (1) W e de�ne the error error(x) resulting from rounding a real n um b er x to a �oating p oin t v alue. Because of the regular w a y in whic h the op erations are de�ned, all the op erations then relate to their abstract mathematical coun terparts according to the same pattern: |- Finite(a) ^ Finite(b) ^ abs(Val(a) + Val(b)) < threshold(float_format) = ) Finite(a + b) ^ (Val(a + b) = (Val(a) + Val(b)) + error(Val(a) + Val(b))) The comparisons are ev en more straigh tforw ard: |- Finite(a) ^ Finite(b) = ) (a < b = Val(a) < Val(b)) John Harrison Univ ersit y of Cam bridge, Decem b er 1997

Floating p oin t v eri�cation in HOL Ligh t: the exp onen tial function 9 Floating p oin t lemmas (2) W e ha v e sev eral lemmas quan tifying the error, e.g. |- abs(x) < threshold(float_format) ^ abs(x) < (&2 pow j / &2 pow 125) = ) abs(error(x)) <= &2 pow j / &2 pow 150 There are man y imp ortan t situations, ho w ev er, where the op erations are exact, b ecause the result is exactly represen table, e.g. subtraction of nearb y v alues with the same sign: |- Finite(a) ^ Finite(b) ^ &2 * abs(Val(a) - Val(b)) <= abs(Val(a)) = ) Finite(a - b) ^ (Val(a - b) = Val(a) - Val(b)) This is a classic result in �oating p oin t error analysis. John Harrison Univ ersit y of Cam bridge, Decem b er 1997

Floating p oin t v eri�cation in HOL Ligh t: the exp onen tial function 10 Informal error analysis T ang's error analysis translates quite directly in to HOL. One needs to: 1. Pro v e that clev er implemen tation tric ks ensure certain remainder terms are calculated exactly . This relies on cancellation, and the fact that pre-stored constan ts ha v e trailing zero es. 2. Pro v e that the p olynomial appro ximation ob eys the appropriate error b ounds. 3. Pro v e that the rounding errors when reconstructing the �nal answ er do not get to o large. In T ang's pap er, 1 is quite brief, 2 is dismissed in a few lines, while 3 is giv en a long and detailed pro of. John Harrison Univ ersit y of Cam bridge, Decem b er 1997

Floating p oin t v eri�cation in HOL Ligh t: the exp onen tial function 11 HOL error analysis In the HOL v ersion, this order of di�cult y is rev ersed! 1. The �rst part is not fundamen tally di�cult, but quite tric ky b ecause it in v olv es a lot of sp ecial cases and lo w-lev el pro ofs. 2. The second part in v olv es n umerical appro ximation, whic h needs a lot of w ork to translate in to a formal pro of (e.g. T a ylor series, Sturm's theorem . . . ). In fact T ang mak es a small mistak e here, though it do esn't a�ect the �nal result. 3. The last part is quite routine, and w e can program HOL to comp ose the rounding errors automatically . Actually , w e deriv e b etter b ounds than T ang do es since w e a v oid making simplifying assumptions to cut do wn the w ork. John Harrison Univ ersit y of Cam bridge, Decem b er 1997

Floating p oin t v eri�cation in HOL Ligh t: the exp onen tial function 12 The �nal result Under the v arious `de�nitional' assumptions, w e con�rm T ang's b ottom-line result: (Isnan(X) = ) Isnan(E)) ^ (X == Plus_infinity _ Finite(X) ^ exp(Val X) >= threshold(float_format) = ) E == Plus_infinity) ^ (X == Minus_infinity = ) E == Plus_zero) ^ (Finite(X) ^ exp(Val X) < threshold(float_format) = ) Isnormal(E) ^ abs(Val(E) - exp(Val X)) < (&54 / &100) * Ulp(E) _ (Isdenormal(E) _ Iszero(E)) ^ abs(Val(E) - exp(Val X)) < (&77 / &100) * Ulp(E)) In fact, this sp eci�cation is a bit more precise than T ang's, e.g. w e are explicit ab out the o v er�o w threshold. John Harrison Univ ersit y of Cam bridge, Decem b er 1997