/rust/registry/src/index.crates.io-1949cf8c6b5b557f/pxfm-0.1.25/src/logs/log2.rs

Source
/*
 * // Copyright (c) Radzivon Bartoshyk 7/2025. All rights reserved.
 * //
 * // Redistribution and use in source and binary forms, with or without modification,
 * // are permitted provided that the following conditions are met:
 * //
 * // 1.  Redistributions of source code must retain the above copyright notice, this
 * // list of conditions and the following disclaimer.
 * //
 * // 2.  Redistributions in binary form must reproduce the above copyright notice,
 * // this list of conditions and the following disclaimer in the documentation
 * // and/or other materials provided with the distribution.
 * //
 * // 3.  Neither the name of the copyright holder nor the names of its
 * // contributors may be used to endorse or promote products derived from
 * // this software without specific prior written permission.
 * //
 * // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * // DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
 * // FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
 * // CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
 * // OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 * // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */
use crate::common::{f_fmla, min_normal_f64};
use crate::double_double::DoubleDouble;
use crate::logs::log2dd::log2_dd;
use crate::logs::log2td::log2_td;
use crate::polyeval::f_polyeval6;

pub(crate) static LOG_RANGE_REDUCTION: [u64; 128] = [
    0x3ff0000000000000,
    0x3fefc00000000000,
    0x3fef800000000000,
    0x3fef400000000000,
    0x3fef000000000000,
    0x3feec00000000000,
    0x3fee800000000000,
    0x3fee400000000000,
    0x3fee000000000000,
    0x3fede00000000000,
    0x3feda00000000000,
    0x3fed600000000000,
    0x3fed400000000000,
    0x3fed000000000000,
    0x3fecc00000000000,
    0x3feca00000000000,
    0x3fec600000000000,
    0x3fec400000000000,
    0x3fec000000000000,
    0x3febe00000000000,
    0x3feba00000000000,
    0x3feb800000000000,
    0x3feb400000000000,
    0x3feb200000000000,
    0x3feae00000000000,
    0x3feac00000000000,
    0x3fea800000000000,
    0x3fea600000000000,
    0x3fea400000000000,
    0x3fea000000000000,
    0x3fe9e00000000000,
    0x3fe9c00000000000,
    0x3fe9800000000000,
    0x3fe9600000000000,
    0x3fe9400000000000,
    0x3fe9200000000000,
    0x3fe9000000000000,
    0x3fe8c00000000000,
    0x3fe8a00000000000,
    0x3fe8800000000000,
    0x3fe8600000000000,
    0x3fe8400000000000,
    0x3fe8000000000000,
    0x3fe7e00000000000,
    0x3fe7c00000000000,
    0x3fe7a00000000000,
    0x3fe7800000000000,
    0x3fe7600000000000,
    0x3fe7400000000000,
    0x3fe7200000000000,
    0x3fe7000000000000,
    0x3fe6e00000000000,
    0x3fe6c00000000000,
    0x3fe6a00000000000,
    0x3fe6800000000000,
    0x3fe6600000000000,
    0x3fe6400000000000,
    0x3fe6200000000000,
    0x3fe6000000000000,
    0x3fe5e00000000000,
    0x3fe5c00000000000,
    0x3fe5a00000000000,
    0x3fe5800000000000,
    0x3fe5600000000000,
    0x3fe5400000000000,
    0x3fe5400000000000,
    0x3fe5200000000000,
    0x3fe5000000000000,
    0x3fe4e00000000000,
    0x3fe4c00000000000,
    0x3fe4a00000000000,
    0x3fe4a00000000000,
    0x3fe4800000000000,
    0x3fe4600000000000,
    0x3fe4400000000000,
    0x3fe4200000000000,
    0x3fe4000000000000,
    0x3fe4000000000000,
    0x3fe3e00000000000,
    0x3fe3c00000000000,
    0x3fe3a00000000000,
    0x3fe3a00000000000,
    0x3fe3800000000000,
    0x3fe3600000000000,
    0x3fe3400000000000,
    0x3fe3400000000000,
    0x3fe3200000000000,
    0x3fe3000000000000,
    0x3fe3000000000000,
    0x3fe2e00000000000,
    0x3fe2c00000000000,
    0x3fe2c00000000000,
    0x3fe2a00000000000,
    0x3fe2800000000000,
    0x3fe2800000000000,
    0x3fe2600000000000,
    0x3fe2400000000000,
    0x3fe2400000000000,
    0x3fe2200000000000,
    0x3fe2000000000000,
    0x3fe2000000000000,
    0x3fe1e00000000000,
    0x3fe1c00000000000,
    0x3fe1c00000000000,
    0x3fe1a00000000000,
    0x3fe1a00000000000,
    0x3fe1800000000000,
    0x3fe1600000000000,
    0x3fe1600000000000,
    0x3fe1400000000000,
    0x3fe1400000000000,
    0x3fe1200000000000,
    0x3fe1000000000000,
    0x3fe1000000000000,
    0x3fe0e00000000000,
    0x3fe0e00000000000,
    0x3fe0c00000000000,
    0x3fe0c00000000000,
    0x3fe0a00000000000,
    0x3fe0a00000000000,
    0x3fe0800000000000,
    0x3fe0800000000000,
    0x3fe0600000000000,
    0x3fe0600000000000,
    0x3fe0400000000000,
    0x3fe0400000000000,
    0x3fe0200000000000,
    0x3fe0000000000000,
];

