Fil-C on Nix

Filnix wiki

mbrock/filnix/nanofil

ieee.pl

IEEE 754 binary32 and binary64 arithmetic with integers only, rounding to nearest, ties to even.

Operands and results are unsigned bit patterns, never host floats. Finite intermediates are exact rationals, and round_ratio is the single rounding boundary, gradual underflow and overflow included. All public inputs must be concrete.

Predicates

K is s (binary32), d (binary64) or x (the x87 80-bit extended format of C's long double on x86-64: a sign, 15 exponent bits and a 64-bit significand whose integer bit is explicit, in the low 80 bits of the pattern). An x NaN operand follows the x87's rule instead of SSE's: when both operands are NaNs, a quiet one wins over a signaling one, else the one with the larger significand, else the first.

PredicateNotes
ieee_binaryK, Op, A, B, ROp is add, sub, mul or div
ieee_negK, A, Bflips the sign bit
ieee_convertFrom, To, A, Bbetween the two formats
ieee_from_intK, N, BN any integer
ieee_to_intK, B, Ntruncates; fails for nonfinite inputs, and the caller checks the destination integer's range
ieee_compareK, Cond, A, B, RR is 1 or 0; Cond is eq, ne, lt, le, gt, ge, o (ordered) or uo (unordered)
ieee_parseK, Chars, Ba decimal or C hexadecimal significand, without a C type suffix: the frontend selects K after removing the suffix
ieee_decimalK, B, Charsan exact decimal, not necessarily the shortest; fails for nonfinite inputs
useslists
usesdcgs
useserror

format(K, fraction bits, exponent bits, bias): the fraction excludes the integer bit, which x stores explicitly (integer_bit/2).

format
39s238127
40d52111023
41x631516383
integer_bitK, Bits

the bits of the explicit integer bit, if any; the exponent starts above it.

integer_bit
45s0
46d0
47x1
48exponent_shiftK, Sh ←
formatK, F, _, _
Sh is F + I
50decodeK, B, S, T ←
must_beinteger, B
formatK, F, E, Bias
S is B >> Sh + E
M is B /\ 1 << F−1
X is B >> Sh /\ 1 << E−1
Max is 1 << E−1
X =:= Max
→
M =:= 0
→
T = inf
else
T = nan
X =:= 0
→
P is 1−Bias−F
Sig is B /\ 1 << Sh−1
T = finiteSig, P
else
N is 1 << F + M
P is X−Bias−F
T = finiteN, P
58sign_bitsK, S, B ←
formatK, _, E, _
B is S << Sh + E
59infinityK, S, B ←
formatK, F, E, _
sign_bitsK, S, SB
B is SB + 1 << E−1 << Sh + I << F
62quietK, A, B ←
formatK, F, _, _
B is A \/ 1 << F−1
63invalidK, B ←
infinityK, 1, I
quietK, I, B

floor(log2(N)), for positive integers, without host logarithms.

66ilogN, L ←
ilogN, -1, L
67ilog0, L, L ←
!
68ilogN, A, L ←
N1 is N >> 1
A1 is A + 1
ilogN1, A1, L
70scaleN, D, E, A, B ←
E ≥ 0
→
A is N << E
B = D
else
A = N
B is D << −E

Rounding an exact positive quotient, including the halfway parity bit.

74round_evenN, D, Q ←
Q0 is N//D
R is N mod D
Twice is 2 × R
Twice > D
or
Twice =:= D
Q0 mod 2 =:= 1
!
Q is Q0 + 1
78round_evenN, D, Q ←
Q is N//D
80round_ratioK, S, 0, _, B ←
!
sign_bitsK, S, B
81round_ratioK, S, N, D, B ←
formatK, F, _, Bias
ilogN, LN
ilogD, LD
E0 is LN−LD
scaleN, D, −E0, A0, B0
A0 < B0
→
E is E0−1
else
E = E0
Emin is 1−Bias
Step is maxE, Emin−F
scaleN, D, −Step, A, C
round_evenA, C, Q
Q ≥ 1 << F + 1
→
Q1 is Q >> 1
E1 is E + 1
else
Q1 = Q
E1 is maxE, Emin
E1 > Bias
→
infinityK, S, B
Q1 < 1 << F
→
Exp = 0
Mant = Q1
else
Exp is E1 + Bias
Mant is Q1−1−I << F
E1 ≤ Bias
→
sign_bitsK, S, SB
B is SB + Exp << Sh + Mant
else
true
96finite_ratiofiniteN, E, A, B ←
scaleN, 1, E, A, B
97signed_ratioS, T, N, D ←
N is 1−2 × S × A
ieee_from_int+K, +N, -Bdet

B is the format-K bit pattern nearest the integer N, ties to even.

103ieee_from_intK, N, B ←
must_beinteger, N
N < 0
→
S = 1
else
S = 0
round_ratioK, S, A, 1, B

Trealla's abs/1 and unary minus wrap a computed native INT64_MIN; binary subtraction promotes it to a bignum on both hosts.

109magnitude_integerN, A ←
N < 0
→
A is 0−N
else
A = N
ieee_to_int+K, +B, -Nsemidet

N is the finite float B truncated toward zero. Fails for an infinity or NaN; the caller checks that N fits the destination integer.

115ieee_to_intK, B, N ←
decodeK, B, S, T
T = finite_, _
N is 1−2 × S × A//D
ieee_neg+K, +A, -Bdet

B is A with its sign bit flipped, NaNs included.

122ieee_negK, A, B ←
must_beinteger, A
sign_bitsK, 1, S
B is A xor S
ieee_binary+K, +Op, +A, +B, -Rdet

R is A Op B in format K, Op add, sub, mul or div, rounded to nearest, ties to even, with SSE's NaN rules (x87's for K = x, header).

