pxfm/logs/log10.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::log2::LOG_COEFFS;
32use crate::logs::log10dd::log10_dd;
33use crate::logs::log10td::log10_td;
34use crate::polyeval::f_polyeval4;
35
36// Generated by SageMath with:
37// print("[")
38// for i in range(128):
39// R = RealField(200)
40// r = R(2)**(-8) * ( R(2)**8 * (R(1) - R(2)**(-8)) / (R(1) + R(i)*R(2)**-7) ).ceil()
41// b = float((r.log10()*R(2)**43).round() * R(2)**-43)
42// c = (r.log10() - R(b))
43// if i == 0 or i == 127:
44// print("(", double_to_hex(0), ",", double_to_hex(0), "),")
45// else:
46// print("(", double_to_hex(-c), ",", double_to_hex(-b), "),")
47// print("];")
48// static LOG10_R_DD: [(u64, u64); 128] = [
49// (0x0000000000000000, 0x0000000000000000),
50// (0xbd16079e5269431b, 0x3f6be76bd77c0000),
51// (0xbcffac7d2ae08e4f, 0x3f7c03a80ae60000),
52// (0x3d2fabf59b5d80b8, 0x3f851824c7580000),
53// (0x3d2310272fe17537, 0x3f8c3d0837780000),
54// (0xbd165201646ebccd, 0x3f91b85d60450000),
55// (0xbd2e22fab794a816, 0x3f9559bd24070000),
56// (0x3d2421bab5f034a4, 0x3f9902c31d628000),
57// (0x3d2fedb4b594a31b, 0x3f9cb38fccd88000),
58// (0xbd026ac877784097, 0x3f9e8eeb09f30000),
59// (0xbd2df23ed2ebe477, 0x3fa125d0432ec000),
60// (0xbd000c12f7a1b586, 0x3fa30838cdc30000),
61// (0x3d083662f181f53f, 0x3fa3faf7c6630000),
62// (0x3d22951bb9cd2fb7, 0x3fa5e3966b7e8000),
63// (0x3d1fae96a708581e, 0x3fa7d070145f4000),
64// (0x3d20744e2ea4f128, 0x3fa8c878eeb04000),
65// (0xbcfb0197d2cb982e, 0x3faabbcebd850000),
66// (0xbd2b1aad42a57f54, 0x3fabb7209d1e4000),
67// (0xbd240bcd23c3e44c, 0x3fadb11ed766c000),
68// (0xbcf626d2c723bf3b, 0x3faeafd05035c000),
69// (0x3cf77c3e779b9fcf, 0x3fb0585283764000),
70// (0x3cef3735158d42c3, 0x3fb0d966cc650000),
71// (0xbd2d227d61f9e88d, 0x3fb1dd5460c8c000),
72// (0xbcdf74be7c4de292, 0x3fb2603072a26000),
73// (0xbd1df5de49ddb160, 0x3fb367ba3aaa2000),
74// (0xbd1e5e3b38ac267a, 0x3fb3ec6ad5408000),
75// (0x3d275da8a5871b9a, 0x3fb4f7aad9bbc000),
76// (0x3d2ef5f3c10990f4, 0x3fb57e3d47c3a000),
77// (0x3d17c3cf23a17d9f, 0x3fb605735ee98000),
78// (0xbd141149840eaa65, 0x3fb715d0ce368000),
79// (0xbd1ff281b9601ce6, 0x3fb79efb57b10000),
80// (0x3d00a581f3edc493, 0x3fb828cfed29a000),
81// (0xbcf80743406505e6, 0x3fb93e7de0fc4000),
82// (0xbce76b169f6b4949, 0x3fb9ca5aa172a000),
83// (0xbd0bc8889d0fe1a0, 0x3fba56e8325f6000),
84// (0xbd25e950adf89934, 0x3fbae4285509a000),
85// (0xbd203ad4133e8c4c, 0x3fbb721cd1716000),
86// (0xbd2ddd18dedb6656, 0x3fbc902a19e66000),
87// (0x3d05e533080ecf32, 0x3fbd204698cb4000),
88// (0x3d27e865b8783768, 0x3fbdb11ed766a000),
89// (0x3d25e50ff38d4de9, 0x3fbe42b4c16ca000),
90// (0x3d25f7ef576ada0c, 0x3fbed50a4a26e000),
91// (0xbd1ff229f20ed3d2, 0x3fbffbfc2bbc8000),
92// (0xbd26f30673aae7ef, 0x3fc0484e4942b000),
93// (0x3d2dae5ed5e3f34c, 0x3fc093025a199000),
94// (0xbd23eea49e637bb3, 0x3fc0de1b56357000),
95// (0x3d182c6326f70b35, 0x3fc1299a4fb3e000),
96// (0x3d2f04d633b79054, 0x3fc175805d158000),
97// (0x3cf8b891b6d05a73, 0x3fc1c1ce9955c000),
98// (0xbcc35ca658049a0a, 0x3fc20e8624039000),
99// (0x3d2ff081a4e81f0b, 0x3fc25ba8215af000),
100// (0x3d21e3f04f63ee01, 0x3fc2a935ba5f1000),
101// (0xbd2e1471e5cb397e, 0x3fc2f7301cf4f000),
102// (0xbd25bd54fd6eb7d9, 0x3fc345987bfef000),
103// (0x3d1fe93f791a7264, 0x3fc394700f795000),
104// (0xbd28aebce3ec8738, 0x3fc3e3b814974000),
105// (0x3d2b0722aa2559f2, 0x3fc43371cde07000),
106// (0xbd1bcb784188a058, 0x3fc4839e83507000),
107// (0x3d220ca9a2bcc728, 0x3fc4d43f8275a000),
108// (0x3d2bb95ec8fd8f68, 0x3fc525561e925000),
109// (0x3cf4e036062e2e73, 0x3fc576e3b0bde000),
110// (0xbcce560b5e7a02b4, 0x3fc5c8e998073000),
111// (0x3ce1e472c751a4c0, 0x3fc61b6939983000),
112// (0xbcf1116ad28239b0, 0x3fc66e6400da4000),
113// (0x3d19ad1d9e405fb9, 0x3fc6c1db5f9bb000),
114// (0x3d19ad1d9e405fb9, 0x3fc6c1db5f9bb000),
115// (0xbd241149840eaa65, 0x3fc715d0ce368000),
116// (0x3d2bfb1de334b1cb, 0x3fc76a45cbb7e000),
117// (0xbcf9fc01708d86e7, 0x3fc7bf3bde09a000),
118// (0x3d24adaf7fe992b6, 0x3fc814b4921bd000),
119// (0xbd1c0b434f5d2b7a, 0x3fc86ab17c10c000),
120// (0xbd1c0b434f5d2b7a, 0x3fc86ab17c10c000),
121// (0x3d24c7acd659652f, 0x3fc8c13437695000),
122// (0x3d23e99da1b485cf, 0x3fc9183e67339000),
123// (0xbd1fb7d8e74e02a0, 0x3fc96fd1b63a0000),
124// (0x3d17c9690248757e, 0x3fc9c7efd734a000),
125// (0xbcb0cee0ed4ca7e9, 0x3fca209a84fbd000),
126// (0xbcb0cee0ed4ca7e9, 0x3fca209a84fbd000),
127// (0x3d0d924eb1fec6c5, 0x3fca79d382bc2000),
128// (0x3cefe6eb5a4df1df, 0x3fcad39c9c2c6000),
129// (0x3d14c9cce6dd603b, 0x3fcb2df7a5c50000),
130// (0x3d14c9cce6dd603b, 0x3fcb2df7a5c50000),
131// (0xbd2a01b9bb0ebf1b, 0x3fcb88e67cf98000),
132// (0x3d22ee896e06dbe8, 0x3fcbe46b08735000),
133// (0xbceff2381071d23e, 0x3fcc4087384f5000),
134// (0xbceff2381071d23e, 0x3fcc4087384f5000),
135// (0xbd22f9fd61140aa6, 0x3fcc9d3d065c6000),
136// (0xbd22386ad8b819b8, 0x3fccfa8e765cc000),
137// (0xbd22386ad8b819b8, 0x3fccfa8e765cc000),
138// (0x3d2d63da3ac9f0d5, 0x3fcd587d96494000),
139// (0x3d2fcc16f3ba09cb, 0x3fcdb70c7e96e000),
140// (0x3d2fcc16f3ba09cb, 0x3fcdb70c7e96e000),
141// (0xbd2cc52c1ba2d838, 0x3fce163d527e7000),
142// (0x3d1f97877007c127, 0x3fce76124046b000),
143// (0x3d1f97877007c127, 0x3fce76124046b000),
144// (0x3d0fbd42fdc335fd, 0x3fced68d81919000),
145// (0xbd2cbc0789055864, 0x3fcf37b15bab1000),
146// (0xbd2cbc0789055864, 0x3fcf37b15bab1000),
147// (0x3d22737df7f29668, 0x3fcf99801fdb7000),
148// (0xbd2ff229f20ed3d2, 0x3fcffbfc2bbc8000),
149// (0xbd2ff229f20ed3d2, 0x3fcffbfc2bbc8000),
150// (0x3d100809c16ae3ec, 0x3fd02f93f4c87000),
151// (0xbd2aa1358e87a6b4, 0x3fd06182e84fd800),
152// (0xbd2aa1358e87a6b4, 0x3fd06182e84fd800),
153// (0xbcff171b2c5f2274, 0x3fd093cc32c91000),
154// (0xbcff171b2c5f2274, 0x3fd093cc32c91000),
155// (0xbd242d6ae59e7d9c, 0x3fd0c6711d6ac000),
156// (0x3d2eac1871dbdbbf, 0x3fd0f972f87ff000),
157// (0x3d2eac1871dbdbbf, 0x3fd0f972f87ff000),
158// (0x3d1feed957c16a44, 0x3fd12cd31b9c9800),
159// (0x3d1feed957c16a44, 0x3fd12cd31b9c9800),
160// (0x3d1a6023a51d15b6, 0x3fd16092e5d3a800),
161// (0x3d2cefc5208422d8, 0x3fd194b3bdef6800),
162// (0x3d2cefc5208422d8, 0x3fd194b3bdef6800),
163// (0xbc81d4f1e8c2daff, 0x3fd1c93712abc800),
164// (0xbc81d4f1e8c2daff, 0x3fd1c93712abc800),
165// (0x3d140eb9b53054c3, 0x3fd1fe1e5af2c000),
166// (0x3d140eb9b53054c3, 0x3fd1fe1e5af2c000),
167// (0x3d29bbc8038401fc, 0x3fd2336b161b3000),
168// (0x3d29bbc8038401fc, 0x3fd2336b161b3000),
169// (0x3cf09ca54daae9f9, 0x3fd2691ecc29f000),
170// (0x3cf09ca54daae9f9, 0x3fd2691ecc29f000),
171// (0x3cf2b528446968a4, 0x3fd29f3b0e155800),
172// (0x3cf2b528446968a4, 0x3fd29f3b0e155800),
173// (0xbd14548507c3dd04, 0x3fd2d5c1760b8800),
174// (0xbd14548507c3dd04, 0x3fd2d5c1760b8800),
175// (0xbd1db59b99249f3a, 0x3fd30cb3a7bb3800),
176// (0x0000000000000000, 0x0000000000000000),
177// ];
178
179// Generated by Sollya with:
180// for i from 0 to 127 do {
181// r = 2^-8 * ceil( 2^8 * (1 - 2^(-8)) / (1 + i*2^-7) );
182// b = nearestint(log(r)*2^43) * 2^-43;
183// c = round(log(r) - b, D, RN);
184// print("{", -c, ",", -b, "},");
185// };
186// We replace LOG_R[0] with log10(1.0) == 0.0
187pub(crate) static LOG_R_DD: [(u64, u64); 128] = [
188 (0x0000000000000000, 0x0000000000000000),
189 (0xbd10c76b999d2be8, 0x3f80101575890000),
190 (0xbd23dc5b06e2f7d2, 0x3f90205658938000),
191 (0xbd2aa0ba325a0c34, 0x3f98492528c90000),
192 (0x3d0111c05cf1d753, 0x3fa0415d89e74000),
193 (0xbd2c167375bdfd28, 0x3fa466aed42e0000),
194 (0xbd197995d05a267d, 0x3fa894aa149fc000),
195 (0xbd1a68f247d82807, 0x3faccb73cdddc000),
196 (0xbd17e5dd7009902c, 0x3fb08598b59e4000),
197 (0xbd25325d560d9e9b, 0x3fb1973bd1466000),
198 (0x3d2cc85ea5db4ed7, 0x3fb3bdf5a7d1e000),
199 (0xbd2c69063c5d1d1e, 0x3fb5e95a4d97a000),
200 (0x3cec1e8da99ded32, 0x3fb700d30aeac000),
201 (0x3d23115c3abd47da, 0x3fb9335e5d594000),
202 (0xbd1390802bf768e5, 0x3fbb6ac88dad6000),
203 (0x3d2646d1c65aacd3, 0x3fbc885801bc4000),
204 (0xbd2dc068afe645e0, 0x3fbec739830a2000),
205 (0xbd2534d64fa10afd, 0x3fbfe89139dbe000),
206 (0x3d21ef78ce2d07f2, 0x3fc1178e8227e000),
207 (0x3d2ca78e44389934, 0x3fc1aa2b7e23f000),
208 (0x3d039d6ccb81b4a1, 0x3fc2d1610c868000),
209 (0x3cc62fa8234b7289, 0x3fc365fcb0159000),
210 (0x3d25837954fdb678, 0x3fc4913d8333b000),
211 (0x3d2633e8e5697dc7, 0x3fc527e5e4a1b000),
212 (0x3d19cf8b2c3c2e78, 0x3fc6574ebe8c1000),
213 (0xbd25118de59c21e1, 0x3fc6f0128b757000),
214 (0x3d1e0ddb9a631e83, 0x3fc823c16551a000),
215 (0xbd073d54aae92cd1, 0x3fc8beafeb390000),
216 (0x3d07f22858a0ff6f, 0x3fc95a5adcf70000),
217 (0xbd28724350562169, 0x3fca93ed3c8ae000),
218 (0xbd0c358d4eace1aa, 0x3fcb31d8575bd000),
219 (0xbd2d4bc4595412b6, 0x3fcbd087383be000),
220 (0xbd084a7e75b6f6e4, 0x3fcd1037f2656000),
221 (0xbd2aff2af715b035, 0x3fcdb13db0d49000),
222 (0x3cc212276041f430, 0x3fce530effe71000),
223 (0xbcca211565bb8e11, 0x3fcef5ade4dd0000),
224 (0x3d1bcbecca0cdf30, 0x3fcf991c6cb3b000),
225 (0x3cf89cdb16ed4e91, 0x3fd07138604d5800),
226 (0x3d27188b163ceae9, 0x3fd0c42d67616000),
227 (0xbd2c210e63a5f01c, 0x3fd1178e8227e800),
228 (0x3d2b9acdf7a51681, 0x3fd16b5ccbacf800),
229 (0x3d2ca6ed5147bdb7, 0x3fd1bf99635a6800),
230 (0x3d2c93c1df5bb3b6, 0x3fd269621134d800),
231 (0x3d2a9cfa4a5004f4, 0x3fd2bef07cdc9000),
232 (0xbd28e27ad3213cb8, 0x3fd314f1e1d36000),
233 (0x3d116ecdb0f177c8, 0x3fd36b6776be1000),
234 (0x3d183b54b606bd5c, 0x3fd3c25277333000),
235 (0x3d08e436ec90e09d, 0x3fd419b423d5e800),
236 (0xbd2f27ce0967d675, 0x3fd4718dc271c800),
237 (0xbd2e20891b0ad8a4, 0x3fd4c9e09e173000),
238 (0x3d2ebe708164c759, 0x3fd522ae0738a000),
239 (0x3d1fadedee5d40ef, 0x3fd57bf753c8d000),
240 (0xbd0a0b2a08a465dc, 0x3fd5d5bddf596000),
241 (0xbd2db623e731ae00, 0x3fd630030b3ab000),
242 (0x3d20a0d32756eba0, 0x3fd68ac83e9c6800),
243 (0x3d1721657c222d87, 0x3fd6e60ee6af1800),
244 (0x3d2d8b0949dc60b3, 0x3fd741d876c67800),
245 (0x3d29ec7d2efd1778, 0x3fd79e26687cf800),
246 (0xbd272090c812566a, 0x3fd7fafa3bd81800),
247 (0x3d2fd56f3333778a, 0x3fd85855776dc800),
248 (0xbd205ae1e5e70470, 0x3fd8b639a88b3000),
249 (0xbd1766b52ee6307d, 0x3fd914a8635bf800),
250 (0xbd152313a502d9f0, 0x3fd973a343135800),
251 (0xbd26279e10d0c0b0, 0x3fd9d32bea15f000),
252 (0x3d23c6457f9d79f5, 0x3fda33440224f800),
253 (0x3d23c6457f9d79f5, 0x3fda33440224f800),
254 (0x3d1e36f2bea77a5d, 0x3fda93ed3c8ad800),
255 (0xbd217cc552774458, 0x3fdaf5295248d000),
256 (0x3d1095252d841995, 0x3fdb56fa04462800),
257 (0x3d27d85bf40a666d, 0x3fdbb9611b80e000),
258 (0x3d2cec807fe8e180, 0x3fdc1c60693fa000),
259 (0x3d2cec807fe8e180, 0x3fdc1c60693fa000),
260 (0xbd29b6ddc15249ae, 0x3fdc7ff9c7455800),
261 (0xbd0797c33ec7a6b0, 0x3fdce42f18064800),
262 (0x3d235bafe9a767a8, 0x3fdd490246def800),
263 (0xbd1ea42d60dc616a, 0x3fddae75484c9800),
264 (0xbd1326b207322938, 0x3fde148a1a272800),
265 (0xbd1326b207322938, 0x3fde148a1a272800),
266 (0xbd2465505372bd08, 0x3fde7b42c3ddb000),
267 (0x3d2f27f45a470251, 0x3fdee2a156b41000),
268 (0x3d12cde56f014a8b, 0x3fdf4aa7ee031800),
269 (0x3d12cde56f014a8b, 0x3fdf4aa7ee031800),
270 (0x3d0085fa3c164935, 0x3fdfb358af7a4800),
271 (0xbd053ba3b1727b1c, 0x3fe00e5ae5b20800),
272 (0xbd04c45fe79539e0, 0x3fe04360be760400),
273 (0xbd04c45fe79539e0, 0x3fe04360be760400),
274 (0x3d26812241edf5fd, 0x3fe078bf0533c400),
275 (0x3d1f486b887e7e27, 0x3fe0ae76e2d05400),
276 (0x3d1f486b887e7e27, 0x3fe0ae76e2d05400),
277 (0x3d1c299807801742, 0x3fe0e4898611cc00),
278 (0xbd258647bb9ddcb2, 0x3fe11af823c75c00),
279 (0xbd258647bb9ddcb2, 0x3fe11af823c75c00),
280 (0xbd2edd97a293ae49, 0x3fe151c3f6f29800),
281 (0x3d14cc4ef8ab4650, 0x3fe188ee40f23c00),
282 (0x3d14cc4ef8ab4650, 0x3fe188ee40f23c00),
283 (0x3cccacdeed70e667, 0x3fe1c07849ae6000),
284 (0xbd2a7242c9fe81d3, 0x3fe1f8635fc61800),
285 (0xbd2a7242c9fe81d3, 0x3fe1f8635fc61800),
286 (0x3d12fc066e48667b, 0x3fe230b0d8bebc00),
287 (0xbd0b61f105226250, 0x3fe269621134dc00),
288 (0xbd0b61f105226250, 0x3fe269621134dc00),
289 (0x3d206d2be797882d, 0x3fe2a2786d0ec000),
290 (0xbd17a6e507b9dc11, 0x3fe2dbf557b0e000),
291 (0xbd17a6e507b9dc11, 0x3fe2dbf557b0e000),
292 (0xbd274e93c5a0ed9c, 0x3fe315da44340800),
293 (0xbd274e93c5a0ed9c, 0x3fe315da44340800),
294 (0x3d10b83f9527e6ac, 0x3fe35028ad9d8c00),
295 (0xbd218b7abb5569a4, 0x3fe38ae217197800),
296 (0xbd218b7abb5569a4, 0x3fe38ae217197800),
297 (0xbd02b7367cfe13c2, 0x3fe3c6080c36c000),
298 (0xbd02b7367cfe13c2, 0x3fe3c6080c36c000),
299 (0xbd26ce7930f0c74c, 0x3fe4019c2125cc00),
300 (0xbcfd984f481051f7, 0x3fe43d9ff2f92400),
301 (0xbcfd984f481051f7, 0x3fe43d9ff2f92400),
302 (0xbd22cb6af94d60aa, 0x3fe47a1527e8a400),
303 (0xbd22cb6af94d60aa, 0x3fe47a1527e8a400),
304 (0x3cef7115ed4c541c, 0x3fe4b6fd6f970c00),
305 (0x3cef7115ed4c541c, 0x3fe4b6fd6f970c00),
306 (0xbd2e6c516d93b8fb, 0x3fe4f45a835a5000),
307 (0xbd2e6c516d93b8fb, 0x3fe4f45a835a5000),
308 (0x3d05ccc45d257531, 0x3fe5322e26867800),
309 (0x3d05ccc45d257531, 0x3fe5322e26867800),
310 (0x3d09980bff3303dd, 0x3fe5707a26bb8c00),
311 (0x3d09980bff3303dd, 0x3fe5707a26bb8c00),
312 (0x3d2dfa63ac10c9fb, 0x3fe5af405c364800),
313 (0x3d2dfa63ac10c9fb, 0x3fe5af405c364800),
314 (0x3d2202380cda46be, 0x3fe5ee82aa241800),
315 (0x0000000000000000, 0x0000000000000000),
316];
317
318/// Logarithm of base 10
319///
320/// Max found ULP 0.5
321pub fn f_log10(x: f64) -> f64 {
322 let mut x_u = x.to_bits();
323
324 const E_BIAS: u64 = (1u64 << (11 - 1u64)) - 1u64;
325
326 let mut x_e: i64 = -(E_BIAS as i64);
327
328 const MAX_NORMAL: u64 = f64::to_bits(f64::MAX);
329
330 if x_u == 1f64.to_bits() {
331 // log2(1.0) = +0.0
332 return 0.0;
333 }
334 if x_u < min_normal_f64().to_bits() || x_u > MAX_NORMAL {
335 if x == 0.0 {
336 return f64::NEG_INFINITY;
337 }
338 if x < 0. || x.is_nan() {
339 return f64::NAN;
340 }
341 if x.is_infinite() || x.is_nan() {
342 return x + x;
343 }
344 // Normalize denormal inputs.
345 x_u = (x * f64::from_bits(0x4330000000000000)).to_bits();
346 x_e -= 52;
347 }
348
349 // log2(x) = log2(2^x_e * x_m)
350 // = x_e + log2(x_m)
351 // Range reduction for log2(x_m):
352 // For each x_m, we would like to find r such that:
353 // -2^-8 <= r * x_m - 1 < 2^-7
354 let shifted = (x_u >> 45) as i64;
355 let index = shifted & 0x7F;
356 let r = f64::from_bits(crate::logs::log2::LOG_RANGE_REDUCTION[index as usize]);
357
358 // Add unbiased exponent. Add an extra 1 if the 8 leading fractional bits are
359 // all 1's.
360 x_e = x_e.wrapping_add(x_u.wrapping_add(1u64 << 45).wrapping_shr(52) as i64);
361 let e_x = x_e as f64;
362
363 const LOG_2_HI: f64 = f64::from_bits(0x3fe62e42fefa3800);
364 const LOG_2_LO: f64 = f64::from_bits(0x3d2ef35793c76730);
365
366 let log_r_dd = LOG_R_DD[index as usize];
367
368 // hi is exact
369 let hi = f_fmla(e_x, LOG_2_HI, f64::from_bits(log_r_dd.1));
370 // lo errors ~ e_x * LSB(LOG_2_LO) + LSB(LOG_R[index].lo) + rounding err
371 // <= 2 * (e_x * LSB(LOG_2_LO) + LSB(LOG_R[index].lo))
372 let lo = f_fmla(e_x, LOG_2_LO, f64::from_bits(log_r_dd.0));
373
374 // Set m = 1.mantissa.
375 let x_m = (x_u & 0x000F_FFFF_FFFF_FFFFu64) | 0x3FF0_0000_0000_0000u64;
376 let m = f64::from_bits(x_m);
377
378 let u;
379 #[cfg(any(
380 all(
381 any(target_arch = "x86", target_arch = "x86_64"),
382 target_feature = "fma"
383 ),
384 target_arch = "aarch64"
385 ))]
386 {
387 u = f_fmla(r, m, -1.0); // exact
388 }
389 #[cfg(not(any(
390 all(
391 any(target_arch = "x86", target_arch = "x86_64"),
392 target_feature = "fma"
393 ),
394 target_arch = "aarch64"
395 )))]
396 {
397 use crate::logs::log2::LOG_CD;
398 let c_m = x_m & 0x3FFF_E000_0000_0000u64;
399 let c = f64::from_bits(c_m);
400 u = f_fmla(r, m - c, f64::from_bits(LOG_CD[index as usize])); // exact
401 }
402
403 let u_sq = u * u;
404 // Degree-7 minimax polynomial
405 let p0 = f_fmla(
406 u,
407 f64::from_bits(LOG_COEFFS[1]),
408 f64::from_bits(LOG_COEFFS[0]),
409 );
410 let p1 = f_fmla(
411 u,
412 f64::from_bits(LOG_COEFFS[3]),
413 f64::from_bits(LOG_COEFFS[2]),
414 );
415 let p2 = f_fmla(
416 u,
417 f64::from_bits(LOG_COEFFS[5]),
418 f64::from_bits(LOG_COEFFS[4]),
419 );
420 let p = f_polyeval4(u_sq, lo, p0, p1, p2);
421
422 // Exact sum:
423 // r1.hi + r1.lo = e_x * log(2)_hi - log(r)_hi + u
424 let mut r1 = DoubleDouble::from_exact_add(hi, u);
425 r1.lo += p;
426
427 // Quick double-double multiplication:
428 // r2.hi + r2.lo ~ r1 * log10(e),
429 // with error bounded by:
430 // 4*ulp( ulp(r2.hi) )
431 const LOG10_E: DoubleDouble = DoubleDouble::new(
432 f64::from_bits(0x3c695355baaafad3),
433 f64::from_bits(0x3fdbcb7b1526e50e),
434 );
435 let r2 = DoubleDouble::quick_mult(r1, LOG10_E);
436
437 const HI_ERR: f64 = f64::from_bits(0x3aa0000000000000);
438
439 // Extra errors from P is from using x^2 to reduce evaluation latency.
440 const P_ERR: f64 = f64::from_bits(0x3cc0000000000000);
441
442 // Technicallly error of r1.lo is bounded by:
443 // |hi|*ulp(log(2)_lo) + C*ulp(u^2)
444 // To simplify the error computation a bit, we replace |hi|*ulp(log(2)_lo)
445 // with the upper bound: 2^11 * ulp(log(2)_lo) = 2^-85.
446 // Total error is bounded by ~ C * ulp(u^2) + 2^-85.
447 let err = f_fmla(u_sq, P_ERR, HI_ERR);
448
449 // Lower bound from the result
450 let left = r2.hi + (r2.lo - err);
451 // Upper bound from the result
452 let right = r2.hi + (r2.lo + err);
453
454 // Ziv's test if fast pass is accurate enough.
455 if left == right {
456 return left;
457 }
458
459 log10_dd_accurate(x)
460}
461
462#[cold]
463#[inline(never)]
464fn log10_dd_accurate(x: f64) -> f64 {
465 let r = log10_dd(x);
466 let err = f_fmla(
467 r.hi,
468 f64::from_bits(0x3b50000000000000), // 2^-74
469 f64::from_bits(0x3990000000000000), // 2^-102
470 );
471 let ub = r.hi + (r.lo + err);
472 let lb = r.hi + (r.lo - err);
473 if ub == lb {
474 return r.to_f64();
475 }
476 log10_dd_accurate_slow(x)
477}
478
479#[cold]
480#[inline(never)]
481fn log10_dd_accurate_slow(x: f64) -> f64 {
482 log10_td(x).to_f64()
483}
484
485#[cfg(test)]
486mod tests {
487 use super::*;
488
489 #[test]
490 fn test_log10d() {
491 assert_eq!(f_log10(0.35), -0.4559319556497244);
492 assert_eq!(f_log10(0.9), -0.045757490560675115);
493 assert_eq!(f_log10(10.), 1.);
494 assert_eq!(f_log10(0.), f64::NEG_INFINITY);
495 assert!(f_log10(-1.).is_nan());
496 assert!(f_log10(f64::NAN).is_nan());
497 assert!(f_log10(f64::NEG_INFINITY).is_nan());
498 assert_eq!(f_log10(f64::INFINITY), f64::INFINITY);
499 }
500}