/**
Generated by SageMath:
```python
values = LOG_RANGE_REDUCTION

R = RealField(150)

def hex_to_float(h):
    return struct.unpack('>d', struct.pack('>Q', h))[0]

real_array = [R(hex_to_float(h)) for h in values]

for r in real_array:
    print_double_double("", -RealField(180)(r).log2())
```
**/
static LOG2_DD: [(u64, u64); 128] = [
    (0x0, 0x0),
    (0x3c26d746128b1857, 0x3f872c7ba20f7327),
    (0x3c2b2a41b08fbe06, 0x3f9743ee861f3556),
    (0x3c3491f06c085bc2, 0x3fa184b8e4c56af8),
    (0x3c477970e03f821c, 0x3fa77394c9d958d5),
    (0x3c0155660710eb2a, 0x3fad6ebd1f1febfe),
    (0x3c298c5452bbce74, 0x3fb1bb32a600549d),
    (0x3c2c141e66faaaad, 0x3fb4c560fe68af88),
    (0x3c59ced1447e30ad, 0x3fb7d60496cfbb4c),
    (0xbc321ba488a62577, 0x3fb960caf9abb7ca),
    (0xbc17936311889913, 0x3fbc7b528b70f1c5),
    (0xbc36cd4d2ae3a2f6, 0x3fbf9c95dc1d1165),
    (0x3c55d243efd93259, 0x3fc097e38ce60649),
    (0xbc5696e2866c718e, 0x3fc22dadc2ab3497),
    (0xbc67a6ed4e1b0936, 0x3fc3c6fb650cde51),
    (0xbc61562d61af73f8, 0x3fc494f863b8df35),
    (0x3c354fae008fbb59, 0x3fc633a8bf437ce1),
    (0xbc60798d1aa21694, 0x3fc7046031c79f85),
    (0x3c699aa6df8b7d83, 0x3fc8a8980abfbd32),
    (0xbc6e95734abd2fcc, 0x3fc97c1cb13c7ec1),
    (0xbc6cc865b3dd0dbb, 0x3fcb2602497d5346),
    (0x3c2bc0af7b82e7d7, 0x3fcbfc67a7fff4cc),
    (0xbc0ba8b1f646ab12, 0x3fcdac22d3e441d3),
    (0xbc6086fce864a1f6, 0x3fce857d3d361368),
    (0x3c7768994400ca0a, 0x3fd01d9bbcfa61d4),
    (0xbc53d56efe4338fe, 0x3fd08bce0d95fa38),
    (0x3c7c8d43e017579b, 0x3fd169c05363f158),
    (0x3c6ae9804237ec8e, 0x3fd1d982c9d52708),
    (0x3c734107c0e54aed, 0x3fd249cd2b13cd6c),
    (0x3c7968925e378d68, 0x3fd32bfee370ee68),
    (0x3c7fcad2f4710e00, 0x3fd39de8e1559f6f),
    (0xbc7c658d602e66b0, 0x3fd4106017c3eca3),
    (0x3c7e393a16b94b52, 0x3fd4f6fbb2cec598),
    (0xbc78d86531d55da2, 0x3fd56b22e6b578e5),
    (0x3c710b5b643a6ecb, 0x3fd5dfdcf1eeae0e),
    (0x3c7ac9080333c605, 0x3fd6552b49986277),
    (0x3c2b6d40900b2502, 0x3fd6cb0f6865c8ea),
    (0x3c68f89e2eb553b2, 0x3fd7b89f02cf2aad),
    (0x3c651d58525aad39, 0x3fd8304d90c11fd3),
    (0x3c799aa6df8b7d83, 0x3fd8a8980abfbd32),
    (0x3c7fdc46af571993, 0x3fd921800924dd3b),
    (0x3c7bca36fd02def0, 0x3fd99b072a96c6b2),
    (0x3c5817fd3b7d7e5d, 0x3fda8ff971810a5e),
    (0xbc45e13b838eba7d, 0x3fdb0b67f4f46810),
    (0xbc501d98c3531027, 0x3fdb877c57b1b070),
    (0x3c7edf515c63dd87, 0x3fdc043859e2fdb3),
    (0x3c6c4aec56233279, 0x3fdc819dc2d45fe4),
    (0x3c78a38b4175d665, 0x3fdcffae611ad12b),
    (0xbc6e15a52a31604a, 0x3fdd7e6c0abc3579),
    (0x3c438c8946414c6a, 0x3fddfdd89d586e2b),
    (0x3c73bed456b24ed1, 0x3fde7df5fe538ab3),
    (0x3c76d261f1753e0b, 0x3fdefec61b011f85),
    (0xbc79ca1a3202b3d7, 0x3fdf804ae8d0cd02),
    (0xbc87398fe685f171, 0x3fe0014332be0033),
    (0xbc89c32630008a1f, 0x3fe042bd4b9a7c99),
    (0xbc7f47806a0e4105, 0x3fe08494c66b8ef0),
    (0xbc48a33c25e8e226, 0x3fe0c6caaf0c5597),
    (0xbc73aec658457c41, 0x3fe1096015dee4da),
    (0xbc8fc7d7c3320aab, 0x3fe14c560fe68af9),
    (0xbc892ba145dcf40b, 0x3fe18fadb6e2d3c2),
    (0xbc7f9fb952bbbccc, 0x3fe1d368296b5255),
    (0xbc8568859c64022e, 0x3fe217868b0c37e8),
    (0xbc84828ddf1fb145, 0x3fe25c0a0463beb0),
    (0xbc8c348e4aab18b8, 0x3fe2a0f3c340705c),
    (0xbc83af881af2f3d9, 0x3fe2e644fac04fd8),
    (0xbc83af881af2f3d9, 0x3fe2e644fac04fd8),
    (0x3c8968925e378d68, 0x3fe32bfee370ee68),
    (0xbc869656a0ad70d4, 0x3fe37222bb70747c),
    (0x3c76d266d6cdc959, 0x3fe3b8b1c68fa6ed),
    (0xbc69575b04fa6fbd, 0x3fe3ffad4e74f1d6),
    (0x3c5b90132aeddb58, 0x3fe44716a2c08262),
    (0x3c5b90132aeddb58, 0x3fe44716a2c08262),
    (0xbc75e35482d13dc1, 0x3fe48eef19317991),
    (0xbc8ca44f1db913d3, 0x3fe4d7380dcc422d),
    (0x3c7817fd3b7d7e5d, 0x3fe51ff2e30214bc),
    (0x3c804613e33c06c9, 0x3fe5692101d9b4a6),
    (0xbc8fc9257edfe9b6, 0x3fe5b2c3da19723b),
    (0xbc8fc9257edfe9b6, 0x3fe5b2c3da19723b),
    (0x3c8149a1977b5b99, 0x3fe5fcdce2727ddb),
    (0xbc8b32266d92c0fe, 0x3fe6476d98ad990a),
    (0x3c821d90b84e7218, 0x3fe6927781d932a8),
    (0x3c821d90b84e7218, 0x3fe6927781d932a8),
    (0x3c7f6e91ad16ecff, 0x3fe6ddfc2a78fc63),
    (0xbc74a31ce1b7e328, 0x3fe729fd26b707c8),
    (0x3c6a7b47d2c352d9, 0x3fe7767c12967a45),
    (0x3c6a7b47d2c352d9, 0x3fe7767c12967a45),
    (0x3c821f9cb2cc5575, 0x3fe7c37a9227e7fb),
    (0x3c8dc572667587b1, 0x3fe810fa51bf65fd),
    (0x3c8dc572667587b1, 0x3fe810fa51bf65fd),
    (0xbc78f93e7aa3bdf8, 0x3fe85efd062c656d),
    (0x3c5b85a54d7ee2fd, 0x3fe8ad846cf369a4),
    (0x3c5b85a54d7ee2fd, 0x3fe8ad846cf369a4),
    (0x3c8bf1d926766301, 0x3fe8fc924c89ac84),
    (0x3c401ee1343fe7ca, 0x3fe94c287492c4db),
    (0x3c401ee1343fe7ca, 0x3fe94c287492c4db),
    (0x3c7fa0a62e6add1b, 0x3fe99c48be2063c8),
    (0x3c8022ddb71189c5, 0x3fe9ecf50bf43f13),
    (0x3c8022ddb71189c5, 0x3fe9ecf50bf43f13),
    (0x3c7ac7fc60a51031, 0x3fea3e2f4ac43f60),
    (0x3c6817fd3b7d7e5d, 0x3fea8ff971810a5e),
    (0x3c6817fd3b7d7e5d, 0x3fea8ff971810a5e),
    (0xbc83138e941643f7, 0x3feae255819f022d),
    (0x3c8e0ae0d3f8a58b, 0x3feb35458761d479),
    (0x3c8e0ae0d3f8a58b, 0x3feb35458761d479),
    (0x3c742b7579f0f8d4, 0x3feb88cb9a2ab521),
    (0x3c742b7579f0f8d4, 0x3feb88cb9a2ab521),
    (0x3c6a7610e40bd6ab, 0x3febdce9dcc96187),
    (0xbc80e5edaecee150, 0x3fec31a27dd00b4a),
    (0xbc80e5edaecee150, 0x3fec31a27dd00b4a),
    (0xbc831d962d3728cc, 0x3fec86f7b7ea4a89),
    (0xbc831d962d3728cc, 0x3fec86f7b7ea4a89),
    (0xbc857391924a6d9d, 0x3fecdcebd2373995),
    (0x3c78333ac7d9ebbb, 0x3fed338120a6dd9d),
    (0x3c78333ac7d9ebbb, 0x3fed338120a6dd9d),
    (0xbc86c0268890da53, 0x3fed8aba045b01c8),
    (0xbc86c0268890da53, 0x3fed8aba045b01c8),
    (0xbc859e7ba5d5ccc9, 0x3fede298ec0bac0d),
    (0xbc859e7ba5d5ccc9, 0x3fede298ec0bac0d),
    (0x3c60b07079619c47, 0x3fee3b20546f554a),
    (0x3c60b07079619c47, 0x3fee3b20546f554a),
    (0xbc8cc4d81bc25adf, 0x3fee9452c8a71028),
    (0xbc8cc4d81bc25adf, 0x3fee9452c8a71028),
    (0xbc776c0a2827d49a, 0x3feeee32e2aeccbf),
    (0xbc776c0a2827d49a, 0x3feeee32e2aeccbf),
    (0xbc8314dc4fc42302, 0x3fef48c34bd1e96f),
    (0xbc8314dc4fc42302, 0x3fef48c34bd1e96f),
    (0xbc817f8e37b00179, 0x3fefa406bd2443df),
    (0x0, 0x0),
];

