Fil-C on Nix

Filnix wiki

the formal campaign

rocq/theories/Ieee.v

IEEE 754 arithmetic with integers only, transcribed from ieee.pl: binary32 (s), binary64 (d) and the x87 80-bit extended format (x), rounding to nearest, ties to even.

Operands and results are unsigned bit patterns, as in interp.pl, never host floats. A finite value is an exact rational, and round_ratio is the one place where rounding happens, gradual underflow and overflow included. The decimal parser and printer (ieee_parse/3, ieee_decimal/3) are the frontend's and stay out.

Parts

Gallinaieee.pl
fmt, frac_bits, exp_bits, bias, integer_bitformat/4, integer_bit/2
decodedecode/4
round_ratioround_ratio/5
ieee_from_int, ieee_to_intthe same
ieee_neg, ieee_binary, ieee_convert, ieee_comparethe same
From Stdlib Require Import ZArith Bool
Open Scope Z_scope

Formats

The three float formats, ieee.pl's K: binary32 (FmtS, s), binary64 (FmtD, d) and x87 extended (FmtX, x).

28Inductive fmt: Type ≔
|FmtS
|FmtD
|FmtX

As format(K, Fraction, _, _): the bits of the fraction of the format k, without the integer bit, which only x stores explicitly.

33Definition frac_bitsk : fmt : Z ≔
match k with
FmtS
⇒
23
FmtD
⇒
52
FmtX
⇒
63

As format(K, _, Exponent, _): the bits of the exponent field.

35Definition exp_bitsk : fmt : Z ≔
match k with
FmtS
⇒
8
FmtD
⇒
11
FmtX
⇒
15

As format(K, _, _, Bias): the exponent bias, the stored exponent of 1.0.

38Definition biask : fmt : Z ≔
match k with
FmtS
⇒
127
FmtD
⇒
1023
FmtX
⇒
16383

As integer_bit(K, Bits): the explicit integer bit, if any; the exponent starts above it.

42Definition integer_bitk : fmt : Z ≔
match k with
FmtX
⇒
1
_
⇒
0

As exponent_shift(K, Sh): the position of the exponent field's lowest bit, above the fraction and the integer bit.

46Definition exponent_shiftk : fmt : Z ≔
frac_bitsk + integer_bitk

Decoding

What a bit pattern is, decode/4's T.

51Inductive fclass: Type ≔
|Inf
an infinity
|NaN
not a number
|Finite
n p : Z
the finite value n × 2^p

As decode(K, B, S, T): the sign S and the value T of the pattern b in the format k. A zero exponent is a subnormal (or zero), whose significand is all the bits below the exponent.

59Definition decodek : fmt b : Z : Z × fclass ≔
let f ≔ frac_bitsk in
let e ≔ exp_bitsk in
let sh ≔ exponent_shiftk in
let s ≔ Z.shiftrb sh + e in
let m ≔ Z.landb 2 ^ f − 1 in
let x ≔ Z.landZ.shiftrb sh 2 ^ e − 1 in
let t ≔
x =? 2 ^ e − 1
→
m =? 0
→
Inf
else
NaN
x =? 0
→
FiniteZ.landb 2 ^ sh − 1 1 − biask − f
else
Finite2 ^ f + m x − biask − f
in
s, t

As sign_bits(K, S, B): the pattern with only the sign s, 0 or 1; also the zero of that sign.

74Definition sign_bitsk : fmt s : Z : Z ≔
Z.shiftls exponent_shiftk + exp_bitsk

