Float32Exact

1module Float32Exact
2
3import BinaryFloatModel
4import Float32Model
5import FloatLiteralSpec
6import Std.Natural
7
8-- binary32 constants from their exact values.  Each is the binary32 nearest
9-- the value, computed by these functions with unbounded naturals -- the
10-- executable's naturals are 64-bit words and these exact values are not --
11-- so a plan takes a constant as `(compile-time (derivation))`: the compiler
12-- evaluates it and hands the program the word.  pi and ln 2 are bounded to
13-- 40 places, enough that both ends of their interval round alike.
14
15def float32NotRounded : Nat =
16  0x7fc0dead
17
18def float32ExactRational =
19  (lambda unrestricted numerator : Nat .
20    (lambda unrestricted denominator : Nat . (modelRoundBits 0 numerator denominator)))
21
22-- a value known to lie in [lo, hi]: its binary32 when both ends round to it
23def float32ExactBetween =
24  (lambda unrestricted loN : Nat .
25    (lambda unrestricted loD : Nat .
26      (lambda unrestricted hiN : Nat .
27        (lambda unrestricted hiD : Nat .
28          (let unrestricted lo =
29            (float32ExactRational loN loD)
30            in
31            (let unrestricted hi =
32              (float32ExactRational hiN hiD)
33              in
34              (naturalSelect (naturalEqual lo hi) lo float32NotRounded)))))))
35
36-- sqrt(n / d) in [floor(sqrt(n d S^2)) / (d S), that + 1 / (d S)]
37def float32RootScale : Nat =
38  1267650600228229401496703205376
39
40def float32ExactRootLow =
41  (lambda unrestricted n : Nat .
42    (lambda unrestricted d : Nat .
43      (modelIntegerRoot
44        (naturalMultiply (naturalMultiply n d) (naturalMultiply float32RootScale float32RootScale)))))
45
46-- sqrt(n / d): the binary32 nearest it, when the root's interval
47-- [floor(sqrt(n d) S), that + 1] / (d S) rounds alike at both ends
48def float32ExactRoot =
49  (lambda unrestricted n : Nat .
50    (lambda unrestricted d : Nat .
51      (let unrestricted low = (float32ExactRootLow n d) in
52        (float32ExactBetween low (naturalMultiply d float32RootScale)
53          (succ low) (naturalMultiply d float32RootScale)))))
54
55-- A natural below 2^24 as its (exact) binary32 word, with naturals that
56-- stay within a machine word, so a program can convert a size at run time
57-- (float32ExactRational's exact rounding needs unbounded naturals and is a
58-- compile-time conversion). The exponent is floor(log2 n); zero is +0.
59def float32SmallNaturalExponent = (lambda unrestricted n : Nat .
60  (nat-eliminate (lambda unrestricted current : Nat . Nat) 0
61    (lambda unrestricted k : Nat . (lambda unrestricted count : Nat .
62      (naturalAdd count (naturalLessOrEqual (naturalPowerOfTwo (succ k)) n)))) 24))
63def float32ExactSmallNatural = (lambda unrestricted n : Nat .
64  (let unrestricted e = (float32SmallNaturalExponent n) in
65    (naturalSelect (naturalIsZero n) 0
66      (naturalAdd (naturalMultiply (naturalAdd e 127) 8388608)
67        (naturalSaturatingSubtract (naturalMultiply n (naturalPowerOfTwo (naturalSaturatingSubtract 23 e))) 8388608)))))
68
69def float32Places : Nat =
70  10000000000000000000000000000000000000000
71
72def float32Ln2Digits : Nat =
73  6931471805599453094172321214581765680755
74
75def float32PiDigits : Nat =
76  31415926535897932384626433832795028841971
77
78def float32ExactInverseRoot =
79  (lambda unrestricted n : Nat .
80    (let unrestricted low =
81      (float32ExactRootLow 1 n)
82      in
83      (float32ExactBetween
84        low
85        (naturalMultiply n float32RootScale)
86        (succ low)
87        (naturalMultiply n float32RootScale))))
88
89-- The EX2 attention path needs score/sqrt(n) in base-two units. Bound both
90-- sqrt(n) and ln(2) from their rational intervals, then require that the
91-- bounds round to the same F32 word. This avoids rounding 1/sqrt(n) before
92-- multiplying by log2(e), which can choose a different final word.
93def float32ExactInverseRootNaturalLogTwo =
94  (lambda unrestricted n : Nat .
95    (let unrestricted rootLow = (float32ExactRootLow n 1) in
96      (float32ExactBetween
97        (naturalMultiply float32RootScale float32Places)
98        (naturalMultiply (succ rootLow) (succ float32Ln2Digits))
99        (naturalMultiply float32RootScale float32Places)
100        (naturalMultiply rootLow float32Ln2Digits))))
101
102def float32Ln2 : Nat =
103  (compile-time (float32ExactBetween float32Ln2Digits float32Places (succ float32Ln2Digits) float32Places))
104
105def float32Log2E : Nat =
106  (compile-time (float32ExactBetween float32Places (succ float32Ln2Digits) float32Places float32Ln2Digits))
107
108
109-- floor((n / d)^(1 / k) S): the largest r with r^k d <= n S^k, by bisection.
110-- The root is below S (n / d + 1), so its bit length is at most the sum of
111-- theirs -- the bound and the fuel; the target n S^k itself can run past
112-- specBitLength's reach (S = 2^96, k = 32 is 3,073 bits)
113def float32ExactRootKLow =
114  (lambda unrestricted k : Nat .
115    (lambda unrestricted n : Nat .
116      (lambda unrestricted d : Nat .
117        (lambda unrestricted scale : Nat .
118          (let unrestricted target = (nat-multiply n (specPower scale k)) in
119          (let unrestricted bits = (nat-add (specBitLength scale) (specBitLength (nat-add (nat-divide n d) 1))) in
120          (app
121            (nat-eliminate
122              (lambda unrestricted current : Nat . (pi unrestricted low : Nat . (pi unrestricted high : Nat . Nat)))
123              (lambda unrestricted low : Nat . (lambda unrestricted high : Nat . low))
124              (lambda unrestricted predecessor : Nat .
125                (lambda unrestricted induction : (pi unrestricted low : Nat . (pi unrestricted high : Nat . Nat)) .
126                  (lambda unrestricted low : Nat . (lambda unrestricted high : Nat .
127                    (specSelect (nat-less-than (nat-add low 1) high)
128                      (specSelect (nat-less-than target (nat-multiply (specPower (nat-divide (nat-add low high) 2) k) d))
129                        (induction low (nat-divide (nat-add low high) 2))
130                        (induction (nat-divide (nat-add low high) 2) high))
131                      low)))))
132              (nat-add bits 2))
133            zero (specPow2 bits))))))))
134
135-- a value in [lo, hi] / d less the value of `word` (a positive binary32
136-- below 2^24, exponent field at most 150): the binary32 nearest the
137-- difference, its sign set when the word is the larger; zero when the
138-- interval reaches the word's value (the difference is below its width
139def float32ExactRemainder =
140  (lambda unrestricted lo : Nat .
141    (lambda unrestricted hi : Nat .
142      (lambda unrestricted d : Nat .
143        (lambda unrestricted word : Nat .
144          (let unrestricted shift = (specPow2 (nat-subtract 150 (modelExponent word))) in
145          (let unrestricted value = (nat-multiply (nat-add (modelFraction word) modelPow2Twenty3) d) in
146          (let unrestricted low = (nat-multiply lo shift) in
147          (let unrestricted high = (nat-multiply hi shift) in
148          (let unrestricted scale = (nat-multiply d shift) in
149          (specSelect (nat-less-than value low)
150            (float32ExactBetween (nat-subtract low value) scale (nat-subtract high value) scale)
151            (specSelect (nat-less-than high value)
152              (nat-add modelPow2Thirty1 (float32ExactBetween (nat-subtract value high) scale (nat-subtract value low) scale))
153              0)))))))))))
154
155-- 2 pi in two words (the second the nearest to what the first misses), and
156-- 1 / (2 pi): pi in [P, P + 1] / 10^40
157def float32TwoPiHigh : Nat =
158  (compile-time (float32ExactBetween (nat-multiply 2 float32PiDigits) float32Places (nat-multiply 2 (succ float32PiDigits)) float32Places))
159
160def float32TwoPiLow : Nat =
161  (compile-time
162    (float32ExactRemainder (nat-multiply 2 float32PiDigits) (nat-multiply 2 (succ float32PiDigits)) float32Places float32TwoPiHigh))
163
164def float32InverseTwoPi : Nat =
165  (compile-time (float32ExactBetween float32Places (nat-multiply 2 (succ float32PiDigits)) float32Places (nat-multiply 2 float32PiDigits)))