#[cfg(not(any(
    all(
        any(target_arch = "x86", target_arch = "x86_64"),
        target_feature = "fma"
    ),
    target_arch = "aarch64"
)))]
pub(crate) static LOG_CD: [u64; 128] = [
    0x0000000000000000,
    0xbf10000000000000,
    0xbf30000000000000,
    0xbf42000000000000,
    0xbf50000000000000,
    0xbf59000000000000,
    0xbf62000000000000,
    0xbf68800000000000,
    0xbf70000000000000,
    0xbf49000000000000,
    0xbf5f000000000000,
    0xbf69c00000000000,
    0xbf30000000000000,
    0xbf5c000000000000,
    0xbf6b000000000000,
    0xbf45000000000000,
    0xbf64000000000000,
    0x3f10000000000000,
    0xbf60000000000000,
    0x3f3a000000000000,
    0xbf5e000000000000,
    0x3f38000000000000,
    0xbf61000000000000,
    0xbf00000000000000,
    0xbf66000000000000,
    0xbf4a000000000000,
    0xbf6e000000000000,
    0xbf5f800000000000,
    0xbf30000000000000,
    0xbf6c000000000000,
    0xbf5f000000000000,
    0xbf3c000000000000,
    0xbf70000000000000,
    0xbf65400000000000,
    0xbf56000000000000,
    0xbf24000000000000,
    0x3f50000000000000,
    0xbf68800000000000,
    0xbf60800000000000,
    0xbf52000000000000,
    0xbf30000000000000,
    0x3f42000000000000,
    0xbf70000000000000,
    0xbf6ac00000000000,
    0xbf66000000000000,
    0xbf61c00000000000,
    0xbf5c000000000000,
    0xbf55800000000000,
    0xbf50000000000000,
    0xbf47000000000000,
    0xbf40000000000000,
    0xbf36000000000000,
    0xbf30000000000000,
    0xbf2c000000000000,
    0xbf30000000000000,
    0xbf36000000000000,
    0xbf40000000000000,
    0xbf47000000000000,
    0xbf50000000000000,
    0xbf55800000000000,
    0xbf5c000000000000,
    0xbf61c00000000000,
    0xbf66000000000000,
    0xbf6ac00000000000,
    0xbf70000000000000,
    0x3f55000000000000,
    0x3f42000000000000,
    0xbf30000000000000,
    0xbf52000000000000,
    0xbf60800000000000,
    0xbf68800000000000,
    0x3f60c00000000000,
    0x3f50000000000000,
    0xbf24000000000000,
    0xbf56000000000000,
    0xbf65400000000000,
    0xbf70000000000000,
    0x3f50000000000000,
    0xbf3c000000000000,
    0xbf5f000000000000,
    0xbf6c000000000000,
    0x3f56800000000000,
    0xbf30000000000000,
    0xbf5f800000000000,
    0xbf6e000000000000,
    0x3f51000000000000,
    0xbf4a000000000000,
    0xbf66000000000000,
    0x3f60000000000000,
    0xbf00000000000000,
    0xbf61000000000000,
    0x3f64800000000000,
    0x3f38000000000000,
    0xbf5e000000000000,
    0x3f66000000000000,
    0x3f3a000000000000,
    0xbf60000000000000,
    0x3f64800000000000,
    0x3f10000000000000,
    0xbf64000000000000,
    0x3f60000000000000,
    0xbf45000000000000,
    0xbf6b000000000000,
    0x3f51000000000000,
    0xbf5c000000000000,
    0x3f65400000000000,
    0xbf30000000000000,
    0xbf69c00000000000,
    0x3f52000000000000,
    0xbf5f000000000000,
    0x3f63000000000000,
    0xbf49000000000000,
    0xbf70000000000000,
    0x3f30000000000000,
    0xbf68800000000000,
    0x3f52800000000000,
    0xbf62000000000000,
    0x3f5f000000000000,
    0xbf59000000000000,
    0x3f64c00000000000,
    0xbf50000000000000,
    0x3f69000000000000,
    0xbf42000000000000,
    0x3f6c400000000000,
    0xbf30000000000000,
    0x3f6e800000000000,
    0xbf10000000000000,
    0xbf70000000000000,
];