As infinity(K, S, B): the infinity of sign s (x's has its integer bit set).

78Definition infinityk : fmt s : Z : Z ≔
sign_bitsk s + Z.shiftl2 ^ exp_bitsk − 1 exponent_shiftk + Z.shiftlinteger_bitk frac_bitsk

As quiet(K, A, B): the pattern a with the top fraction bit set, which makes a NaN quiet.

84Definition quietk : fmt a : Z : Z ≔
Z.lora 2 ^ frac_bitsk − 1

As invalid(K, B): the default NaN, a quiet negative one, the result of an invalid operation.

88Definition invalidk : fmt : Z ≔
quietk infinityk 1

Rounding

As ilog(N, L): floor(log2(N)) for a positive n, and -1 for 0.

93Definition ilogn : Z : Z ≔
n ≤? 0
→
−1
else
Z.log2n

As scale(N, D, E, A, B): the ratio A/B that is n/d × 2^e, both parts integers.

97Definition scalen d e : Z : Z × Z ≔
0 ≤? e
→
Z.shiftln e, d
else
n, Z.shiftld −e

As round_even(N, D, Q): the positive quotient N/D rounded to nearest, a tie to the even neighbour. N and D are positive, so Prolog's // is Z.div (Z.quot is the same here, but its extraction divides bit by bit, which takes minutes on the 16000-bit ratios of an x denormal).

105Definition round_evenn d : Z : Z ≔
let q0 ≔ n / d in
let twice ≔ 2 × n mod d in
d <? twice || twice =? d && q0 mod 2 =? 1
→
q0 + 1
else
q0

As round_ratio(K, S, N, D, B): the pattern of format k and sign s nearest n/d, for n ≥ 0, ties to even. The exponent is the one of n/d, raised to the least normal one (so small values round as subnormals), and an overflow is an infinity.

114Definition round_ratiok : fmt s n d : Z : Z ≔
n =? 0
→
sign_bitsk s
else
let f ≔ frac_bitsk in
let e0 ≔ ilogn − ilogd in
let 'a0, b0 ≔ scalen d −e0 in
let e ≔
a0 <? b0
→
e0 − 1
else
e0
in
let emin ≔ 1 − biask in
let step ≔ Z.maxe emin − f in
let 'a, c ≔ scalen d −step in
let q ≔ round_evena c in
let 'q1, e1 ≔
2 ^ f + 1 ≤? q
→
Z.shiftrq 1, e + 1
else
q, Z.maxe emin
in
biask <? e1
→
infinityk s
else
let 'exp, mant ≔
q1 <? 2 ^ f
→
0, q1
else
e1 + biask, q1 − Z.shiftl1 − integer_bitk f
in
sign_bitsk s + Z.shiftlexp exponent_shiftk + mant

As finite_ratio(T, A, B): the finite value n × 2^p as a ratio of non-negative integers.

134Definition finite_ration p : Z : Z × Z ≔
scalen 1 p

As signed_ratio(S, T, N, D): the same ratio, with the sign s applied to its numerator.

138Definition signed_ratios n p : Z : Z × Z ≔
let 'a, d ≔ finite_ration p in
1 − 2 × s × a, d

Integers

As ieee_from_int(K, N, B): the pattern of format k nearest the integer n.

145Definition ieee_from_intk : fmt n : Z : Z ≔
round_ratiok
n <? 0
→
1
else
0
Z.absn
1

As ieee_to_int(K, B, N): the finite pattern b truncated toward zero; none for an infinity or NaN. The ratio's parts are non-negative, so // is Z.div here too.

151Definition ieee_to_intk : fmt b : Z : optionZ ≔
match decodek b with
s, Finiten p
⇒
let 'a, d ≔ finite_ration p in
Some1 − 2 × s × a / d
_
⇒
None

Arithmetic

As ieee_neg(K, A, B): the pattern a with its sign bit flipped, NaNs included.

161Definition ieee_negk : fmt a : Z : Z ≔
Z.lxora sign_bitsk 1

As x87_nan(A, B, R): the result of an x op on two NaNs a and b, by the x87's rule (Intel SDM vol. 1, 4.8.3.5): a quiet one wins over a signaling one, else the one with the larger significand, else the first, and the winner is quieted.

167Definition x87_nana b : Z : Z ≔
let qa ≔ Z.landZ.shiftra 62 1 in
let qb ≔ Z.landZ.shiftrb 62 1 in
let ma ≔ Z.landa 2 ^ 64 − 1 in
let mb ≔ Z.landb 2 ^ 64 − 1 in
let w ≔
qb <? qa
→
a
qa <? qb
→
b
ma <? mb
→
b
else
a
in
quietFmtX w

The four arithmetic operations: add, subtract, multiply, divide.

176Inductive fop: Type ≔
|FAdd
|FSub
|FMul
|FDiv

Whether a decoded value t is a zero, of either sign.

179Definition is_zerot : fclass : bool ≔
match t with
Finite0 _
⇒
true
_
⇒
false

Whether a decoded value t is an infinity.

183Definition is_inft : fclass : bool ≔
match t with
Inf
⇒
true
_
⇒
false

As binary(add, K, SA, TA, SB, TB, R): the sum of two non-NaN values, ta of sign sa and tb of sign sb. An exact zero sum is negative only when both operands are; infinities of opposite signs give the default NaN.

189Definition binary_addk : fmt sa : Z ta : fclass sb : Z tb : fclass : Z ≔
match ta, tb with
Inf, Inf
⇒
negbsa =? sb
→
invalidk
else
infinityk sa
Inf, _
⇒
infinityk sa
_, Inf
⇒
infinityk sb
Finitena pa, Finitenb pb
⇒
let 'a, c ≔ signed_ratiosa na pa in
let 'b, d ≔ signed_ratiosb nb pb in
let n ≔ a × d + b × c in
let s ≔
n <? 0
→
1
n =? 0 && sa =? 1 && sb =? 1
→
1
else
0
in
round_ratiok s Z.absn c × d
_, _
⇒
invalidk

As binary(Op, K, SA, TA, SB, TB, R): the op o on two non-NaN values. The invalid cases, zero times infinity, zero over zero and infinity over infinity, give the default NaN invalid.

206Definition binaryo : fop k : fmt sa : Z ta : fclass sb : Z tb : fclass : Z ≔
let s ≔ Z.lxorsa sb in
match o with
FAdd
⇒
binary_addk sa ta sb tb
FSub
⇒
binary_addk sa ta 1 − sb tb
FMul
⇒
is_infta && is_zerotb || is_inftb && is_zerota
→
invalidk
is_infta || is_inftb
→
infinityk s
else
match ta, tb with
Finitena pa, Finitenb pb
⇒
let 'a, c ≔ finite_rationa pa in
let 'b, d ≔ finite_rationb pb in
round_ratiok s a × b c × d
_, _
⇒
invalidk
FDiv
⇒
is_infta && is_inftb || is_zerota && is_zerotb
→
invalidk
is_infta
→
infinityk s
is_inftb
→
sign_bitsk s
is_zerotb
→
infinityk s
else
match ta, tb with
Finitena pa, Finitenb pb
⇒
let 'a, c ≔ finite_rationa pa in
let 'b, d ≔ finite_rationb pb in
round_ratiok s a × d c × b
_, _
⇒
invalidk

As ieee_binary(K, Op, A, B, R): the op o on the patterns a and b of format k. A NaN operand gives that NaN quieted, the first if both are, as x86's SSE does (x87_nan for x); otherwise binary on the values.

239Definition ieee_binaryk : fmt o : fop a b : Z : Z ≔
let 'sa, ta ≔ decodek a in
let 'sb, tb ≔ decodek b in
match ta, tb, k with
NaN, NaN, FmtX
⇒
x87_nana b
NaN, _, _
⇒
quietk a
_, NaN, _
⇒
quietk b
_, _, _
⇒
binaryo k sa ta sb tb

As ieee_convert(From, To, A, B): the pattern a of format from in the format to, rounded to nearest; a NaN keeps as much of its payload as fits, aligned at the top, and is quieted.

252Definition ieee_convertfrom to : fmt a : Z : Z ≔
match decodefrom a with
s, Inf
⇒
infinityto s
s, NaN
⇒
let payload ≔ Z.landa 2 ^ frac_bitsfrom − 1 in
let shift ≔ frac_bitsto − frac_bitsfrom in
let p ≔
0 ≤? shift
→
Z.shiftlpayload shift
else
Z.shiftrpayload −shift
in
quietto Z.lorinfinityto s p
s, Finiten p
⇒
let 'nn, d ≔ finite_ration p in
round_ratioto s nn d

As ieee_sqrt(K, A, B): the square root of the pattern a of format k, rounded to nearest, ties to even, as sqrtss, sqrtsd and fsqrt compute it. A NaN is quieted; +0, -0 and +inf are themselves; any other negative number, -inf included, gives the default NaN. For a finite n × 2^p, p made even, the root of n × 4^g with g = f + 5 has at least g bits, two more than rounding needs, and an inexact root lying strictly between q and q + 1 rounds as q + 1/2 does. Z.sqrt is isqrt/2, the floor of the root.

271Definition ieee_sqrtk : fmt a : Z : Z ≔
let 's, t ≔ decodek a in
match t with
NaN
⇒
quietk a
Finite0 _
⇒
a
Inf
⇒
s =? 1
→
invalidk
else
a
Finiten0 p0
⇒
s =? 1
→
invalidk
else
let 'n1, p ≔
Z.evenp0
→
n0, p0
else
n0 × 2, p0 − 1
in
let g ≔ frac_bitsk + 5 in
let m ≔ Z.shiftln1 2 × g in
let q ≔ Z.sqrtm in
let 'num, den0 ≔
q × q =? m
→
q, 1
else
2 × q + 1, 2
in
let 'num1, den ≔ scalenum den0 p / 2 − g in
round_ratiok 0 num1 den

Comparison

A float comparison's condition (condition/2): equal, not equal, the four orders, ordered (neither is a NaN) and unordered (one is).

293Inductive fcond: Type ≔
|FEq
|FNe
|FLt
|FLe
|FGt
|FGe
|FOrd
|FUnord

How two values compare: the standard order, or unordered when one is a NaN.

297Inductive order: Type ≔
|Below
|Same
|Above
|Unordered

The order of a Rocq comparison.

300Definition order_ofc : comparison : order ≔
match c with
Lt
⇒
Below
Eq
⇒
Same
Gt
⇒
Above

As condition(Cond, Order): whether the condition c holds of two values that compare as o.

305Definition conditionc : fcond o : order : bool ≔
match c, o with
FEq, Same
⇒
true
FNe, Below | Above | Unordered
⇒
true
FLt, Below
⇒
true
FLe, Below | Same
⇒
true
FGt, Above
⇒
true
FGe, Above | Same
⇒
true
FOrd, Below | Same | Above
⇒
true
FUnord, Unordered
⇒
true
_, _
⇒
false

As ieee_compare(K, Cond, A, B, R): 1 when a c b holds of the patterns a and b of format k, else 0. With a NaN only ne and uo hold; the two zeros are equal.

321Definition ieee_comparek : fmt c : fcond a b : Z : Z ≔
let 'sa, ta ≔ decodek a in
let 'sb, tb ≔ decodek b in
let o ≔
match ta, tb with
NaN, _ | _, NaN
⇒
Unordered
Inf, Inf
⇒
order_ofZ.comparesb sa
Inf, _
⇒
sa =? 1
→
Below
else
Above
_, Inf
⇒
sb =? 1
→
Above
else
Below
Finitena pa, Finitenb pb
⇒
let 'n, d ≔ signed_ratiosa na pa in
let 'm, e ≔ signed_ratiosb nb pb in
order_ofZ.comparen × e m × d
in
conditionc o
→
1
else
0