pxfm/logs/
log2.rs

1/*
2 * // Copyright (c) Radzivon Bartoshyk 7/2025. All rights reserved.
3 * //
4 * // Redistribution and use in source and binary forms, with or without modification,
5 * // are permitted provided that the following conditions are met:
6 * //
7 * // 1.  Redistributions of source code must retain the above copyright notice, this
8 * // list of conditions and the following disclaimer.
9 * //
10 * // 2.  Redistributions in binary form must reproduce the above copyright notice,
11 * // this list of conditions and the following disclaimer in the documentation
12 * // and/or other materials provided with the distribution.
13 * //
14 * // 3.  Neither the name of the copyright holder nor the names of its
15 * // contributors may be used to endorse or promote products derived from
16 * // this software without specific prior written permission.
17 * //
18 * // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
19 * // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
20 * // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
21 * // DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
22 * // FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
23 * // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
24 * // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
25 * // CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
26 * // OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
27 * // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
28 */
29use crate::common::{f_fmla, min_normal_f64};
30use crate::double_double::DoubleDouble;
31use crate::logs::log2dd::log2_dd;
32use crate::logs::log2td::log2_td;
33use crate::polyeval::f_polyeval6;
34
35pub(crate) static LOG_RANGE_REDUCTION: [u64; 128] = [
36    0x3ff0000000000000,
37    0x3fefc00000000000,
38    0x3fef800000000000,
39    0x3fef400000000000,
40    0x3fef000000000000,
41    0x3feec00000000000,
42    0x3fee800000000000,
43    0x3fee400000000000,
44    0x3fee000000000000,
45    0x3fede00000000000,
46    0x3feda00000000000,
47    0x3fed600000000000,
48    0x3fed400000000000,
49    0x3fed000000000000,
50    0x3fecc00000000000,
51    0x3feca00000000000,
52    0x3fec600000000000,
53    0x3fec400000000000,
54    0x3fec000000000000,
55    0x3febe00000000000,
56    0x3feba00000000000,
57    0x3feb800000000000,
58    0x3feb400000000000,
59    0x3feb200000000000,
60    0x3feae00000000000,
61    0x3feac00000000000,
62    0x3fea800000000000,
63    0x3fea600000000000,
64    0x3fea400000000000,
65    0x3fea000000000000,
66    0x3fe9e00000000000,
67    0x3fe9c00000000000,
68    0x3fe9800000000000,
69    0x3fe9600000000000,
70    0x3fe9400000000000,
71    0x3fe9200000000000,
72    0x3fe9000000000000,
73    0x3fe8c00000000000,
74    0x3fe8a00000000000,
75    0x3fe8800000000000,
76    0x3fe8600000000000,
77    0x3fe8400000000000,
78    0x3fe8000000000000,
79    0x3fe7e00000000000,
80    0x3fe7c00000000000,
81    0x3fe7a00000000000,
82    0x3fe7800000000000,
83    0x3fe7600000000000,
84    0x3fe7400000000000,
85    0x3fe7200000000000,
86    0x3fe7000000000000,
87    0x3fe6e00000000000,
88    0x3fe6c00000000000,
89    0x3fe6a00000000000,
90    0x3fe6800000000000,
91    0x3fe6600000000000,
92    0x3fe6400000000000,
93    0x3fe6200000000000,
94    0x3fe6000000000000,
95    0x3fe5e00000000000,
96    0x3fe5c00000000000,
97    0x3fe5a00000000000,
98    0x3fe5800000000000,
99    0x3fe5600000000000,
100    0x3fe5400000000000,
101    0x3fe5400000000000,
102    0x3fe5200000000000,
103    0x3fe5000000000000,
104    0x3fe4e00000000000,
105    0x3fe4c00000000000,
106    0x3fe4a00000000000,
107    0x3fe4a00000000000,
108    0x3fe4800000000000,
109    0x3fe4600000000000,
110    0x3fe4400000000000,
111    0x3fe4200000000000,
112    0x3fe4000000000000,
113    0x3fe4000000000000,
114    0x3fe3e00000000000,
115    0x3fe3c00000000000,
116    0x3fe3a00000000000,
117    0x3fe3a00000000000,
118    0x3fe3800000000000,
119    0x3fe3600000000000,
120    0x3fe3400000000000,
121    0x3fe3400000000000,
122    0x3fe3200000000000,
123    0x3fe3000000000000,
124    0x3fe3000000000000,
125    0x3fe2e00000000000,
126    0x3fe2c00000000000,
127    0x3fe2c00000000000,
128    0x3fe2a00000000000,
129    0x3fe2800000000000,
130    0x3fe2800000000000,
131    0x3fe2600000000000,
132    0x3fe2400000000000,
133    0x3fe2400000000000,
134    0x3fe2200000000000,
135    0x3fe2000000000000,
136    0x3fe2000000000000,
137    0x3fe1e00000000000,
138    0x3fe1c00000000000,
139    0x3fe1c00000000000,
140    0x3fe1a00000000000,
141    0x3fe1a00000000000,
142    0x3fe1800000000000,
143    0x3fe1600000000000,
144    0x3fe1600000000000,
145    0x3fe1400000000000,
146    0x3fe1400000000000,
147    0x3fe1200000000000,
148    0x3fe1000000000000,
149    0x3fe1000000000000,
150    0x3fe0e00000000000,
151    0x3fe0e00000000000,
152    0x3fe0c00000000000,
153    0x3fe0c00000000000,
154    0x3fe0a00000000000,
155    0x3fe0a00000000000,
156    0x3fe0800000000000,
157    0x3fe0800000000000,
158    0x3fe0600000000000,
159    0x3fe0600000000000,
160    0x3fe0400000000000,
161    0x3fe0400000000000,
162    0x3fe0200000000000,
163    0x3fe0000000000000,
164];
165
166/**
167Generated by SageMath:
168```python
169values = LOG_RANGE_REDUCTION
170
171R = RealField(150)
172
173def hex_to_float(h):
174    return struct.unpack('>d', struct.pack('>Q', h))[0]
175
176real_array = [R(hex_to_float(h)) for h in values]
177
178for r in real_array:
179    print_double_double("", -RealField(180)(r).log2())
180```
181**/
182static LOG2_DD: [(u64, u64); 128] = [
183    (0x0, 0x0),
184    (0x3c26d746128b1857, 0x3f872c7ba20f7327),
185    (0x3c2b2a41b08fbe06, 0x3f9743ee861f3556),
186    (0x3c3491f06c085bc2, 0x3fa184b8e4c56af8),
187    (0x3c477970e03f821c, 0x3fa77394c9d958d5),
188    (0x3c0155660710eb2a, 0x3fad6ebd1f1febfe),
189    (0x3c298c5452bbce74, 0x3fb1bb32a600549d),
190    (0x3c2c141e66faaaad, 0x3fb4c560fe68af88),
191    (0x3c59ced1447e30ad, 0x3fb7d60496cfbb4c),
192    (0xbc321ba488a62577, 0x3fb960caf9abb7ca),
193    (0xbc17936311889913, 0x3fbc7b528b70f1c5),
194    (0xbc36cd4d2ae3a2f6, 0x3fbf9c95dc1d1165),
195    (0x3c55d243efd93259, 0x3fc097e38ce60649),
196    (0xbc5696e2866c718e, 0x3fc22dadc2ab3497),
197    (0xbc67a6ed4e1b0936, 0x3fc3c6fb650cde51),
198    (0xbc61562d61af73f8, 0x3fc494f863b8df35),
199    (0x3c354fae008fbb59, 0x3fc633a8bf437ce1),
200    (0xbc60798d1aa21694, 0x3fc7046031c79f85),
201    (0x3c699aa6df8b7d83, 0x3fc8a8980abfbd32),
202    (0xbc6e95734abd2fcc, 0x3fc97c1cb13c7ec1),
203    (0xbc6cc865b3dd0dbb, 0x3fcb2602497d5346),
204    (0x3c2bc0af7b82e7d7, 0x3fcbfc67a7fff4cc),
205    (0xbc0ba8b1f646ab12, 0x3fcdac22d3e441d3),
206    (0xbc6086fce864a1f6, 0x3fce857d3d361368),
207    (0x3c7768994400ca0a, 0x3fd01d9bbcfa61d4),
208    (0xbc53d56efe4338fe, 0x3fd08bce0d95fa38),
209    (0x3c7c8d43e017579b, 0x3fd169c05363f158),
210    (0x3c6ae9804237ec8e, 0x3fd1d982c9d52708),
211    (0x3c734107c0e54aed, 0x3fd249cd2b13cd6c),
212    (0x3c7968925e378d68, 0x3fd32bfee370ee68),
213    (0x3c7fcad2f4710e00, 0x3fd39de8e1559f6f),
214    (0xbc7c658d602e66b0, 0x3fd4106017c3eca3),
215    (0x3c7e393a16b94b52, 0x3fd4f6fbb2cec598),
216    (0xbc78d86531d55da2, 0x3fd56b22e6b578e5),
217    (0x3c710b5b643a6ecb, 0x3fd5dfdcf1eeae0e),
218    (0x3c7ac9080333c605, 0x3fd6552b49986277),
219    (0x3c2b6d40900b2502, 0x3fd6cb0f6865c8ea),
220    (0x3c68f89e2eb553b2, 0x3fd7b89f02cf2aad),
221    (0x3c651d58525aad39, 0x3fd8304d90c11fd3),
222    (0x3c799aa6df8b7d83, 0x3fd8a8980abfbd32),
223    (0x3c7fdc46af571993, 0x3fd921800924dd3b),
224    (0x3c7bca36fd02def0, 0x3fd99b072a96c6b2),
225    (0x3c5817fd3b7d7e5d, 0x3fda8ff971810a5e),
226    (0xbc45e13b838eba7d, 0x3fdb0b67f4f46810),
227    (0xbc501d98c3531027, 0x3fdb877c57b1b070),
228    (0x3c7edf515c63dd87, 0x3fdc043859e2fdb3),
229    (0x3c6c4aec56233279, 0x3fdc819dc2d45fe4),
230    (0x3c78a38b4175d665, 0x3fdcffae611ad12b),
231    (0xbc6e15a52a31604a, 0x3fdd7e6c0abc3579),
232    (0x3c438c8946414c6a, 0x3fddfdd89d586e2b),
233    (0x3c73bed456b24ed1, 0x3fde7df5fe538ab3),
234    (0x3c76d261f1753e0b, 0x3fdefec61b011f85),
235    (0xbc79ca1a3202b3d7, 0x3fdf804ae8d0cd02),
236    (0xbc87398fe685f171, 0x3fe0014332be0033),
237    (0xbc89c32630008a1f, 0x3fe042bd4b9a7c99),
238    (0xbc7f47806a0e4105, 0x3fe08494c66b8ef0),
239    (0xbc48a33c25e8e226, 0x3fe0c6caaf0c5597),
240    (0xbc73aec658457c41, 0x3fe1096015dee4da),
241    (0xbc8fc7d7c3320aab, 0x3fe14c560fe68af9),
242    (0xbc892ba145dcf40b, 0x3fe18fadb6e2d3c2),
243    (0xbc7f9fb952bbbccc, 0x3fe1d368296b5255),
244    (0xbc8568859c64022e, 0x3fe217868b0c37e8),
245    (0xbc84828ddf1fb145, 0x3fe25c0a0463beb0),
246    (0xbc8c348e4aab18b8, 0x3fe2a0f3c340705c),
247    (0xbc83af881af2f3d9, 0x3fe2e644fac04fd8),
248    (0xbc83af881af2f3d9, 0x3fe2e644fac04fd8),
249    (0x3c8968925e378d68, 0x3fe32bfee370ee68),
250    (0xbc869656a0ad70d4, 0x3fe37222bb70747c),
251    (0x3c76d266d6cdc959, 0x3fe3b8b1c68fa6ed),
252    (0xbc69575b04fa6fbd, 0x3fe3ffad4e74f1d6),
253    (0x3c5b90132aeddb58, 0x3fe44716a2c08262),
254    (0x3c5b90132aeddb58, 0x3fe44716a2c08262),
255    (0xbc75e35482d13dc1, 0x3fe48eef19317991),
256    (0xbc8ca44f1db913d3, 0x3fe4d7380dcc422d),
257    (0x3c7817fd3b7d7e5d, 0x3fe51ff2e30214bc),
258    (0x3c804613e33c06c9, 0x3fe5692101d9b4a6),
259    (0xbc8fc9257edfe9b6, 0x3fe5b2c3da19723b),
260    (0xbc8fc9257edfe9b6, 0x3fe5b2c3da19723b),
261    (0x3c8149a1977b5b99, 0x3fe5fcdce2727ddb),
262    (0xbc8b32266d92c0fe, 0x3fe6476d98ad990a),
263    (0x3c821d90b84e7218, 0x3fe6927781d932a8),
264    (0x3c821d90b84e7218, 0x3fe6927781d932a8),
265    (0x3c7f6e91ad16ecff, 0x3fe6ddfc2a78fc63),
266    (0xbc74a31ce1b7e328, 0x3fe729fd26b707c8),
267    (0x3c6a7b47d2c352d9, 0x3fe7767c12967a45),
268    (0x3c6a7b47d2c352d9, 0x3fe7767c12967a45),
269    (0x3c821f9cb2cc5575, 0x3fe7c37a9227e7fb),
270    (0x3c8dc572667587b1, 0x3fe810fa51bf65fd),
271    (0x3c8dc572667587b1, 0x3fe810fa51bf65fd),
272    (0xbc78f93e7aa3bdf8, 0x3fe85efd062c656d),
273    (0x3c5b85a54d7ee2fd, 0x3fe8ad846cf369a4),
274    (0x3c5b85a54d7ee2fd, 0x3fe8ad846cf369a4),
275    (0x3c8bf1d926766301, 0x3fe8fc924c89ac84),
276    (0x3c401ee1343fe7ca, 0x3fe94c287492c4db),
277    (0x3c401ee1343fe7ca, 0x3fe94c287492c4db),
278    (0x3c7fa0a62e6add1b, 0x3fe99c48be2063c8),
279    (0x3c8022ddb71189c5, 0x3fe9ecf50bf43f13),
280    (0x3c8022ddb71189c5, 0x3fe9ecf50bf43f13),
281    (0x3c7ac7fc60a51031, 0x3fea3e2f4ac43f60),
282    (0x3c6817fd3b7d7e5d, 0x3fea8ff971810a5e),
283    (0x3c6817fd3b7d7e5d, 0x3fea8ff971810a5e),
284    (0xbc83138e941643f7, 0x3feae255819f022d),
285    (0x3c8e0ae0d3f8a58b, 0x3feb35458761d479),
286    (0x3c8e0ae0d3f8a58b, 0x3feb35458761d479),
287    (0x3c742b7579f0f8d4, 0x3feb88cb9a2ab521),
288    (0x3c742b7579f0f8d4, 0x3feb88cb9a2ab521),
289    (0x3c6a7610e40bd6ab, 0x3febdce9dcc96187),
290    (0xbc80e5edaecee150, 0x3fec31a27dd00b4a),
291    (0xbc80e5edaecee150, 0x3fec31a27dd00b4a),
292    (0xbc831d962d3728cc, 0x3fec86f7b7ea4a89),
293    (0xbc831d962d3728cc, 0x3fec86f7b7ea4a89),
294    (0xbc857391924a6d9d, 0x3fecdcebd2373995),
295    (0x3c78333ac7d9ebbb, 0x3fed338120a6dd9d),
296    (0x3c78333ac7d9ebbb, 0x3fed338120a6dd9d),
297    (0xbc86c0268890da53, 0x3fed8aba045b01c8),
298    (0xbc86c0268890da53, 0x3fed8aba045b01c8),
299    (0xbc859e7ba5d5ccc9, 0x3fede298ec0bac0d),
300    (0xbc859e7ba5d5ccc9, 0x3fede298ec0bac0d),
301    (0x3c60b07079619c47, 0x3fee3b20546f554a),
302    (0x3c60b07079619c47, 0x3fee3b20546f554a),
303    (0xbc8cc4d81bc25adf, 0x3fee9452c8a71028),
304    (0xbc8cc4d81bc25adf, 0x3fee9452c8a71028),
305    (0xbc776c0a2827d49a, 0x3feeee32e2aeccbf),
306    (0xbc776c0a2827d49a, 0x3feeee32e2aeccbf),
307    (0xbc8314dc4fc42302, 0x3fef48c34bd1e96f),
308    (0xbc8314dc4fc42302, 0x3fef48c34bd1e96f),
309    (0xbc817f8e37b00179, 0x3fefa406bd2443df),
310    (0x0, 0x0),
311];
312
313#[cfg(not(any(
314    all(
315        any(target_arch = "x86", target_arch = "x86_64"),
316        target_feature = "fma"
317    ),
318    all(target_arch = "aarch64", target_feature = "neon")
319)))]
320pub(crate) static LOG_CD: [u64; 128] = [
321    0x0000000000000000,
322    0xbf10000000000000,
323    0xbf30000000000000,
324    0xbf42000000000000,
325    0xbf50000000000000,
326    0xbf59000000000000,
327    0xbf62000000000000,
328    0xbf68800000000000,
329    0xbf70000000000000,
330    0xbf49000000000000,
331    0xbf5f000000000000,
332    0xbf69c00000000000,
333    0xbf30000000000000,
334    0xbf5c000000000000,
335    0xbf6b000000000000,
336    0xbf45000000000000,
337    0xbf64000000000000,
338    0x3f10000000000000,
339    0xbf60000000000000,
340    0x3f3a000000000000,
341    0xbf5e000000000000,
342    0x3f38000000000000,
343    0xbf61000000000000,
344    0xbf00000000000000,
345    0xbf66000000000000,
346    0xbf4a000000000000,
347    0xbf6e000000000000,
348    0xbf5f800000000000,
349    0xbf30000000000000,
350    0xbf6c000000000000,
351    0xbf5f000000000000,
352    0xbf3c000000000000,
353    0xbf70000000000000,
354    0xbf65400000000000,
355    0xbf56000000000000,
356    0xbf24000000000000,
357    0x3f50000000000000,
358    0xbf68800000000000,
359    0xbf60800000000000,
360    0xbf52000000000000,
361    0xbf30000000000000,
362    0x3f42000000000000,
363    0xbf70000000000000,
364    0xbf6ac00000000000,
365    0xbf66000000000000,
366    0xbf61c00000000000,
367    0xbf5c000000000000,
368    0xbf55800000000000,
369    0xbf50000000000000,
370    0xbf47000000000000,
371    0xbf40000000000000,
372    0xbf36000000000000,
373    0xbf30000000000000,
374    0xbf2c000000000000,
375    0xbf30000000000000,
376    0xbf36000000000000,
377    0xbf40000000000000,
378    0xbf47000000000000,
379    0xbf50000000000000,
380    0xbf55800000000000,
381    0xbf5c000000000000,
382    0xbf61c00000000000,
383    0xbf66000000000000,
384    0xbf6ac00000000000,
385    0xbf70000000000000,
386    0x3f55000000000000,
387    0x3f42000000000000,
388    0xbf30000000000000,
389    0xbf52000000000000,
390    0xbf60800000000000,
391    0xbf68800000000000,
392    0x3f60c00000000000,
393    0x3f50000000000000,
394    0xbf24000000000000,
395    0xbf56000000000000,
396    0xbf65400000000000,
397    0xbf70000000000000,
398    0x3f50000000000000,
399    0xbf3c000000000000,
400    0xbf5f000000000000,
401    0xbf6c000000000000,
402    0x3f56800000000000,
403    0xbf30000000000000,
404    0xbf5f800000000000,
405    0xbf6e000000000000,
406    0x3f51000000000000,
407    0xbf4a000000000000,
408    0xbf66000000000000,
409    0x3f60000000000000,
410    0xbf00000000000000,
411    0xbf61000000000000,
412    0x3f64800000000000,
413    0x3f38000000000000,
414    0xbf5e000000000000,
415    0x3f66000000000000,
416    0x3f3a000000000000,
417    0xbf60000000000000,
418    0x3f64800000000000,
419    0x3f10000000000000,
420    0xbf64000000000000,
421    0x3f60000000000000,
422    0xbf45000000000000,
423    0xbf6b000000000000,
424    0x3f51000000000000,
425    0xbf5c000000000000,
426    0x3f65400000000000,
427    0xbf30000000000000,
428    0xbf69c00000000000,
429    0x3f52000000000000,
430    0xbf5f000000000000,
431    0x3f63000000000000,
432    0xbf49000000000000,
433    0xbf70000000000000,
434    0x3f30000000000000,
435    0xbf68800000000000,
436    0x3f52800000000000,
437    0xbf62000000000000,
438    0x3f5f000000000000,
439    0xbf59000000000000,
440    0x3f64c00000000000,
441    0xbf50000000000000,
442    0x3f69000000000000,
443    0xbf42000000000000,
444    0x3f6c400000000000,
445    0xbf30000000000000,
446    0x3f6e800000000000,
447    0xbf10000000000000,
448    0xbf70000000000000,
449];
450
451pub(crate) const LOG_COEFFS: [u64; 6] = [
452    0xbfdfffffffffffff,
453    0x3fd5555555554a9b,
454    0xbfd0000000094567,
455    0x3fc99999dcc9823c,
456    0xbfc55550ac2e537a,
457    0x3fc21a02c4e624d7,
458];
459
460/// Log2(x)
461///
462/// Max found ULP 0.5
463pub fn f_log2(x: f64) -> f64 {
464    let mut x_u = x.to_bits();
465
466    const E_BIAS: u64 = (1u64 << (11 - 1u64)) - 1u64;
467
468    let mut x_e: i64 = -(E_BIAS as i64);
469
470    const MAX_NORMAL: u64 = f64::to_bits(f64::MAX);
471
472    if x_u == 1f64.to_bits() {
473        // log2(1.0) = +0.0
474        return 0.0;
475    }
476    if x_u < min_normal_f64().to_bits() || x_u > MAX_NORMAL {
477        if x == 0.0 {
478            return f64::NEG_INFINITY;
479        }
480        if x < 0. || x.is_nan() {
481            return f64::NAN;
482        }
483        if x.is_infinite() || x.is_nan() {
484            return x + x;
485        }
486        // Normalize denormal inputs.
487        x_u = (x * f64::from_bits(0x4330000000000000)).to_bits();
488        x_e -= 52;
489    }
490
491    // log2(x) = log2(2^x_e * x_m)
492    //         = x_e + log2(x_m)
493    // Range reduction for log2(x_m):
494    // For each x_m, we would like to find r such that:
495    //   -2^-8 <= r * x_m - 1 < 2^-7
496    let shifted = (x_u >> 45) as i64;
497    let index = shifted & 0x7F;
498
499    // Add unbiased exponent. Add an extra 1 if the 8 leading fractional bits are
500    // all 1's.
501    x_e = x_e.wrapping_add(x_u.wrapping_add(1u64 << 45).wrapping_shr(52) as i64);
502    let e_x = x_e as f64;
503
504    // Set m = 1.mantissa.
505    let x_m = (x_u & 0x000F_FFFF_FFFF_FFFFu64) | 0x3FF0_0000_0000_0000u64;
506    let m = f64::from_bits(x_m);
507
508    let u;
509
510    let r = f64::from_bits(LOG_RANGE_REDUCTION[index as usize]);
511
512    #[cfg(any(
513        all(
514            any(target_arch = "x86", target_arch = "x86_64"),
515            target_feature = "fma"
516        ),
517        all(target_arch = "aarch64", target_feature = "neon")
518    ))]
519    {
520        u = f_fmla(r, m, -1.0); // exact   
521    }
522    #[cfg(not(any(
523        all(
524            any(target_arch = "x86", target_arch = "x86_64"),
525            target_feature = "fma"
526        ),
527        all(target_arch = "aarch64", target_feature = "neon")
528    )))]
529    {
530        let c_m = x_m & 0x3FFF_E000_0000_0000u64;
531        let c = f64::from_bits(c_m);
532        u = f_fmla(r, m - c, f64::from_bits(LOG_CD[index as usize])); // exact
533    }
534
535    // Exact sum:
536    let log_vals = DoubleDouble::from_bit_pair(LOG2_DD[index as usize]);
537
538    /*
539       Poly generated in Sollya;
540       d = [-2^-8, 2^-7];
541       f = (log2(1 + x))/x;
542       pf = fpminimax(f, 5, [|D...|], [-2^-8, 2^-7]);
543
544       See ./notes/log2.sollya
545    */
546    const C: [u64; 6] = [
547        0x3ff71547652b82fe,
548        0xbfe71547652b7a07,
549        0x3fdec709dc458db1,
550        0xbfd715479c2266c9,
551        0x3fd2776ae1ddf8f0,
552        0xbfce7b2178870157,
553    ];
554
555    // Degree-7 minimax polynomial
556    let p = f_fmla(
557        f_polyeval6(
558            u,
559            f64::from_bits(C[0]),
560            f64::from_bits(C[1]),
561            f64::from_bits(C[2]),
562            f64::from_bits(C[3]),
563            f64::from_bits(C[4]),
564            f64::from_bits(C[5]),
565        ),
566        u,
567        log_vals.lo,
568    );
569
570    // log2(x) = e + log2(r) + log2(index)
571
572    let rr = DoubleDouble::from_full_exact_add(log_vals.hi, e_x);
573    let r3 = DoubleDouble::add_f64(rr, p);
574
575    let err = f_fmla(
576        r3.hi,
577        f64::from_bits(0x3bf0000000000000), // 2^-64
578        f64::from_bits(0x3c60000000000000), // 2^-57
579    );
580
581    // Lower bound from the result
582    let left = r3.hi + (r3.lo - err);
583    // Upper bound from the result
584    let right = r3.hi + (r3.lo + err);
585
586    // Ziv's test if fast pass is accurate enough.
587    if left == right {
588        return left;
589    }
590
591    log2_hard(x)
592}
593
594#[cold]
595#[inline(never)]
596fn log2_hard(x: f64) -> f64 {
597    let r = log2_dd(x);
598    let err = f_fmla(
599        r.hi,
600        f64::from_bits(0x3b50000000000000), // 2^-74
601        f64::from_bits(0x3990000000000000), // 2^-102
602    );
603    let ub = r.hi + (r.lo + err);
604    let lb = r.hi + (r.lo - err);
605    if ub == lb {
606        return r.to_f64();
607    }
608    log2_hard_slow(x)
609}
610
611#[cold]
612#[inline(never)]
613fn log2_hard_slow(x: f64) -> f64 {
614    log2_td(x).to_f64()
615}
616
617#[cfg(test)]
618mod tests {
619    use super::*;
620
621    #[test]
622    fn test_log2d() {
623        assert_eq!(f_log2(0.7500000000000461), -0.4150374992787552);
624        assert_eq!(f_log2(1.99999999779061), 0.999999998406262);
625        assert_eq!(f_log2(0.7812499981882864), -0.35614381357087554);
626        assert_eq!(f_log2(0.5937499997089616), -0.7520724872635803);
627        assert_eq!(f_log2(0.9983015060407214), -0.0024524921736992287);
628        assert_eq!(f_log2(0.2284926469437778), -2.129780355811612);
629        assert_eq!(f_log2(24.), 4.584962500721156181453738943);
630        assert!((f_log2(0.35) - 0.35f64.log2()).abs() < 1e-8);
631        assert!((f_log2(0.9) - 0.9f64.log2()).abs() < 1e-8);
632        assert_eq!(f_log2(0.), f64::NEG_INFINITY);
633        assert!(f_log2(-1.).is_nan());
634        assert!(f_log2(f64::NAN).is_nan());
635        assert!(f_log2(f64::NEG_INFINITY).is_nan());
636        assert_eq!(f_log2(f64::INFINITY), f64::INFINITY);
637    }
638
639    #[test]
640    fn test_log2_control_values() {
641        assert_eq!(
642            f_log2(f64::from_bits(0x3ff1406d79e1b574)),
643            0.10866415915666149
644        );
645        assert_eq!(
646            f_log2(f64::from_bits(0x3ffb4ebe40c95a01)),
647            0.7712301192941611
648        );
649        assert_eq!(
650            f_log2(f64::from_bits(0x3ff0b541b6746bd1)),
651            0.062470074690080965
652        );
653        assert_eq!(
654            f_log2(f64::from_bits(0x3ff68d778f076021)),
655            0.4952222169409832
656        );
657        assert_eq!(
658            f_log2(f64::from_bits(0x3ff67eaf07ce24d1)),
659            0.4915233709893074
660        );
661        assert_eq!(
662            f_log2(f64::from_bits(0x3ff0b53197bd66c8)),
663            0.062448835548542664
664        );
665        assert_eq!(
666            f_log2(f64::from_bits(0x3fe6b015f8d9a784)),
667            -0.496152942566177
668        );
669        assert_eq!(
670            f_log2(f64::from_bits(0x3ff160732376ad7f)),
671            0.11908694357502585
672        );
673        assert_eq!(
674            f_log2(f64::from_bits(0x3ff0b53f741fb8c4)),
675            0.062467098188572705
676        );
677    }
678}
pxfm/logs/log2.rs

pxfm/logs/
log2.rs