pub(crate) const LOG_COEFFS: [u64; 6] = [
    0xbfdfffffffffffff,
    0x3fd5555555554a9b,
    0xbfd0000000094567,
    0x3fc99999dcc9823c,
    0xbfc55550ac2e537a,
    0x3fc21a02c4e624d7,
];

/// Log2(x)
///
/// Max found ULP 0.5
pub fn f_log2(x: f64) -> f64 {
    let mut x_u = x.to_bits();

    const E_BIAS: u64 = (1u64 << (11 - 1u64)) - 1u64;

    let mut x_e: i64 = -(E_BIAS as i64);

    const MAX_NORMAL: u64 = f64::to_bits(f64::MAX);

    if x_u == 1f64.to_bits() {
        // log2(1.0) = +0.0
        return 0.0;
    }
    if x_u < min_normal_f64().to_bits() || x_u > MAX_NORMAL {
        if x == 0.0 {
            return f64::NEG_INFINITY;
        }
        if x < 0. || x.is_nan() {
            return f64::NAN;
        }
        if x.is_infinite() || x.is_nan() {
            return x + x;
        }
        // Normalize denormal inputs.
        x_u = (x * f64::from_bits(0x4330000000000000)).to_bits();
        x_e -= 52;
    }

    // log2(x) = log2(2^x_e * x_m)
    //         = x_e + log2(x_m)
    // Range reduction for log2(x_m):
    // For each x_m, we would like to find r such that:
    //   -2^-8 <= r * x_m - 1 < 2^-7
    let shifted = (x_u >> 45) as i64;
    let index = shifted & 0x7F;

    // Add unbiased exponent. Add an extra 1 if the 8 leading fractional bits are
    // all 1's.
    x_e = x_e.wrapping_add(x_u.wrapping_add(1u64 << 45).wrapping_shr(52) as i64);
    let e_x = x_e as f64;

    // Set m = 1.mantissa.
    let x_m = (x_u & 0x000F_FFFF_FFFF_FFFFu64) | 0x3FF0_0000_0000_0000u64;
    let m = f64::from_bits(x_m);

    let u;

    let r = f64::from_bits(LOG_RANGE_REDUCTION[index as usize]);

    #[cfg(any(
        all(
            any(target_arch = "x86", target_arch = "x86_64"),
            target_feature = "fma"
        ),
        target_arch = "aarch64"
    ))]
    {
        u = f_fmla(r, m, -1.0); // exact   
    }
    #[cfg(not(any(
        all(
            any(target_arch = "x86", target_arch = "x86_64"),
            target_feature = "fma"
        ),
        target_arch = "aarch64"
    )))]
    {
        let c_m = x_m & 0x3FFF_E000_0000_0000u64;
        let c = f64::from_bits(c_m);
        u = f_fmla(r, m - c, f64::from_bits(LOG_CD[index as usize])); // exact
    }

    // Exact sum:
    let log_vals = DoubleDouble::from_bit_pair(LOG2_DD[index as usize]);

    /*
       Poly generated in Sollya;
       d = [-2^-8, 2^-7];
       f = (log2(1 + x))/x;
       pf = fpminimax(f, 5, [|D...|], [-2^-8, 2^-7]);

       See ./notes/log2.sollya
    */
    const C: [u64; 6] = [
        0x3ff71547652b82fe,
        0xbfe71547652b7a07,
        0x3fdec709dc458db1,
        0xbfd715479c2266c9,
        0x3fd2776ae1ddf8f0,
        0xbfce7b2178870157,
    ];

    // Degree-7 minimax polynomial
    let p = f_fmla(
        f_polyeval6(
            u,
            f64::from_bits(C[0]),
            f64::from_bits(C[1]),
            f64::from_bits(C[2]),
            f64::from_bits(C[3]),
            f64::from_bits(C[4]),
            f64::from_bits(C[5]),
        ),
        u,
        log_vals.lo,
    );

    // log2(x) = e + log2(r) + log2(index)

    let rr = DoubleDouble::from_full_exact_add(log_vals.hi, e_x);
    let r3 = DoubleDouble::add_f64(rr, p);

    let err = f_fmla(
        r3.hi,
        f64::from_bits(0x3bf0000000000000), // 2^-64
        f64::from_bits(0x3c60000000000000), // 2^-57
    );

    // Lower bound from the result
    let left = r3.hi + (r3.lo - err);
    // Upper bound from the result
    let right = r3.hi + (r3.lo + err);

    // Ziv's test if fast pass is accurate enough.
    if left == right {
        return left;
    }

    log2_hard(x)
}

