/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 | | } |