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
46def float32Places : Nat =
47  10000000000000000000000000000000000000000
48
49def float32Ln2Digits : Nat =
50  6931471805599453094172321214581765680755
51
52def float32PiDigits : Nat =
53  31415926535897932384626433832795028841971
54
55def float32ExactInverseRoot =
56  (lambda unrestricted n : Nat .
57    (let unrestricted low =
58      (float32ExactRootLow 1 n)
59      in
60      (float32ExactBetween
61        low
62        (naturalMultiply n float32RootScale)
63        (succ low)
64        (naturalMultiply n float32RootScale))))
65
66-- The EX2 attention path needs score/sqrt(n) in base-two units. Bound both
67-- sqrt(n) and ln(2) from their rational intervals, then require that the
68-- bounds round to the same F32 word. This avoids rounding 1/sqrt(n) before
69-- multiplying by log2(e), which can choose a different final word.
70def float32ExactInverseRootNaturalLogTwo =
71  (lambda unrestricted n : Nat .
72    (let unrestricted rootLow = (float32ExactRootLow n 1) in
73      (float32ExactBetween
74        (naturalMultiply float32RootScale float32Places)
75        (naturalMultiply (succ rootLow) (succ float32Ln2Digits))
76        (naturalMultiply float32RootScale float32Places)
77        (naturalMultiply rootLow float32Ln2Digits))))
78
79def float32Ln2 : Nat =
80  (compile-time (float32ExactBetween float32Ln2Digits float32Places (succ float32Ln2Digits) float32Places))
81
82def float32Log2E : Nat =
83  (compile-time (float32ExactBetween float32Places (succ float32Ln2Digits) float32Places float32Ln2Digits))
84
85
86-- floor((n / d)^(1 / k) S): the largest r with r^k d <= n S^k, by bisection.
87-- The root is below S (n / d + 1), so its bit length is at most the sum of
88-- theirs -- the bound and the fuel; the target n S^k itself can run past
89-- specBitLength's reach (S = 2^96, k = 32 is 3,073 bits)
90def float32ExactRootKLow =
91  (lambda unrestricted k : Nat .
92    (lambda unrestricted n : Nat .
93      (lambda unrestricted d : Nat .
94        (lambda unrestricted scale : Nat .
95          (let unrestricted target = (nat-multiply n (specPower scale k)) in
96          (let unrestricted bits = (nat-add (specBitLength scale) (specBitLength (nat-add (nat-divide n d) 1))) in
97          (app
98            (nat-eliminate
99              (lambda unrestricted current : Nat . (pi unrestricted low : Nat . (pi unrestricted high : Nat . Nat)))
100              (lambda unrestricted low : Nat . (lambda unrestricted high : Nat . low))
101              (lambda unrestricted predecessor : Nat .
102                (lambda unrestricted induction : (pi unrestricted low : Nat . (pi unrestricted high : Nat . Nat)) .
103                  (lambda unrestricted low : Nat . (lambda unrestricted high : Nat .
104                    (specSelect (nat-less-than (nat-add low 1) high)
105                      (specSelect (nat-less-than target (nat-multiply (specPower (nat-divide (nat-add low high) 2) k) d))
106                        (induction low (nat-divide (nat-add low high) 2))
107                        (induction (nat-divide (nat-add low high) 2) high))
108                      low)))))
109              (nat-add bits 2))
110            zero (specPow2 bits))))))))
111
112-- a value in [lo, hi] / d less the value of `word` (a positive binary32
113-- below 2^24, exponent field at most 150): the binary32 nearest the
114-- difference, its sign set when the word is the larger; zero when the
115-- interval reaches the word's value (the difference is below its width
116def float32ExactRemainder =
117  (lambda unrestricted lo : Nat .
118    (lambda unrestricted hi : Nat .
119      (lambda unrestricted d : Nat .
120        (lambda unrestricted word : Nat .
121          (let unrestricted shift = (specPow2 (nat-subtract 150 (modelExponent word))) in
122          (let unrestricted value = (nat-multiply (nat-add (modelFraction word) modelPow2Twenty3) d) in
123          (let unrestricted low = (nat-multiply lo shift) in
124          (let unrestricted high = (nat-multiply hi shift) in
125          (let unrestricted scale = (nat-multiply d shift) in
126          (specSelect (nat-less-than value low)
127            (float32ExactBetween (nat-subtract low value) scale (nat-subtract high value) scale)
128            (specSelect (nat-less-than high value)
129              (nat-add modelPow2Thirty1 (float32ExactBetween (nat-subtract value high) scale (nat-subtract value low) scale))
130              0)))))))))))
131
132-- 2 pi in two words (the second the nearest to what the first misses), and
133-- 1 / (2 pi): pi in [P, P + 1] / 10^40
134def float32TwoPiHigh : Nat =
135  (compile-time (float32ExactBetween (nat-multiply 2 float32PiDigits) float32Places (nat-multiply 2 (succ float32PiDigits)) float32Places))
136
137def float32TwoPiLow : Nat =
138  (compile-time
139    (float32ExactRemainder (nat-multiply 2 float32PiDigits) (nat-multiply 2 (succ float32PiDigits)) float32Places float32TwoPiHigh))
140
141def float32InverseTwoPi : Nat =
142  (compile-time (float32ExactBetween float32Places (nat-multiply 2 (succ float32PiDigits)) float32Places (nat-multiply 2 float32PiDigits)))