#[cold]
#[inline(never)]
fn log2_hard(x: f64) -> f64 {
    let r = log2_dd(x);
    let err = f_fmla(
        r.hi,
        f64::from_bits(0x3b50000000000000), // 2^-74
        f64::from_bits(0x3990000000000000), // 2^-102
    );
    let ub = r.hi + (r.lo + err);
    let lb = r.hi + (r.lo - err);
    if ub == lb {
        return r.to_f64();
    }
    log2_hard_slow(x)
}

#[cold]
#[inline(never)]
fn log2_hard_slow(x: f64) -> f64 {
    log2_td(x).to_f64()
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_log2d() {
        assert_eq!(f_log2(0.7500000000000461), -0.4150374992787552);
        assert_eq!(f_log2(1.99999999779061), 0.999999998406262);
        assert_eq!(f_log2(0.7812499981882864), -0.35614381357087554);
        assert_eq!(f_log2(0.5937499997089616), -0.7520724872635803);
        assert_eq!(f_log2(0.9983015060407214), -0.0024524921736992287);
        assert_eq!(f_log2(0.2284926469437778), -2.129780355811612);
        assert_eq!(f_log2(24.), 4.584962500721156181453738943);
        assert!((f_log2(0.35) - 0.35f64.log2()).abs() < 1e-8);
        assert!((f_log2(0.9) - 0.9f64.log2()).abs() < 1e-8);
        assert_eq!(f_log2(0.), f64::NEG_INFINITY);
        assert!(f_log2(-1.).is_nan());
        assert!(f_log2(f64::NAN).is_nan());
        assert!(f_log2(f64::NEG_INFINITY).is_nan());
        assert_eq!(f_log2(f64::INFINITY), f64::INFINITY);
    }

    #[test]
    fn test_log2_control_values() {
        assert_eq!(
            f_log2(f64::from_bits(0x3ff1406d79e1b574)),
            0.10866415915666149
        );
        assert_eq!(
            f_log2(f64::from_bits(0x3ffb4ebe40c95a01)),
            0.7712301192941611
        );
        assert_eq!(
            f_log2(f64::from_bits(0x3ff0b541b6746bd1)),
            0.062470074690080965
        );
        assert_eq!(
            f_log2(f64::from_bits(0x3ff68d778f076021)),
            0.4952222169409832
        );
        assert_eq!(
            f_log2(f64::from_bits(0x3ff67eaf07ce24d1)),
            0.4915233709893074
        );
        assert_eq!(
            f_log2(f64::from_bits(0x3ff0b53197bd66c8)),
            0.062448835548542664
        );
        assert_eq!(
            f_log2(f64::from_bits(0x3fe6b015f8d9a784)),
            -0.496152942566177
        );
        assert_eq!(
            f_log2(f64::from_bits(0x3ff160732376ad7f)),
            0.11908694357502585
        );
        assert_eq!(
            f_log2(f64::from_bits(0x3ff0b53f741fb8c4)),
            0.062467098188572705
        );
    }
}