129ieee_binaryK, Op, A, B, R ←
decodeK, A, SA, TA
decodeK, B, SB, TB
TA ≡ nan
TB ≡ nan
K ≡ x
→
x87_nanA, B, R
TA ≡ nan
→
quietK, A, R
TB ≡ nan
→
quietK, B, R
else
binaryOp, K, SA, TA, SB, TB, R

The x87 picks between two NaN operands (Intel SDM vol. 1, 4.8.3.5).

137x87_nanA, B, R ←
QA is A >> 62 /\ 1
QB is B >> 62 /\ 1
MA is A /\ 1 << 64−1
MB is B /\ 1 << 64−1
QA > QB
→
W = A
QB > QA
→
W = B
MB > MA
→
W = B
else
W = A
quietx, W, R
143binarysub, K, SA, TA, SB, TB, R ←
SC is 1−SB
binaryadd, K, SA, TA, SC, TB, R
145binaryadd, K, SA, TA, SB, TB, R ←
TA ≡ inf
TB ≡ inf
SA ≠ SB
→
TA ≡ inf
→
infinityK, SA, R
TB ≡ inf
→
infinityK, SB, R
else
signed_ratioSA, TA, A, C
signed_ratioSB, TB, B, D
N is A × D + B × C
Den is C × D
N < 0
→
S = 1
N =:= 0
SA =:= 1
SB =:= 1
→
S = 1
else
S = 0
round_ratioK, S, Abs, Den, R
155binarymul, K, SA, TA, SB, TB, R ←
S is SA xor SB
TA ≡ inf
TB = finite0, _
or
TB ≡ inf
TA = finite0, _
→
TA ≡ inf
or
TB ≡ inf
→
infinityK, S, R
else
finite_ratioTA, A, C
finite_ratioTB, B, D
N is A × B
Den is C × D
round_ratioK, S, N, Den, R
162binarydiv, K, SA, TA, SB, TB, R ←
S is SA xor SB
TA ≡ inf
TB ≡ inf
or
TA = finite0, _
TB = finite0, _
→
TA ≡ inf
→
infinityK, S, R
TB ≡ inf
→
sign_bitsK, S, R
TB = finite0, _
→
infinityK, S, R
else
finite_ratioTA, A, C
finite_ratioTB, B, D
N is A × D
Den is C × B
round_ratioK, S, N, Den, R
ieee_convert+From, +To, +A, -Bdet

B is the format-From float A converted to format To, rounded to nearest; a NaN keeps as much of its payload as fits, aligned at the top, and is quieted.

