Float32Exact

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