storage_element_format.bend
stelliferous@0.0.3.0 · 145 lines · 0 definitions
1# Stored element formats and the FP32 value of each. An element arrives as its
2# bytes, least significant first, gathered into a low word and a high word.
3import Base
4import stelliferous@0.0.3.0/math_fp32.bend as Fp32
5
6type Format is Data:
7 Fp32{}
8 Fp64{}
9 Fp16{}
10 BFloat16{}
11 Int64{}
12 Int32{}
13
14# Bytes per element.
15def width(format: Format) -> U32:
16 match format:
17 case Fp32{}: 4
18 case Fp64{}: 8
19 case Fp16{}: 2
20 case BFloat16{}: 2
21 case Int64{}: 8
22 case Int32{}: 4
23
24# IEEE half: a sign, five exponent bits and ten fraction bits. A normal value
25# keeps its fraction and gains 112 on the exponent; infinities and NaNs keep
26# their payload; a subnormal is normalized by doubling the fraction until its
27# leading bit is at position ten, so widening is exact.
28
29def normalized(fuel: Nat,done: Bool,+sign: U32,+exponent: U32,+mantissa: U32) -> F32:
30 match fuel done:
31 case fuel True{}: Fp32.from_bits((sign + exponent * 8388608 + (mantissa - 1024) * 8192 : U32))
32 case 0n False{}: Fp32.from_bits(sign)
33 case 1n+rest False{}:
34 +doubled = (mantissa * 2 : U32)
35 normalized(rest,U32.is_le(1024,doubled),sign,(exponent - 1 : U32),doubled)
36
37def subnormal(zero: Bool,+sign: U32,+mantissa: U32) -> F32:
38 match zero:
39 case True{}: Fp32.from_bits(sign)
40 case False{}: normalized(11n,False{},sign,113,mantissa)
41
42def half_classified(lowest: Bool,highest: Bool,+sign: U32,+exponent: U32,+mantissa: U32) -> F32:
43 match lowest highest:
44 case True{} highest: subnormal(U32.is_zero(mantissa),sign,mantissa)
45 case False{} True{}: Fp32.from_bits((sign + 2139095040 + mantissa * 8192 : U32))
46 case False{} False{}: Fp32.from_bits((sign + (exponent + 112) * 8388608 + mantissa * 8192 : U32))
47
48def half(+word: U32) -> F32:
49 +exponent = (word / 1024 % 32 : U32)
50 half_classified(U32.is_zero(exponent),U32.is_eq(exponent,31),(word / 32768 % 2 * 2147483648 : U32),exponent,(word % 1024 : U32))
51
52# A two's complement integer rounds to the nearest FP32, ties to even, as a C
53# conversion does. The magnitude is scaled until its leading bit is at position
54# 23 of the low word: a small one is doubled, which is exact, and a large one is
55# halved while the dropped bits are kept as a guard bit and a sticky bit.
56
57def carried(full: Bool,+sign: U32,+exponent: U32,+low: U32) -> F32:
58 match full:
59 case True{}: Fp32.from_bits((sign + (exponent + 1) * 8388608 : U32))
60 case False{}: Fp32.from_bits((sign + exponent * 8388608 + (low - 8388608) : U32))
61
62def rounded(up: Bool,+sign: U32,+exponent: U32,+low: U32) -> F32:
63 match up:
64 case False{}: Fp32.from_bits((sign + exponent * 8388608 + (low - 8388608) : U32))
65 case True{}: carried(U32.is_eq((low + 1 : U32),16777216),sign,exponent,(low + 1 : U32))
66
67def halved(fuel: Nat,done: Bool,+sign: U32,+exponent: U32,+low: U32,+high: U32,+guard: U32,+sticky: U32) -> F32:
68 match fuel done:
69 case fuel True{}: rounded(U32.is_eq(guard,1) && (U32.is_eq(sticky,1) || U32.is_eq((low % 2 : U32),1)),sign,exponent,low)
70 case 0n False{}: Fp32.from_bits(sign)
71 case 1n+rest False{}:
72 +next_low = (low / 2 + high % 2 * 2147483648 : U32)
73 +next_high = (high / 2 : U32)
74 halved(rest,U32.is_zero(next_high) && U32.is_lt(next_low,16777216),sign,(exponent + 1 : U32),next_low,next_high,(low % 2 : U32),U32.max(sticky,guard))
75
76def doubled(fuel: Nat,done: Bool,+sign: U32,+exponent: U32,+low: U32,+high: U32) -> F32:
77 match fuel done:
78 case fuel True{}: halved(41n,U32.is_zero(high) && U32.is_lt(low,16777216),sign,exponent,low,high,0,0)
79 case 0n False{}: Fp32.from_bits(sign)
80 case 1n+rest False{}:
81 +twice = (low * 2 : U32)
82 doubled(rest,U32.is_le(8388608,twice),sign,(exponent - 1 : U32),twice,high)
83
84def scaled(zero: Bool,+sign: U32,+low: U32,+high: U32) -> F32:
85 match zero:
86 case True{}: Fp32.from_bits(0)
87 case False{}: doubled(24n,U32.is_le(8388608,low) || Bool.not(U32.is_zero(high)),sign,150,low,high)
88
89def magnitude(negative: Bool,+low: U32,+high: U32) -> F32:
90 match negative:
91 case False{}: scaled(U32.is_zero(low) && U32.is_zero(high),0,low,high)
92 case True{}: scaled(False{},2147483648,(0 - low : U32),(0 - high - Bool.pick(U32,U32.is_zero(low),0,1) : U32))
93
94def integer(+low: U32,+high: U32) -> F32:
95 magnitude(U32.is_le(2147483648,high),low,high)
96
97def extended(+low: U32) -> U32:
98 Bool.pick(U32,U32.is_le(2147483648,low),4294967295,0)
99
100# IEEE double: a sign, eleven exponent bits and 52 fraction bits, the top 20 of
101# them in the high word. It rounds to the nearest FP32, ties to even, as a C
102# conversion does. With its leading bit restored, the fraction is a 53-bit
103# integer significand: a normal FP32 is that integer halved as above, starting
104# from the exponent that keeps the value; an FP32 subnormal is the significand
105# shifted to the subnormal scale with the same guard and sticky bits, which
106# carries into the smallest normal when it rounds up; anything smaller is a
107# signed zero. An infinity stays infinite, and so does a value above the FP32
108# range. A NaN stays a quiet NaN that keeps the top of its payload.
109
110def shifted(fuel: Nat,done: Bool,+sign: U32,+remaining: U32,+low: U32,+high: U32,+guard: U32,+sticky: U32) -> F32:
111 match fuel done:
112 case fuel True{}:
113 Fp32.from_bits((sign + low + Bool.pick(U32,U32.is_eq(guard,1) && (U32.is_eq(sticky,1) || U32.is_eq((low % 2 : U32),1)),1,0) : U32))
114 case 0n False{}: Fp32.from_bits(sign)
115 case 1n+rest False{}:
116 shifted(rest,U32.is_eq(remaining,1),sign,(remaining - 1 : U32),(low / 2 + high % 2 * 2147483648 : U32),(high / 2 : U32),(low % 2 : U32),U32.max(sticky,guard))
117
118def special(infinite: Bool,+sign: U32,+low: U32,+fraction: U32) -> F32:
119 match infinite:
120 case True{}: Fp32.from_bits((sign + 2139095040 : U32))
121 case False{}: Fp32.from_bits((sign + 2143289344 + (fraction * 8 + low / 536870912) % 4194304 : U32))
122
123def double_classified(nonfinite: Bool,large: Bool,normal: Bool,small: Bool,+sign: U32,+exponent: U32,+low: U32,+fraction: U32) -> F32:
124 match nonfinite large normal small:
125 case True{} large normal small: special(U32.is_zero(fraction) && U32.is_zero(low),sign,low,fraction)
126 case False{} True{} normal small: Fp32.from_bits((sign + 2139095040 : U32))
127 case False{} False{} True{} small: halved(41n,False{},sign,(exponent - 925 : U32),low,(fraction + 1048576 : U32),0,0)
128 case False{} False{} False{} True{}: shifted(55n,False{},sign,(926 - exponent : U32),low,(fraction + 1048576 : U32),0,0)
129 case False{} False{} False{} False{}: Fp32.from_bits(sign)
130
131def double(+low: U32,+high: U32) -> F32:
132 +exponent = (high / 1048576 % 2048 : U32)
133 double_classified(U32.is_eq(exponent,2047),U32.is_le(1151,exponent),U32.is_le(897,exponent),U32.is_le(872,exponent),
134 (high / 2147483648 * 2147483648 : U32),exponent,low,(high % 1048576 : U32))
135
136# The value of the element whose first four bytes form `low` and whose later
137# bytes form `high`.
138def widen(format: Format,+low: U32,+high: U32) -> F32:
139 match format:
140 case Fp32{}: Fp32.from_bits(low)
141 case Fp64{}: double(low,high)
142 case Fp16{}: half(low)
143 case BFloat16{}: Fp32.from_bits((low * 65536 : U32))
144 case Int64{}: integer(low,high)
145 case Int32{}: integer(low,extended(low))