Coverage Report

Created: 2025-11-05 08:08

next uncovered line (L), next uncovered region (R), next uncovered branch (B)
/rust/registry/src/index.crates.io-1949cf8c6b5b557f/pxfm-0.1.25/src/logs/log10.rs
Line
Count
Source
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
 */
29
use crate::common::{f_fmla, min_normal_f64};
30
use crate::double_double::DoubleDouble;
31
use crate::logs::log2::LOG_COEFFS;
32
use crate::logs::log10dd::log10_dd;
33
use crate::logs::log10td::log10_td;
34
use 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
187
pub(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
321
0
pub fn f_log10(x: f64) -> f64 {
322
0
    let mut x_u = x.to_bits();
323
324
    const E_BIAS: u64 = (1u64 << (11 - 1u64)) - 1u64;
325
326
0
    let mut x_e: i64 = -(E_BIAS as i64);
327
328
    const MAX_NORMAL: u64 = f64::to_bits(f64::MAX);
329
330
0
    if x_u == 1f64.to_bits() {
331
        // log2(1.0) = +0.0
332
0
        return 0.0;
333
0
    }
334
0
    if x_u < min_normal_f64().to_bits() || x_u > MAX_NORMAL {
335
0
        if x == 0.0 {
336
0
            return f64::NEG_INFINITY;
337
0
        }
338
0
        if x < 0. || x.is_nan() {
339
0
            return f64::NAN;
340
0
        }
341
0
        if x.is_infinite() || x.is_nan() {
342
0
            return x + x;
343
0
        }
344
        // Normalize denormal inputs.
345
0
        x_u = (x * f64::from_bits(0x4330000000000000)).to_bits();
346
0
        x_e -= 52;
347
0
    }
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
0
    let shifted = (x_u >> 45) as i64;
355
0
    let index = shifted & 0x7F;
356
0
    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
0
    x_e = x_e.wrapping_add(x_u.wrapping_add(1u64 << 45).wrapping_shr(52) as i64);
361
0
    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
0
    let log_r_dd = LOG_R_DD[index as usize];
367
368
    // hi is exact
369
0
    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
0
    let lo = f_fmla(e_x, LOG_2_LO, f64::from_bits(log_r_dd.0));
373
374
    // Set m = 1.mantissa.
375
0
    let x_m = (x_u & 0x000F_FFFF_FFFF_FFFFu64) | 0x3FF0_0000_0000_0000u64;
376
0
    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
0
        let c_m = x_m & 0x3FFF_E000_0000_0000u64;
399
0
        let c = f64::from_bits(c_m);
400
0
        u = f_fmla(r, m - c, f64::from_bits(LOG_CD[index as usize])); // exact
401
    }
402
403
0
    let u_sq = u * u;
404
    // Degree-7 minimax polynomial
405
0
    let p0 = f_fmla(
406
0
        u,
407
0
        f64::from_bits(LOG_COEFFS[1]),
408
0
        f64::from_bits(LOG_COEFFS[0]),
409
    );
410
0
    let p1 = f_fmla(
411
0
        u,
412
0
        f64::from_bits(LOG_COEFFS[3]),
413
0
        f64::from_bits(LOG_COEFFS[2]),
414
    );
415
0
    let p2 = f_fmla(
416
0
        u,
417
0
        f64::from_bits(LOG_COEFFS[5]),
418
0
        f64::from_bits(LOG_COEFFS[4]),
419
    );
420
0
    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
0
    let mut r1 = DoubleDouble::from_exact_add(hi, u);
425
0
    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
0
    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
0
    let err = f_fmla(u_sq, P_ERR, HI_ERR);
448
449
    // Lower bound from the result
450
0
    let left = r2.hi + (r2.lo - err);
451
    // Upper bound from the result
452
0
    let right = r2.hi + (r2.lo + err);
453
454
    // Ziv's test if fast pass is accurate enough.
455
0
    if left == right {
456
0
        return left;
457
0
    }
458
459
0
    log10_dd_accurate(x)
460
0
}
461
462
#[cold]
463
#[inline(never)]
464
0
fn log10_dd_accurate(x: f64) -> f64 {
465
0
    let r = log10_dd(x);
466
0
    let err = f_fmla(
467
0
        r.hi,
468
0
        f64::from_bits(0x3b50000000000000), // 2^-74
469
0
        f64::from_bits(0x3990000000000000), // 2^-102
470
    );
471
0
    let ub = r.hi + (r.lo + err);
472
0
    let lb = r.hi + (r.lo - err);
473
0
    if ub == lb {
474
0
        return r.to_f64();
475
0
    }
476
0
    log10_dd_accurate_slow(x)
477
0
}
478
479
#[cold]
480
#[inline(never)]
481
0
fn log10_dd_accurate_slow(x: f64) -> f64 {
482
0
    log10_td(x).to_f64()
483
0
}
484
485
#[cfg(test)]
486
mod 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
}