Repository navigation
Expand file tree
/
Copy pathProgram.lean
More file actions
242 lines (233 loc) · 16.8 KB
/
Copy pathProgram.lean
File metadata and controls
242 lines (233 loc) · 16.8 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
import LeanExe.Dialect.Loop
import LeanExe.Dialect.Build
import LeanExe.Dialect.RepeatWhile
import LeanExe.Dialect.InsertErase
import LeanExe.Dialect.Float32
/-! The thirtieth program of the verified compiler: `log`, ported from Arm's optimized-routines,
`math/log.c` and `math/log_data.c` at commit `503fafe311c177de0e571c458c7c337b1ca5f522`
(Copyright (c) 2018-2025, Arm Limited, SPDX-License-Identifier: MIT OR Apache-2.0 WITH
LLVM-exception). The configuration is the one for a target without fused multiply-add:
`LOG_TABLE_BITS` 7, so `N` is 128, `LOG_POLY_ORDER` 6, `LOG_POLY1_ORDER` 12, and the reduction
through `tab2`. The code writes `x = 2 ^ k z` with `z` in `[0x1.6p-1, 0x1.6p0)`, picks the
subinterval `i` of `z` among 128, and computes `log x = k ln 2 + log c + log1p r` with
`r = z/c - 1` from the table entry `c` of the subinterval, `|r| < 1/256`, and a polynomial of
degree 6. Arguments within about `2 ^ -4` of 1 take a polynomial of degree 12 in `x - 1` instead,
with `x - 1` split so that its square loses no bits. Arm's comments bound the error of the first
path, whose results exceed `2 ^ -4` in magnitude, by `0.5 + 2.06/N` units in the last place, plus
the polynomial's error and 0.001 without fused multiply-add, about 0.52 in all, and give about
0.507 units for the second. The module holds the tables in data segments: `logWith` takes them
as parameters, and `log` passes them. Each constant is written as the bit pattern of the
hexadecimal float in Arm's source, and the results of `__math_divzero` and `__math_invalid` from
`math/math_err.c` as the expressions that compute them. `LogAccuracy.lean` checks the results
against a reference computed with Lean's integers. -/
namespace Library.Numerics.Float64.Log
/-- Arm's `__log_data.tab` for `N = 128`: word `2 i` is the bits of `invc`, near `1/c`, and word
`2 i + 1` the bits of `logc`, near `log c`, for the center `c` of subinterval `i`. -/
def logTable : Array UInt64 :=
#[
0x3FF734F0C3E0DE9F, 0xBFD7CC7F79E69000, 0x3FF713786A2CE91F, 0xBFD76FEEC20D0000,
0x3FF6F26008FAB5A0, 0xBFD713E31351E000, 0x3FF6D1A61F138C7D, 0xBFD6B85B38287800,
0x3FF6B1490BC5B4D1, 0xBFD65D5590807800, 0x3FF69147332F0CBA, 0xBFD602D076180000,
0x3FF6719F18224223, 0xBFD5A8CA86909000, 0x3FF6524F99A51ED9, 0xBFD54F4356035000,
0x3FF63356AA8F24C4, 0xBFD4F637C36B4000, 0x3FF614B36B9DDC14, 0xBFD49DA7FDA85000,
0x3FF5F66452C65C4C, 0xBFD445923989A800, 0x3FF5D867B5912C4F, 0xBFD3EDF439B0B800,
0x3FF5BABCCB5B90DE, 0xBFD396CE448F7000, 0x3FF59D61F2D91A78, 0xBFD3401E17BDA000,
0x3FF5805612465687, 0xBFD2E9E2EF468000, 0x3FF56397CEE76BD3, 0xBFD2941B3830E000,
0x3FF54725E2A77F93, 0xBFD23EC58CDA8800, 0x3FF52AFF42064583, 0xBFD1E9E129279000,
0x3FF50F22DBB2BDDF, 0xBFD1956D2B48F800, 0x3FF4F38F4734DED7, 0xBFD141679AB9F800,
0x3FF4D843CFDE2840, 0xBFD0EDD094EF9800, 0x3FF4BD3EC078A3C8, 0xBFD09AA518DB1000,
0x3FF4A27FC3E0258A, 0xBFD047E65263B800, 0x3FF4880524D48434, 0xBFCFEB224586F000,
0x3FF46DCE1B192D0B, 0xBFCF474A7517B000, 0x3FF453D9D3391854, 0xBFCEA4443D103000,
0x3FF43A2744B4845A, 0xBFCE020D44E9B000, 0x3FF420B54115F8FB, 0xBFCD60A22977F000,
0x3FF40782DA3EF4B1, 0xBFCCC00104959000, 0x3FF3EE8F5D57FE8F, 0xBFCC202956891000,
0x3FF3D5D9A00B4CE9, 0xBFCB81178D811000, 0x3FF3BD60C010C12B, 0xBFCAE2C9CCD3D000,
0x3FF3A5242B75DAB8, 0xBFCA45402E129000, 0x3FF38D22CD9FD002, 0xBFC9A877681DF000,
0x3FF3755BC5847A1C, 0xBFC90C6D69483000, 0x3FF35DCE49AD36E2, 0xBFC87120A645C000,
0x3FF34679984DD440, 0xBFC7D68FB4143000, 0x3FF32F5CCEFFCB24, 0xBFC73CB83C627000,
0x3FF3187775A10D49, 0xBFC6A39A9B376000, 0x3FF301C8373E3990, 0xBFC60B3154B7A000,
0x3FF2EB4EBB95F841, 0xBFC5737D76243000, 0x3FF2D50A0219A9D1, 0xBFC4DC7B8FC23000,
0x3FF2BEF9A8B7FD2A, 0xBFC4462C51D20000, 0x3FF2A91C7A0C1BAB, 0xBFC3B08ABC830000,
0x3FF293726014B530, 0xBFC31B996B490000, 0x3FF27DFA5757A1F5, 0xBFC2875490A44000,
0x3FF268B39B1D3BBF, 0xBFC1F3B9F879A000, 0x3FF2539D838FF5BD, 0xBFC160C8252CA000,
0x3FF23EB7AAC9083B, 0xBFC0CE7F57F72000, 0x3FF22A012BA940B6, 0xBFC03CDC49FEA000,
0x3FF2157996CC4132, 0xBFBF57BDBC4B8000, 0x3FF201201DD2FC9B, 0xBFBE370896404000,
0x3FF1ECF4494D480B, 0xBFBD17983EF94000, 0x3FF1D8F5528F6569, 0xBFBBF9674ED8A000,
0x3FF1C52311577E7C, 0xBFBADC79202F6000, 0x3FF1B17C74CB26E9, 0xBFB9C0C3E7288000,
0x3FF19E010C2C1AB6, 0xBFB8A646B372C000, 0x3FF18AB07BB670BD, 0xBFB78D01B3AC0000,
0x3FF1778A25EFBCB6, 0xBFB674F145380000, 0x3FF1648D354C31DA, 0xBFB55E0E6D878000,
0x3FF151B990275FDD, 0xBFB4485CDEA1E000, 0x3FF13F0EA432D24C, 0xBFB333D94D6AA000,
0x3FF12C8B7210F9DA, 0xBFB22079F8C56000, 0x3FF11A3028ECB531, 0xBFB10E4698622000,
0x3FF107FBDA8434AF, 0xBFAFFA6C6AD20000, 0x3FF0F5EE0F4E6BB3, 0xBFADDA8D4A774000,
0x3FF0E4065D2A9FCE, 0xBFABBCECE4850000, 0x3FF0D244632CA521, 0xBFA9A1894012C000,
0x3FF0C0A77CE2981A, 0xBFA788583302C000, 0x3FF0AF2F83C636D1, 0xBFA5715E67D68000,
0x3FF09DDB98A01339, 0xBFA35C8A49658000, 0x3FF08CABAF52E7DF, 0xBFA149E364154000,
0x3FF07B9F2F4E28FB, 0xBF9E72C082EB8000, 0x3FF06AB58C358F19, 0xBF9A55F152528000,
0x3FF059EEA5ECF92C, 0xBF963D62CF818000, 0x3FF04949CDD12C90, 0xBF9228FB8CAA0000,
0x3FF038C6C6F0ADA9, 0xBF8C317B20F90000, 0x3FF02865137932A9, 0xBF8419355DAA0000,
0x3FF0182427EA7348, 0xBF781203C2EC0000, 0x3FF008040614B195, 0xBF60040979240000,
0x3FEFE01FF726FA1A, 0x3F6FEFF384900000, 0x3FEFA11CC261EA74, 0x3F87DC41353D0000,
0x3FEF6310B081992E, 0x3F93CEA3C4C28000, 0x3FEF25F63CEEADCD, 0x3F9B9FC114890000,
0x3FEEE9C8039113E7, 0x3FA1B0D8CE110000, 0x3FEEAE8078CBB1AB, 0x3FA58A5BD001C000,
0x3FEE741AA29D0C9B, 0x3FA95C8340D88000, 0x3FEE3A91830A99B5, 0x3FAD276AEF578000,
0x3FEE01E009609A56, 0x3FB07598E598C000, 0x3FEDCA01E577BB98, 0x3FB253F5E30D2000,
0x3FED92F20B7C9103, 0x3FB42EDD8B380000, 0x3FED5CAC66FB5CCE, 0x3FB606598757C000,
0x3FED272CAA5EDE9D, 0x3FB7DA76356A0000, 0x3FECF26E3E6B2CCD, 0x3FB9AB434E1C6000,
0x3FECBE6DA2A77902, 0x3FBB78C7BB0D6000, 0x3FEC8B266D37086D, 0x3FBD431332E72000,
0x3FEC5894BD5D5804, 0x3FBF0A3171DE6000, 0x3FEC26B533BB9F8C, 0x3FC067152B914000,
0x3FEBF583EEECE73F, 0x3FC147858292B000, 0x3FEBC4FD75DB96C1, 0x3FC2266ECDCA3000,
0x3FEB951E0C864A28, 0x3FC303D7A6C55000, 0x3FEB65E2C5EF3E2C, 0x3FC3DFC33C331000,
0x3FEB374867C9888B, 0x3FC4BA366B7A8000, 0x3FEB094B211D304A, 0x3FC5933928D1F000,
0x3FEADBE885F2EF7E, 0x3FC66ACD2418F000, 0x3FEAAF1D31603DA2, 0x3FC740F8EC669000,
0x3FEA82E63FD358A7, 0x3FC815C0F51AF000, 0x3FEA5740EF09738B, 0x3FC8E92954F68000,
0x3FEA2C2A90AB4B27, 0x3FC9BB3602F84000, 0x3FEA01A01393F2D1, 0x3FCA8BED1C2C0000,
0x3FE9D79F24DB3C1B, 0x3FCB5B515C01D000, 0x3FE9AE2505C7B190, 0x3FCC2967CCBCC000,
0x3FE9852EF297CE2F, 0x3FCCF635D5486000, 0x3FE95CBAEEA44B75, 0x3FCDC1BD3446C000,
0x3FE934C69DE74838, 0x3FCE8C01B8CFE000, 0x3FE90D4F2F6752E6, 0x3FCF5509C0179000,
0x3FE8E6528EFFD79D, 0x3FD00E6C121FB800, 0x3FE8BFCE9FCC007C, 0x3FD071B80E93D000,
0x3FE899C0DABEC30E, 0x3FD0D46B9E867000, 0x3FE87427AA2317FB, 0x3FD13687334BD000,
0x3FE84F00ACB39A08, 0x3FD1980D67234800, 0x3FE82A49E8653E55, 0x3FD1F8FFE0CC8000,
0x3FE8060195F40260, 0x3FD2595FD7636800, 0x3FE7E22563E0A329, 0x3FD2B9300914A800,
0x3FE7BEB377DCB5AD, 0x3FD3187210436000, 0x3FE79BAA679725C2, 0x3FD377266DEC1800,
0x3FE77907F2170657, 0x3FD3D54FFBAF3000, 0x3FE756CADBD6130C, 0x3FD432EEE32FE000]
/-- Arm's `__log_data.tab2` for `N = 128`: word `2 i` is the bits of `chi`, `c` rounded, and word
`2 i + 1` the bits of `clo`, `c - chi` rounded. -/
def logTable2 : Array UInt64 :=
#[
0x3FE61000014FB66B, 0x3C7E026C91425B3C, 0x3FE63000034DB495, 0x3C8DBFEA48005D41,
0x3FE650000D94D478, 0x3C8E7FA786D6A5B7, 0x3FE67000074E6FAD, 0x3C61FCEA6B54254C,
0x3FE68FFFFEDF0FAE, 0xBC7C7E274C590EFD, 0x3FE6B0000763C5BC, 0xBC8AC16848DCDA01,
0x3FE6D0001E5CC1F6, 0x3C833F1C9D499311, 0x3FE6EFFFEB05F63E, 0xBC7E80041AE22D53,
0x3FE710000E869780, 0x3C7BFF6671097952, 0x3FE72FFFFC67E912, 0x3C8C00E226BD8724,
0x3FE74FFFDF81116A, 0xBC6E02916EF101D2, 0x3FE770000F679C90, 0xBC67FC71CD549C74,
0x3FE78FFFFA7EC835, 0x3C81BEC19EF50483, 0x3FE7AFFFFE20C2E6, 0xBC707E1729CC6465,
0x3FE7CFFFED3FC900, 0xBC808072087B8B1C, 0x3FE7EFFFE9261A76, 0x3C8DC0286D9DF9AE,
0x3FE81000049CA3E8, 0x3C897FD251E54C33, 0x3FE8300017932C8F, 0xBC8AFEE9B630F381,
0x3FE850000633739C, 0x3C89BFBF6B6535BC, 0x3FE87000204289C6, 0xBC8BBF65F3117B75,
0x3FE88FFFEBF57904, 0xBC89006EA23DCB57, 0x3FE8B00022BC04DF, 0xBC7D00DF38E04B0A,
0x3FE8CFFFE50C1B8A, 0xBC88007146FF9F05, 0x3FE8EFFFFC918E43, 0x3C83817BD07A7038,
0x3FE910001EFA5FC7, 0x3C893E9176DFB403, 0x3FE9300013467BB9, 0x3C7F804E4B980276,
0x3FE94FFFE6EE076F, 0xBC8F7EF0D9FF622E, 0x3FE96FFFDE3C12D1, 0xBC7082AA962638BA,
0x3FE98FFFF4458A0D, 0xBC87801B9164A8EF, 0x3FE9AFFFDD982E3E, 0xBC8740E08A5A9337,
0x3FE9CFFFED49FB66, 0x3C3FCE08C19BE000, 0x3FE9F00020F19C51, 0xBC8A3FAA27885B0A,
0x3FEA10001145B006, 0x3C74FF489958DA56, 0x3FEA300007BBF6FA, 0x3C8CBEAB8A2B6D18,
0x3FEA500010971D79, 0x3C88FECADD787930, 0x3FEA70001DF52E48, 0xBC8F41763DD8ABDB,
0x3FEA90001C593352, 0xBC8EBF0284C27612, 0x3FEAB0002A4F3E4B, 0xBC69FD043CFF3F5F,
0x3FEACFFFD7AE1ED1, 0xBC823EE7129070B4, 0x3FEAEFFFEE510478, 0x3C6A063EE00EDEA3,
0x3FEB0FFFDB650D5B, 0x3C5A06C8381F0AB9, 0x3FEB2FFFFEAACA57, 0xBC79011E74233C1D,
0x3FEB4FFFD995BADC, 0xBC79FF1068862A9F, 0x3FEB7000249E659C, 0x3C8AFF45D0864F3E,
0x3FEB8FFFF9871640, 0x3C7CFE7796C2C3F9, 0x3FEBAFFFD204CB4F, 0xBC63FF27EEF22BC4,
0x3FEBCFFFD2415C45, 0xBC6CFFB7EE3BEA21, 0x3FEBEFFFF86309DF, 0xBC814103972E0B5C,
0x3FEC0FFFE1B57653, 0x3C8BC16494B76A19, 0x3FEC2FFFF1FA57E3, 0xBC64FEEF8D30C6ED,
0x3FEC4FFFDCBFE424, 0xBC843F68BCEC4775, 0x3FEC6FFFED54B9F7, 0x3C847EA3F053E0EC,
0x3FEC8FFFEB998FD5, 0x3C7383068DF992F1, 0x3FECB0002125219A, 0xBC68FD8E64180E04,
0x3FECCFFFDD94469C, 0x3C8E7EBE1CC7EA72, 0x3FECEFFFEAFDC476, 0x3C8EBE39AD9F88FE,
0x3FED1000169AF82B, 0x3C757D91A8B95A71, 0x3FED30000D0FF71D, 0x3C89C1906970C7DA,
0x3FED4FFFEA790FC4, 0xBC580E37C558FE0C, 0x3FED70002EDC87E5, 0xBC7F80D64DC10F44,
0x3FED900021DC82AA, 0xBC747C8F94FD5C5C, 0x3FEDAFFFD86B0283, 0x3C8C7F1DC521617E,
0x3FEDD000296C4739, 0x3C88019EB2FFB153, 0x3FEDEFFFE54490F5, 0x3C6E00D2C652CC89,
0x3FEE0FFFCDABF694, 0xBC7F8340202D69D2, 0x3FEE2FFFDB52C8DD, 0x3C7B00C1CA1B0864,
0x3FEE4FFFF24216EF, 0x3C72FFA8B094AB51, 0x3FEE6FFFE88A5E11, 0xBC57F673B1EFBE59,
0x3FEE9000119EFF0D, 0xBC84808D5E0BC801, 0x3FEEAFFFDFA51744, 0x3C780006D54320B5,
0x3FEED0001A127FA1, 0xBC5002F860565C92, 0x3FEEF00007BABCC4, 0xBC8540445D35E611,
0x3FEF0FFFF57A8D02, 0xBC4FFB3139EF9105, 0x3FEF30001EE58AC7, 0x3C8A81ACF2731155,
0x3FEF4FFFF5823494, 0x3C8A3F41D4D7C743, 0x3FEF6FFFFCA94C6B, 0xBC6202F41C987875,
0x3FEF8FFFE1F9C441, 0x3C777DD1F477E74B, 0x3FEFAFFFD2E0E37E, 0xBC6F01199A7CA331,
0x3FEFD0001C77E49E, 0x3C7181EE4BCEACB1, 0x3FEFEFFFF7E0C331, 0xBC6E05370170875A,
0x3FF00FFFF465606E, 0xBC8A7EAD491C0ADA, 0x3FF02FFFF3867A58, 0xBC977F69C3FCB2E0,
0x3FF04FFFFDFC0D17, 0x3C97BFFE34CB945B, 0x3FF0700003CD4D82, 0x3C820083C0E456CB,
0x3FF08FFFF9F2CBE8, 0xBC6DFFDFBE37751A, 0x3FF0B000010CDA65, 0xBC913F7FAEE626EB,
0x3FF0D00001A4D338, 0x3C807DFA79489FF7, 0x3FF0EFFFFADAFDFD, 0xBC77040570D66BC0,
0x3FF110000BBAFD96, 0x3C8E80D4846D0B62, 0x3FF12FFFFAE5F45D, 0x3C9DBFFA64FD36EF,
0x3FF150000DD59AD9, 0x3C9A0077701250AE, 0x3FF170000F21559A, 0x3C8DFDF9E2E3DEEE,
0x3FF18FFFFC275426, 0x3C910030DC3B7273, 0x3FF1B000123D3C59, 0x3C997F7980030188,
0x3FF1CFFFF8299EB7, 0xBC65F932AB9F8C67, 0x3FF1EFFFF48AD400, 0x3C937FBF9DA75BEB,
0x3FF210000C8B86A4, 0x3C9F806B91FD5B22, 0x3FF2300003854303, 0x3C93FFC2EB9FBF33,
0x3FF24FFFFFBCF684, 0x3C7601E77E2E2E72, 0x3FF26FFFF52921D9, 0x3C7FFCBB767F0C61,
0x3FF2900014933A3C, 0xBC7202CA3C02412B, 0x3FF2B00014556313, 0xBC92808233F21F02,
0x3FF2CFFFEBFE523B, 0xBC88FF7E384FDCF2, 0x3FF2F0000BB8AD96, 0xBC85FF51503041C5,
0x3FF30FFFFB7AE2AF, 0xBC810071885E289D, 0x3FF32FFFFEAC5F7F, 0xBC91FF5D3FB7B715,
0x3FF350000CA66756, 0x3C957F82228B82BD, 0x3FF3700011FBF721, 0x3C8000BAC40DD5CC,
0x3FF38FFFF9592FB9, 0xBC943F9D2DB2A751, 0x3FF3B00004DDD242, 0x3C857F6B707638E1,
0x3FF3CFFFF5B2C957, 0x3C7A023A10BF1231, 0x3FF3EFFFEAB0B418, 0x3C987F6D66B152B0,
0x3FF410001532AFF4, 0x3C67F8375F198524, 0x3FF4300017478B29, 0x3C8301E672DC5143,
0x3FF44FFFE795B463, 0x3C89FF69B8B2895A, 0x3FF46FFFE80475E0, 0xBC95C0B19BC2F254,
0x3FF48FFFEF6FC1E7, 0x3C9B4009F23A2A72, 0x3FF4AFFFE5BEA704, 0xBC94FFB7BF0D7D45,
0x3FF4D000171027DE, 0xBC99C06471DC6A3D, 0x3FF4F0000FF03EE2, 0x3C977F890B85531C,
0x3FF5100012DC4BD1, 0x3C6004657166A436, 0x3FF530001605277A, 0xBC96BFCECE233209,
0x3FF54FFFECDB704C, 0xBC8902720505A1D7, 0x3FF56FFFEF5F54A9, 0x3C9BBFE60EC96412,
0x3FF5900017E61012, 0x3C887EC581AFEF90, 0x3FF5B00003C93E92, 0xBC9F41080ABF0CC0,
0x3FF5D0001D4919BC, 0xBC98812AFB254729, 0x3FF5EFFFE7B87A89, 0xBC947EB780ED6904]
/-- The path of Arm's `log` for `x` within about `2 ^ -4` of 1, other than 1: the polynomial of
degree 12 in `r = x - 1`, whose term `-r²/2` after `r` is computed from `r = rhi + rlo` with `rhi`
of 26 bits, so that `rhi²` is exact. -/
def logNear1 (x : Float) : Float :=
let B0 := Float.ofBits 0xBFE0000000000000 -- -0x1p-1
let B1 := Float.ofBits 0x3FD5555555555577 -- 0x1.5555555555577p-2
let B2 := Float.ofBits 0xBFCFFFFFFFFFFDCB -- -0x1.ffffffffffdcbp-3
let B3 := Float.ofBits 0x3FC999999995DD0C -- 0x1.999999995dd0cp-3
let B4 := Float.ofBits 0xBFC55555556745A7 -- -0x1.55555556745a7p-3
let B5 := Float.ofBits 0x3FC24924A344DE30 -- 0x1.24924a344de3p-3
let B6 := Float.ofBits 0xBFBFFFFFA4423D65 -- -0x1.fffffa4423d65p-4
let B7 := Float.ofBits 0x3FBC7184282AD6CA -- 0x1.c7184282ad6cap-4
let B8 := Float.ofBits 0xBFB999EB43B068FF -- -0x1.999eb43b068ffp-4
let B9 := Float.ofBits 0x3FB78182F7AFD085 -- 0x1.78182f7afd085p-4
let B10 := Float.ofBits 0xBFB5521375D145CD -- -0x1.5521375d145cdp-4
let r := x - 1.0
let r2 := r * r
let r3 := r * r2
let y := r3 * (B1 + r * B2 + r2 * B3
+ r3 * (B4 + r * B5 + r2 * B6 + r3 * (B7 + r * B8 + r2 * B9 + r3 * B10)))
let w := r * Float.ofBits 0x41A0000000000000 -- 0x1p27
let rhi := r + w - w
let rlo := r - rhi
let w := rhi * rhi * B0
let hi := r + w
let lo := r - hi + w
let lo := lo + B0 * rlo * (rhi + r)
let y := y + lo
y + hi
/-- The main path of Arm's `log` for the bits `ix` of a positive normal `x`, or of a subnormal
`x` scaled by `2 ^ 52` with 52 taken from the exponent field. `tmp >>> 52`, extended by the sign
of `tmp`, is C's arithmetic shift of `(int64_t) tmp`, and `kd` is its conversion to a double. -/
def logCore (tab tab2 : Array UInt64) (ix : UInt64) : Float :=
let ln2hi := Float.ofBits 0x3FE62E42FEFA3800 -- 0x1.62e42fefa3800p-1
let ln2lo := Float.ofBits 0x3D2EF35793C76730 -- 0x1.ef35793c76730p-45
let A0 := Float.ofBits 0xBFE0000000000001 -- -0x1.0000000000001p-1
let A1 := Float.ofBits 0x3FD555555551305B -- 0x1.555555551305bp-2
let A2 := Float.ofBits 0xBFCFFFFFFFEB4590 -- -0x1.fffffffeb459p-3
let A3 := Float.ofBits 0x3FC999B324F10111 -- 0x1.999b324f10111p-3
let A4 := Float.ofBits 0xBFC55575E506C89F -- -0x1.55575e506c89fp-3
let tmp := ix - 0x3fe6000000000000
let i := (tmp >>> 45) % 128
let k := (tmp >>> 52) ||| (if tmp >>> 63 == 1 then 0xFFFFFFFFFFFFF000 else 0)
let iz := ix - (tmp &&& ((0xfff : UInt64) <<< 52))
let invc := Float.ofBits tab[(2 * i).toNat]!
let logc := Float.ofBits tab[(2 * i + 1).toNat]!
let z := Float.ofBits iz
let r := (z - Float.ofBits tab2[(2 * i).toNat]! - Float.ofBits tab2[(2 * i + 1).toNat]!) * invc
let kd := if k >>> 63 == 1 then -(0 - k).toFloat else k.toFloat
let w := kd * ln2hi + logc
let hi := w + r
let lo := w - hi + r + kd * ln2lo
let r2 := r * r
lo + r2 * A0 + r * r2 * (A1 + r * A2 + r2 * (A3 + r * A4)) + hi
/-- Arm's `log` with the tables `tab` and `tab2`. The bits from `1 - 2 ^ -4` to below
`1 + 0x1.09p-4` take `logNear1`, and 1 itself gives 0. `top` is the top 16 bits of `x`: `±0`
gives `-∞` as `__math_divzero` computes it, `∞` gives `∞`, a negative number or NaN gives NaN as
`__math_invalid` computes it, and a subnormal is scaled by `2 ^ 52`. Arm compares 32-bit words,
where `top - 0x0010` wraps for `top < 0x0010`; the comparison on 64-bit words gives the same
result, since the wrapped difference exceeds `0x7ff0 - 0x0010` either way. -/
def logWith (tab tab2 : Array UInt64) (x : Float) : Float :=
let ix := x.toBits
let top := ix >>> 48
if ix - 0x3FEE000000000000 < 0x3FF1090000000000 - 0x3FEE000000000000 then
if ix == 0x3FF0000000000000 then 0.0 else logNear1 x
else if top - 0x0010 ≥ 0x7ff0 - 0x0010 then
if ix * 2 == 0 then -1.0 / 0.0
else if ix == 0x7FF0000000000000 then x
else if top &&& 0x8000 != 0 || top &&& 0x7ff0 == 0x7ff0 then (x - x) / (x - x)
else logCore tab tab2 ((x * Float.ofBits 0x4330000000000000).toBits - ((52 : UInt64) <<< 52))
else logCore tab tab2 ix
/-- `log x`, from the tables `logTable` and `logTable2`. -/
def log (x : Float) : Float := logWith logTable logTable2 x
end Library.Numerics.Float64.Log