177ieee_convertFrom, To, A, B ←
decodeFrom, A, S, T
T ≡ inf
→
infinityTo, S, B
T ≡ nan
→
formatFrom, FF, _, _
formatTo, FT, _, _
Payload is A /\ 1 << FF−1
Shift is FT−FF
Shift ≥ 0
→
P is Payload << Shift
else
P is Payload >> −Shift
infinityTo, S, I
U is I \/ P
quietTo, U, B
else
round_ratioTo, S, N, D, B
ieee_compare+K, +Cond, +A, +B, -Rdet

R is 1 when A Cond B holds for format-K floats, else 0. When either operand is a NaN only ne and uo hold, as C's != and SSE's unordered compare give.

192ieee_compareK, Cond, A, B, R ←
decodeK, A, SA, TA
decodeK, B, SB, TB
TA ≡ nan
or
TB ≡ nan
→
Order = unordered
TA ≡ inf
TB ≡ inf
→
compareOrder, SB, SA
TA ≡ inf
→
SA =:= 1
→
Order = (<)
else
Order = (>)
TB ≡ inf
→
SB =:= 1
→
Order = (>)
else
Order = (<)
else
signed_ratioSA, TA, N, D
signed_ratioSB, TB, M, C
X is N × C
Y is M × D
compareOrder, X, Y
conditionCond, Order
→
R = 1
else
R = 0
201conditioneq, =
202conditionne, O ←
O ≢ (=)
203conditionlt, <
204conditionle, O ←
memberchkO, <, =
205conditiongt, >
206conditionge, O ←
memberchkO, >, =
207conditiono, O ←
O ≢ unordered
208conditionuo, unordered
ieee_parse+K, +Cs, -Bsemidet

B is the format-K float nearest the decimal or C hexadecimal literal Cs, without a type suffix. Literals are parsed into a rational before rounding: no double rounding through binary64 for a binary32 literal.

215ieee_parseK, Cs, B ←
phraseliteralS, N, Radix, E, Cs
E ≥ 0
→
A is N × Radix^E
D = 1
else
A = N
D is Radix^−E
round_ratioK, S, A, D, B
221literalS, N, R, E ⟶
222sign1 ⟶
"-"
!
223sign0 ⟶
"+"
!
224sign0 ⟶
[]
225magnitudeN, 2, E ⟶
"0"
X
memberchkX, x, 'X'
!
significand16, N, Places
P
memberchkP, p, 'P'
E is EX−4 × Places
229magnitudeN, 10, E ⟶
significand10, N, Places
E is EX−Places
231decimal_exponentE ⟶
C
memberchkC, e, 'E'
!
232decimal_exponent0 ⟶
[]
233exponent_valueE ⟶
Ds ≢ []
digits_numberDs, 10, 0, N
E is 1−2 × S × N
236significandBase, N, Places ⟶
digit_listBase, Before
fractionBase, After
appendBefore, After, All
All ≢ []
lengthAfter, Places
digits_numberAll, Base, 0, N
240fractionBase, Ds ⟶
"."
!
digit_listBase, Ds
241fraction_, [] ⟶
[]
242digit_listBase, N|Ns ⟶
C
N < Base
!
digit_listBase, Ns
244digit_list_, [] ⟶
[]
245digit_valueC, N ←
char_codeC, V
V ≥ 48
V ≤ 57
→
N is V−48
V ≥ 97
V ≤ 102
→
N is V−87
V ≥ 65
V ≤ 70
→
N is V−55
250digits_number[], _, N, N
251digits_numberD|Ds, Base, A, N ←
A1 is A × Base + D
digits_numberDs, Base, A1, N
ieee_decimal+K, +B, -Cssemidet

Cs is an exact decimal for the finite float B, not necessarily the shortest; fails for an infinity or NaN. Every finite binary float has a terminating decimal expansion. Writing N*5^k e-k avoids padding and is exactly representable by strtof/strtod.

260ieee_decimalK, B, Cs ←
decodeK, B, S, finiteN, E
E ≥ 0
→
M is N << E
P = 0
else
M is N × 5^−E
P = E
number_charsM, MC
number_charsP, PC
S =:= 1
→
Prefix = -
else
Prefix = []
appendPrefix, MC, e, PC, Cs