Coverage Report

Created: 2025-10-14 06:57

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::log2dd::log2_dd;
32		use crate::logs::log2td::log2_td;
33		use crate::polyeval::f_polyeval6;
34
35		pub(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		/**
167		Generated by SageMath:
168		```python
169		values = LOG_RANGE_REDUCTION
170
171		R = RealField(150)
172
173		def hex_to_float(h):
174		return struct.unpack('>d', struct.pack('>Q', h))[0]
175
176		real_array = [R(hex_to_float(h)) for h in values]
177
178		for r in real_array:
179		print_double_double("", -RealField(180)(r).log2())
180		```
181		**/
182		static 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		target_arch = "aarch64"
319		)))]
320		pub(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
451		pub(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
463	0	pub fn f_log2(x: f64) -> f64 {
464	0	let mut x_u = x.to_bits();
465
466		const E_BIAS: u64 = (1u64 << (11 - 1u64)) - 1u64;
467
468	0	let mut x_e: i64 = -(E_BIAS as i64);
469
470		const MAX_NORMAL: u64 = f64::to_bits(f64::MAX);
471
472	0	if x_u == 1f64.to_bits() {
473		// log2(1.0) = +0.0
474	0	return 0.0;
475	0	}
476	0	if x_u < min_normal_f64().to_bits() \|\| x_u > MAX_NORMAL {
477	0	if x == 0.0 {
478	0	return f64::NEG_INFINITY;
479	0	}
480	0	if x < 0. \|\| x.is_nan() {
481	0	return f64::NAN;
482	0	}
483	0	if x.is_infinite() \|\| x.is_nan() {
484	0	return x + x;
485	0	}
486		// Normalize denormal inputs.
487	0	x_u = (x * f64::from_bits(0x4330000000000000)).to_bits();
488	0	x_e -= 52;
489	0	}
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	0	let shifted = (x_u >> 45) as i64;
497	0	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	0	x_e = x_e.wrapping_add(x_u.wrapping_add(1u64 << 45).wrapping_shr(52) as i64);
502	0	let e_x = x_e as f64;
503
504		// Set m = 1.mantissa.
505	0	let x_m = (x_u & 0x000F_FFFF_FFFF_FFFFu64) \| 0x3FF0_0000_0000_0000u64;
506	0	let m = f64::from_bits(x_m);
507
508		let u;
509
510	0	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		target_arch = "aarch64"
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		target_arch = "aarch64"
528		)))]
529	0	{
530	0	let c_m = x_m & 0x3FFF_E000_0000_0000u64;
531	0	let c = f64::from_bits(c_m);
532	0	u = f_fmla(r, m - c, f64::from_bits(LOG_CD[index as usize])); // exact
533	0	}
534
535		// Exact sum:
536	0	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	0	let p = f_fmla(
557	0	f_polyeval6(
558	0	u,
559	0	f64::from_bits(C[0]),
560	0	f64::from_bits(C[1]),
561	0	f64::from_bits(C[2]),
562	0	f64::from_bits(C[3]),
563	0	f64::from_bits(C[4]),
564	0	f64::from_bits(C[5]),
565		),
566	0	u,
567	0	log_vals.lo,
568		);
569
570		// log2(x) = e + log2(r) + log2(index)
571
572	0	let rr = DoubleDouble::from_full_exact_add(log_vals.hi, e_x);
573	0	let r3 = DoubleDouble::add_f64(rr, p);
574
575	0	let err = f_fmla(
576	0	r3.hi,
577	0	f64::from_bits(0x3bf0000000000000), // 2^-64
578	0	f64::from_bits(0x3c60000000000000), // 2^-57
579		);
580
581		// Lower bound from the result
582	0	let left = r3.hi + (r3.lo - err);
583		// Upper bound from the result
584	0	let right = r3.hi + (r3.lo + err);
585
586		// Ziv's test if fast pass is accurate enough.
587	0	if left == right {
588	0	return left;
589	0	}
590
591	0	log2_hard(x)
592	0	}
593
594		#[cold]
595		#[inline(never)]
596	0	fn log2_hard(x: f64) -> f64 {
597	0	let r = log2_dd(x);
598	0	let err = f_fmla(
599	0	r.hi,
600	0	f64::from_bits(0x3b50000000000000), // 2^-74
601	0	f64::from_bits(0x3990000000000000), // 2^-102
602		);
603	0	let ub = r.hi + (r.lo + err);
604	0	let lb = r.hi + (r.lo - err);
605	0	if ub == lb {
606	0	return r.to_f64();
607	0	}
608	0	log2_hard_slow(x)
609	0	}
610
611		#[cold]
612		#[inline(never)]
613	0	fn log2_hard_slow(x: f64) -> f64 {
614	0	log2_td(x).to_f64()
615	0	}
616
617		#[cfg(test)]
618		mod 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		}