/rust/registry/src/index.crates.io-1949cf8c6b5b557f/pxfm-0.1.27/src/tangent/atan2.rs

Source
/*
 * // Copyright (c) Radzivon Bartoshyk 6/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;
use crate::double_double::DoubleDouble;
use crate::dyadic_float::{DyadicFloat128, DyadicSign};
use crate::polyeval::f_polyeval4;
use crate::rounding::CpuRound;

// atan(i/64) with i = 0..64, generated by Sollya with:
// > for i from 0 to 64 do {
//     a = round(atan(i/64), D, RN);
//     b = round(atan(i/64) - a, D, RN);
//     print("{", b, ",", a, "},");
//   };
pub(crate) static ATAN_I: [(u64, u64); 65] = [
    (0x0000000000000000, 0x0000000000000000),
    (0xbc2220c39d4dff50, 0x3f8fff555bbb729b),
    (0xbc35ec431444912c, 0x3f9ffd55bba97625),
    (0xbc086ef8f794f105, 0x3fa7fb818430da2a),
    (0xbc3c934d86d23f1d, 0x3faff55bb72cfdea),
    (0x3c5ac4ce285df847, 0x3fb3f59f0e7c559d),
    (0xbc5cfb654c0c3d98, 0x3fb7ee182602f10f),
    (0x3c5f7b8f29a05987, 0x3fbbe39ebe6f07c3),
    (0xbc4cd37686760c17, 0x3fbfd5ba9aac2f6e),
    (0xbc4b485914dacf8c, 0x3fc1e1fafb043727),
    (0x3c661a3b0ce9281b, 0x3fc3d6eee8c6626c),
    (0xbc5054ab2c010f3d, 0x3fc5c9811e3ec26a),
    (0x3c5347b0b4f881ca, 0x3fc7b97b4bce5b02),
    (0x3c4cf601e7b4348e, 0x3fc9a6a8e96c8626),
    (0x3c217b10d2e0e5ab, 0x3fcb90d7529260a2),
    (0x3c6c648d1534597e, 0x3fcd77d5df205736),
    (0x3c68ab6e3cf7afbd, 0x3fcf5b75f92c80dd),
    (0x3c762e47390cb865, 0x3fd09dc597d86362),
    (0x3c630ca4748b1bf9, 0x3fd18bf5a30bf178),
    (0xbc7077cdd36dfc81, 0x3fd278372057ef46),
    (0xbc6963a544b672d8, 0x3fd362773707ebcc),
    (0xbc75d5e43c55b3ba, 0x3fd44aa436c2af0a),
    (0xbc62566480884082, 0x3fd530ad9951cd4a),
    (0xbc7a725715711f00, 0x3fd614840309cfe2),
    (0xbc7c63aae6f6e918, 0x3fd6f61941e4def1),
    (0x3c769c885c2b249a, 0x3fd7d5604b63b3f7),
    (0x3c7b6d0ba3748fa8, 0x3fd8b24d394a1b25),
    (0x3c79e6c988fd0a77, 0x3fd98cd5454d6b18),
    (0xbc724dec1b50b7ff, 0x3fda64eec3cc23fd),
    (0x3c7ae187b1ca5040, 0x3fdb3a911da65c6c),
    (0xbc7cc1ce70934c34, 0x3fdc0db4c94ec9f0),
    (0xbc7a2cfa4418f1ad, 0x3fdcde53432c1351),
    (0x3c7a2b7f222f65e2, 0x3fddac670561bb4f),
    (0x3c70e53dc1bf3435, 0x3fde77eb7f175a34),
    (0xbc6a3992dc382a23, 0x3fdf40dd0b541418),
    (0xbc8b32c949c9d593, 0x3fe0039c73c1a40c),
    (0xbc7d5b495f6349e6, 0x3fe0657e94db30d0),
    (0x3c5974fa13b5404f, 0x3fe0c6145b5b43da),
    (0xbc52bdaee1c0ee35, 0x3fe1255d9bfbd2a9),
    (0x3c8c621cec00c301, 0x3fe1835a88be7c13),
    (0xbc5928df287a668f, 0x3fe1e00babdefeb4),
    (0x3c6c421c9f38224e, 0x3fe23b71e2cc9e6a),
    (0xbc709e73b0c6c087, 0x3fe2958e59308e31),
    (0x3c8c5d5e9ff0cf8d, 0x3fe2ee628406cbca),
    (0x3c81021137c71102, 0x3fe345f01cce37bb),
    (0xbc82304331d8bf46, 0x3fe39c391cd4171a),
    (0x3c7ecf8b492644f0, 0x3fe3f13fb89e96f4),
    (0xbc7f76d0163f79c8, 0x3fe445065b795b56),
    (0x3c72419a87f2a458, 0x3fe4978fa3269ee1),
    (0x3c84a33dbeb3796c, 0x3fe4e8de5bb6ec04),
    (0xbc81bb74abda520c, 0x3fe538f57b89061f),
    (0xbc75e5c9d8c5a950, 0x3fe587d81f732fbb),
    (0x3c60028e4bc5e7ca, 0x3fe5d58987169b18),
    (0xbc62b785350ee8c1, 0x3fe6220d115d7b8e),
    (0xbc76ea6febe8bbba, 0x3fe66d663923e087),
    (0xbc8a80386188c50e, 0x3fe6b798920b3d99),
    (0xbc78c34d25aadef6, 0x3fe700a7c5784634),
    (0x3c47b2a6165884a1, 0x3fe748978fba8e0f),
    (0x3c8406a089803740, 0x3fe78f6bbd5d315e),
    (0x3c8560821e2f3aa9, 0x3fe7d528289fa093),
    (0xbc7bf76229d3b917, 0x3fe819d0b7158a4d),
    (0x3c66b66e7fc8b8c3, 0x3fe85d69576cc2c5),
    (0xbc855b9a5e177a1b, 0x3fe89ff5ff57f1f8),
    (0xbc7ec182ab042f61, 0x3fe8e17aa99cc05e),
    (0x3c81a62633145c07, 0x3fe921fb54442d18),
];

// Sage math:
// print("pub(crate) static ATAN_RATIONAL_128: [DyadicFloat128; 65] = [")
// for i in range(65):
//     x = D3(i)/D3(64)
//     v = atan(x)
//     print_dyadic(v)
//
// print("];")
pub(crate) static ATAN_RATIONAL_128: [DyadicFloat128; 65] = [
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: 0,
        mantissa: 0x0_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -134,
        mantissa: 0xfffaaadd_db94d5bb_e78c5640_15f76048_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -133,
        mantissa: 0xffeaaddd_4bb12542_779d776d_da8c6214_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -132,
        mantissa: 0xbfdc0c21_86d14fcf_220e10d6_1df56ec7_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -132,
        mantissa: 0xffaaddb9_67ef4e36_cb2792dc_0e2e0d51_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -131,
        mantissa: 0x9facf873_e2aceb58_99c50bbf_08e6cdf6_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -131,
        mantissa: 0xbf70c130_17887460_93567e78_4cf83676_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -131,
        mantissa: 0xdf1cf5f3_783e1bef_71e5340b_30e5d9ef_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -131,
        mantissa: 0xfeadd4d5_617b6e32_c897989f_3e888ef8_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -130,
        mantissa: 0x8f0fd7d8_21b93725_bd375929_83a0af9a_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -130,
        mantissa: 0x9eb77746_331362c3_47619d25_0360fe85_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -130,
        mantissa: 0xae4c08f1_f6134efa_b54d3fef_0c2de994_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -130,
        mantissa: 0xbdcbda5e_72d81134_7b0b4f88_1c9c7488_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -130,
        mantissa: 0xcd35474b_643130e7_b00f3da1_a46eeb3b_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -130,
        mantissa: 0xdc86ba94_93051022_f621a5c1_cb552f03_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -130,
        mantissa: 0xebbeaef9_02b9b38c_91a2a68b_2fbd78e8_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -130,
        mantissa: 0xfadbafc9_6406eb15_6dc79ef5_f7a217e6_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -129,
        mantissa: 0x84ee2cbe_c31b12c5_c8e72197_0cabd3a3_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -129,
        mantissa: 0x8c5fad18_5f8bc130_ca4748b1_bf88298d_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -129,
        mantissa: 0x93c1b902_bf7a2df1_06459240_6fe1447a_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -129,
        mantissa: 0x9b13b9b8_3f5e5e69_c5abb498_d27af328_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -129,
        mantissa: 0xa25521b6_15784d45_43787549_88b8d9e3_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -129,
        mantissa: 0xa9856cca_8e6a4eda_99b7f77b_f7d9e8c1_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -129,
        mantissa: 0xb0a42018_4e7f0cb1_b51d51dc_200a0fc3_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -129,
        mantissa: 0xb7b0ca0f_26f78473_8aa32122_dcfe4483_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -129,
        mantissa: 0xbeab025b_1d9fbad3_910b8564_93411026_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -129,
        mantissa: 0xc59269ca_50d92b6d_a1746e91_f50a28de_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -129,
        mantissa: 0xcc66aa2a_6b58c33c_d9311fa1_4ed9b7c4_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -129,
        mantissa: 0xd327761e_611fe5b6_427c95e9_001e7136_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -129,
        mantissa: 0xd9d488ed_32e3635c_30f6394a_0806345d_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -129,
        mantissa: 0xe06da64a_764f7c67_c631ed96_798cb804_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -129,
        mantissa: 0xe6f29a19_609a84ba_60b77ce1_ca6dc2c8_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -129,
        mantissa: 0xed63382b_0dda7b45_6fe445ec_bc3a8d03_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -129,
        mantissa: 0xf3bf5bf8_bad1a21c_a7b837e6_86adf3fa_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -129,
        mantissa: 0xfa06e85a_a0a0be5c_66d23c7d_5dc8ecc2_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -128,
        mantissa: 0x801ce39e_0d205c99_a6d6c6c5_4d938596_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -128,
        mantissa: 0x832bf4a6_d9867e2a_4b6a09cb_61a515c1_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -128,
        mantissa: 0x8630a2da_da1ed065_d3e84ed5_013ca37e_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -128,
        mantissa: 0x892aecdf_de9547b5_094478fc_472b4afc_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -128,
        mantissa: 0x8c1ad445_f3e09b8c_439d8018_60205921_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -128,
        mantissa: 0x8f005d5e_f7f59f9b_5c835e16_65c43748_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -128,
        mantissa: 0x91db8f16_64f350e2_10e4f9c1_126e0220_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -128,
        mantissa: 0x94ac72c9_847186f6_18c4f393_f78a32f9_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -128,
        mantissa: 0x97731420_365e538b_abd3fe19_f1aeb6b3_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -128,
        mantissa: 0x9a2f80e6_71bdda20_4226f8e2_204ff3bd_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -128,
        mantissa: 0x9ce1c8e6_a0b8cdb9_f799c4e8_174cf11c_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -128,
        mantissa: 0x9f89fdc4_f4b7a1ec_f8b49264_4f0701e0_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -128,
        mantissa: 0xa22832db_cadaae08_92fe9c08_637af0e6_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -128,
        mantissa: 0xa4bc7d19_34f70924_19a87f2a_457dac9f_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -128,
        mantissa: 0xa746f2dd_b7602294_67b7d66f_2d74e019_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -128,
        mantissa: 0xa9c7abdc_4830f5c8_916a84b5_be7933f6_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -128,
        mantissa: 0xac3ec0fb_997dd6a1_a36273a5_6afa8ef4_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -128,
        mantissa: 0xaeac4c38_b4d8c080_14725e2f_3e52070a_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -128,
        mantissa: 0xb110688a_ebdc6f6a_43d65788_b9f6a7b5_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -128,
        mantissa: 0xb36b31c9_1f043691_59014174_4462f93a_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -128,
        mantissa: 0xb5bcc490_59ecc4af_f8f3cee7_5e3907d5_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -128,
        mantissa: 0xb8053e2b_c2319e73_cb2da552_10a4443d_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -128,
        mantissa: 0xba44bc7d_d470782f_654c2cb1_0942e386_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -128,
        mantissa: 0xbc7b5dea_e98af280_d4113006_e80fb290_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -128,
        mantissa: 0xbea94144_fd049aac_1043c5e7_55282e7d_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -128,
        mantissa: 0xc0ce85b8_ac526640_89dd62c4_6e92fa25_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -128,
        mantissa: 0xc2eb4abb_661628b5_b373fe45_c61bb9fb_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -128,
        mantissa: 0xc4ffaffa_bf8fbd54_8cb43d10_bc9e0221_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -128,
        mantissa: 0xc70bd54c_e602ee13_e7d54fbd_09f2be38_u128,
    },
    DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -128,
        mantissa: 0xc90fdaa2_2168c234_c4c6628b_80dc1cd1_u128,
    },
];

// Approximate atan(x) for |x| <= 2^-7.
// Using degree-9 Taylor polynomial:
//  P = x - x^3/3 + x^5/5 -x^7/7 + x^9/9;
// Then the absolute error is bounded by:
//   |atan(x) - P(x)| < |x|^11/11 < 2^(-7*11) / 11 < 2^-80.
// And the relative error is bounded by:
//   |(atan(x) - P(x))/atan(x)| < |x|^10 / 10 < 2^-73.
// For x = x_hi + x_lo, fully expand the polynomial and drop any terms less than
//   ulp(x_hi^3 / 3) gives us:
// P(x) ~ x_hi - x_hi^3/3 + x_hi^5/5 - x_hi^7/7 + x_hi^9/9 +
//        + x_lo * (1 - x_hi^2 + x_hi^4)
// Since p.lo is ~ x^3/3, the relative error from rounding is bounded by:
//   |(atan(x) - P(x))/atan(x)| < ulp(x^2) <= 2^(-14-52) = 2^-66.
#[inline]
pub(crate) fn atan_eval(x: DoubleDouble) -> DoubleDouble {
    let p_hi = x.hi;
    let x_hi_sq = x.hi * x.hi;
    // c0 ~ x_hi^2 * 1/5 - 1/3
    let c0 = f_fmla(
        x_hi_sq,
        f64::from_bits(0x3fc999999999999a),
        f64::from_bits(0xbfd5555555555555),
    );
    // c1 ~ x_hi^2 * 1/9 - 1/7
    let c1 = f_fmla(
        x_hi_sq,
        f64::from_bits(0x3fbc71c71c71c71c),
        f64::from_bits(0xbfc2492492492492),
    );
    // x_hi^3
    let x_hi_3 = x_hi_sq * x.hi;
    // x_hi^4
    let x_hi_4 = x_hi_sq * x_hi_sq;
    // d0 ~ 1/3 - x_hi^2 / 5 + x_hi^4 / 7 - x_hi^6 / 9
    let d0 = f_fmla(x_hi_4, c1, c0);
    // x_lo - x_lo * x_hi^2 + x_lo * x_hi^4
    let d1 = f_fmla(x_hi_4 - x_hi_sq, x.lo, x.lo);
    // p.lo ~ -x_hi^3/3 + x_hi^5/5 - x_hi^7/7 + x_hi^9/9 +
    //        + x_lo * (1 - x_hi^2 + x_hi^4)
    let p_lo = f_fmla(x_hi_3, d0, d1);
    DoubleDouble::new(p_lo, p_hi)
}

#[inline(always)]
#[allow(unused)]
pub(crate) fn atan_eval_fma(x: DoubleDouble) -> DoubleDouble {
    let p_hi = x.hi;
    let x_hi_sq = x.hi * x.hi;
    // c0 ~ x_hi^2 * 1/5 - 1/3
    let c0 = f64::mul_add(
        x_hi_sq,
        f64::from_bits(0x3fc999999999999a),
        f64::from_bits(0xbfd5555555555555),
    );
    // c1 ~ x_hi^2 * 1/9 - 1/7
    let c1 = f64::mul_add(
        x_hi_sq,
        f64::from_bits(0x3fbc71c71c71c71c),
        f64::from_bits(0xbfc2492492492492),
    );
    // x_hi^3
    let x_hi_3 = x_hi_sq * x.hi;
    // x_hi^4
    let x_hi_4 = x_hi_sq * x_hi_sq;
    // d0 ~ 1/3 - x_hi^2 / 5 + x_hi^4 / 7 - x_hi^6 / 9
    let d0 = f64::mul_add(x_hi_4, c1, c0);
    // x_lo - x_lo * x_hi^2 + x_lo * x_hi^4
    let d1 = f64::mul_add(x_hi_4 - x_hi_sq, x.lo, x.lo);
    // p.lo ~ -x_hi^3/3 + x_hi^5/5 - x_hi^7/7 + x_hi^9/9 +
    //        + x_lo * (1 - x_hi^2 + x_hi^4)
    let p_lo = f64::mul_add(x_hi_3, d0, d1);
    DoubleDouble::new(p_lo, p_hi)
}

#[inline]
fn atan_eval_hard(x: DyadicFloat128) -> DyadicFloat128 {
    // let x_hi_sq = x * x;
    // c0 ~ x_hi^2 * 1/5 - 1/3
    const C: [DyadicFloat128; 4] = [
        DyadicFloat128 {
            sign: DyadicSign::Neg,
            exponent: -129,
            mantissa: 0xaaaaaaaa_aaaaaaaa_aaaaaaaa_aaaaaaab_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -130,
            mantissa: 0xcccccccc_cccccccc_cccccccc_cccccccd_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Neg,
            exponent: -130,
            mantissa: 0x92492492_49249249_24924924_92492492_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -131,
            mantissa: 0xe38e38e3_8e38e38e_38e38e38_e38e38e4_u128,
        },
    ];
    let dx2 = x * x;
    // Taylor polynomial
    // P = x - x^3/3 + x^5/5 -x^7/7 + x^9/9;
    let p = f_polyeval4(dx2, C[0], C[1], C[2], C[3]);
    x + dx2 * x * p
}

#[cold]
#[inline(never)]
pub(crate) fn atan2_hard(y: f64, x: f64) -> DyadicFloat128 {
    /*
       Sage math:
       from sage.all import *

       def format_dyadic_hex(value):
           l = hex(value)[2:]
           n = 8
           x = [l[i:i + n] for i in range(0, len(l), n)]
           return "0x" + "_".join(x) + "_u128"

       def print_dyadic(value):
           (s, m, e) = RealField(128)(value).sign_mantissa_exponent();
           print("DyadicFloat128 {")
           print(f"    sign: DyadicSign::{'Pos' if s >= 0 else 'Neg'},")
           print(f"    exponent: {e},")
           print(f"    mantissa: {format_dyadic_hex(m)},")
           print("},")

       D3 = RealField(157)

       print("Minus Pi")
       print_dyadic(-D3.pi())

       print("\nPI over 2")
       print_dyadic(D3.pi() / 2)

       print("\nMinus PI over 2")
       print_dyadic(-D3.pi() / 2)
    */
    static ZERO: DyadicFloat128 = DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: 0,
        mantissa: 0_u128,
    };
    static MINUS_PI: DyadicFloat128 = DyadicFloat128 {
        sign: DyadicSign::Neg,
        exponent: -126,
        mantissa: 0xc90fdaa2_2168c234_c4c6628b_80dc1cd1_u128,
    };
    static PI_OVER_2: DyadicFloat128 = DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -127,
        mantissa: 0xc90fdaa2_2168c234_c4c6628b_80dc1cd1_u128,
    };
    static MPI_OVER_2: DyadicFloat128 = DyadicFloat128 {
        sign: DyadicSign::Neg,
        exponent: -127,
        mantissa: 0xc90fdaa2_2168c234_c4c6628b_80dc1cd1_u128,
    };
    static CONST_ADJ: [[[DyadicFloat128; 2]; 2]; 2] = [
        [[ZERO, MPI_OVER_2], [ZERO, MPI_OVER_2]],
        [[MINUS_PI, PI_OVER_2], [MINUS_PI, PI_OVER_2]],
    ];

    let x_sign = x.is_sign_negative() as usize;
    let y_sign = y.is_sign_negative() as usize;
    let x_bits = x.to_bits() & 0x7fff_ffff_ffff_ffff;
    let y_bits = y.to_bits() & 0x7fff_ffff_ffff_ffff;
    let x_abs = x_bits;
    let y_abs = y_bits;
    let recip = x_abs < y_abs;

    let min_abs = if recip { x_abs } else { y_abs };
    let max_abs = if !recip { x_abs } else { y_abs };
    let min_exp = min_abs.wrapping_shr(52);
    let max_exp = max_abs.wrapping_shr(52);

    let mut num = f64::from_bits(min_abs);
    let mut den = f64::from_bits(max_abs);

    if max_exp > 0x7ffu64 - 128u64 || min_exp < 128u64 {
        let scale_up = min_exp < 128u64;
        let scale_down = max_exp > 0x7ffu64 - 128u64;
        // At least one input is denormal, multiply both numerator and denominator
        // by some large enough power of 2 to normalize denormal inputs.
        if scale_up {
            num *= f64::from_bits(0x43f0000000000000);
            if !scale_down {
                den *= f64::from_bits(0x43f0000000000000)
            }
        } else if scale_down {
            den *= f64::from_bits(0x3bf0000000000000);
            if !scale_up {
                num *= f64::from_bits(0x3bf0000000000000);
            }
        }
    }

    static IS_NEG: [DyadicSign; 2] = [DyadicSign::Pos, DyadicSign::Neg];

    let final_sign = IS_NEG[((x_sign != y_sign) != recip) as usize];
    let const_term = CONST_ADJ[x_sign][y_sign][recip as usize];

    let num = DyadicFloat128::new_from_f64(num);
    let den = DyadicFloat128::new_from_f64(den);

    let den_recip0 = den.reciprocal();

    let mut k_product = num * den_recip0;
    k_product.exponent += 6;

    let mut k = k_product.round_to_nearest_f64();
    let idx = k as u64;
    // k = idx / 64
    k *= f64::from_bits(0x3f90000000000000);

    // Range reduction:
    // atan(n/d) - atan(k/64) = atan((n/d - k/64) / (1 + (n/d) * (k/64)))
    //                        = atan((n - d * k/64)) / (d + n * k/64))
    let k_rational128 = DyadicFloat128::new_from_f64(k);
    let num_k = num * k_rational128;
    let den_k = den * k_rational128;

    // num_dd = n - d * k
    let num_rational128 = num - den_k;
    // den_dd = d + n * k
    let den_rational128 = den + num_k;

    // q = (n - d * k) / (d + n * k)
    let den_rational128_recip = den_rational128.reciprocal();
    let q = num_rational128 * den_rational128_recip;

    let p = atan_eval_hard(q);

    let vl = ATAN_RATIONAL_128[idx as usize];
    let mut r = p + vl + const_term;
    r.sign = r.sign.mult(final_sign);

    r
}

static IS_NEG: [f64; 2] = [1.0, -1.0];
const ZERO: DoubleDouble = DoubleDouble::new(0.0, 0.0);
const MZERO: DoubleDouble = DoubleDouble::new(-0.0, -0.0);
const PI: DoubleDouble = DoubleDouble::new(
    f64::from_bits(0x3ca1a62633145c07),
    f64::from_bits(0x400921fb54442d18),
);
const MPI: DoubleDouble = DoubleDouble::new(
    f64::from_bits(0xbca1a62633145c07),
    f64::from_bits(0xc00921fb54442d18),
);
const PI_OVER_2: DoubleDouble = DoubleDouble::new(
    f64::from_bits(0x3c91a62633145c07),
    f64::from_bits(0x3ff921fb54442d18),
);
const MPI_OVER_2: DoubleDouble = DoubleDouble::new(
    f64::from_bits(0xbc91a62633145c07),
    f64::from_bits(0xbff921fb54442d18),
);
const PI_OVER_4: DoubleDouble = DoubleDouble::new(
    f64::from_bits(0x3c81a62633145c07),
    f64::from_bits(0x3fe921fb54442d18),
);
const THREE_PI_OVER_4: DoubleDouble = DoubleDouble::new(
    f64::from_bits(0x3c9a79394c9e8a0a),
    f64::from_bits(0x4002d97c7f3321d2),
);

// Adjustment for constant term:
//   CONST_ADJ[x_sign][y_sign][recip]
static CONST_ADJ: [[[DoubleDouble; 2]; 2]; 2] = [
    [[ZERO, MPI_OVER_2], [MZERO, MPI_OVER_2]],
    [[MPI, PI_OVER_2], [MPI, PI_OVER_2]],
];

// Exceptional cases:
//   EXCEPT[y_except][x_except][x_is_neg]
// with x_except & y_except:
//   0: zero
//   1: finite, non-zero
//   2: infinity
static EXCEPTS: [[[DoubleDouble; 2]; 3]; 3] = [
    [[ZERO, PI], [ZERO, PI], [ZERO, PI]],
    [[PI_OVER_2, PI_OVER_2], [ZERO, ZERO], [ZERO, PI]],
    [
        [PI_OVER_2, PI_OVER_2],
        [PI_OVER_2, PI_OVER_2],
        [PI_OVER_4, THREE_PI_OVER_4],
    ],
];

#[inline(always)]
fn atan2_gen_fma(y: f64, x: f64) -> f64 {
    let x_sign = x.is_sign_negative() as usize;
    let y_sign = y.is_sign_negative() as usize;
    let x_bits = x.to_bits() & 0x7fff_ffff_ffff_ffff;
    let y_bits = y.to_bits() & 0x7fff_ffff_ffff_ffff;
    let x_abs = x_bits;
    let y_abs = y_bits;
    let recip = x_abs < y_abs;
    let mut min_abs = if recip { x_abs } else { y_abs };
    let mut max_abs = if !recip { x_abs } else { y_abs };
    let mut min_exp = min_abs.wrapping_shr(52);
    let mut max_exp = max_abs.wrapping_shr(52);

    let mut num = f64::from_bits(min_abs);
    let mut den = f64::from_bits(max_abs);

    // Check for exceptional cases, whether inputs are 0, inf, nan, or close to
    // overflow, or close to underflow.
    if max_exp > 0x7ffu64 - 128u64 || min_exp < 128u64 {
        if x.is_nan() || y.is_nan() {
            return f64::NAN;
        }
        let x_except = if x == 0.0 {
            0
        } else if x.is_infinite() {
            2
        } else {
            1
        };
        let y_except = if y == 0.0 {
            0
        } else if y.is_infinite() {
            2
        } else {
            1
        };

        if (x_except != 1) || (y_except != 1) {
            let r = EXCEPTS[y_except][x_except][x_sign];
            return f_fmla(IS_NEG[y_sign], r.hi, IS_NEG[y_sign] * r.lo);
        }
        let scale_up = min_exp < 128u64;
        let scale_down = max_exp > 0x7ffu64 - 128u64;
        // At least one input is denormal, multiply both numerator and denominator
        // by some large enough power of 2 to normalize denormal inputs.
        // if scale_up || scale_down {
        //     return atan2_hard(y, x).fast_as_f64();
        // }
        if scale_up {
            num *= f64::from_bits(0x43f0000000000000);
            if !scale_down {
                den *= f64::from_bits(0x43f0000000000000);
            }
        } else if scale_down {
            den *= f64::from_bits(0x3bf0000000000000);
            if !scale_up {
                num *= f64::from_bits(0x3bf0000000000000);
            }
        }

        min_abs = num.to_bits();
        max_abs = den.to_bits();
        min_exp = min_abs.wrapping_shr(52);
        max_exp = max_abs.wrapping_shr(52);
    }
    let final_sign = IS_NEG[((x_sign != y_sign) != recip) as usize];
    let const_term = CONST_ADJ[x_sign][y_sign][recip as usize];
    let exp_diff = max_exp - min_exp;
    // We have the following bound for normalized n and d:
    //   2^(-exp_diff - 1) < n/d < 2^(-exp_diff + 1).
    if exp_diff > 54 {
        return f_fmla(
            final_sign,
            const_term.hi,
            final_sign * (const_term.lo + num / den),
        );
    }

    let mut k = (64.0 * num / den).cpu_round();
    let idx = k as u64;
    // k = idx / 64
    k *= f64::from_bits(0x3f90000000000000);

    // Range reduction:
    // atan(n/d) - atan(k/64) = atan((n/d - k/64) / (1 + (n/d) * (k/64)))
    //                        = atan((n - d * k/64)) / (d + n * k/64))
    let num_k = DoubleDouble::from_exact_mult(num, k);
    let den_k = DoubleDouble::from_exact_mult(den, k);

    // num_dd = n - d * k
    let num_dd = DoubleDouble::from_exact_add(num - den_k.hi, -den_k.lo);
    // den_dd = d + n * k
    let mut den_dd = DoubleDouble::from_exact_add(den, num_k.hi);
    den_dd.lo += num_k.lo;

    // q = (n - d * k) / (d + n * k)
    let q = DoubleDouble::div(num_dd, den_dd);
    // p ~ atan(q)
    let p = atan_eval(q);

    let vl = ATAN_I[idx as usize];
    let vlo = DoubleDouble::from_bit_pair(vl);
    let mut r = DoubleDouble::add(const_term, DoubleDouble::add(vlo, p));

    let err = f_fmla(
        p.hi,
        f64::from_bits(0x3bd0000000000000),
        f64::from_bits(0x3c00000000000000),
    );

    let ub = r.hi + (r.lo + err);
    let lb = r.hi + (r.lo - err);

    if ub == lb {
        r.hi *= final_sign;
        r.lo *= final_sign;

        return r.to_f64();
    }
    atan2_hard(y, x).fast_as_f64()
}

#[cfg(any(target_arch = "x86", target_arch = "x86_64"))]
#[target_feature(enable = "avx", enable = "fma")]
unsafe fn atan2_fma_impl(y: f64, x: f64) -> f64 {
    let x_sign = x.is_sign_negative() as usize;
    let y_sign = y.is_sign_negative() as usize;
    let x_bits = x.to_bits() & 0x7fff_ffff_ffff_ffff;
    let y_bits = y.to_bits() & 0x7fff_ffff_ffff_ffff;
    let x_abs = x_bits;
    let y_abs = y_bits;
    let recip = x_abs < y_abs;
    let mut min_abs = if recip { x_abs } else { y_abs };
    let mut max_abs = if !recip { x_abs } else { y_abs };
    let mut min_exp = min_abs.wrapping_shr(52);
    let mut max_exp = max_abs.wrapping_shr(52);

    let mut num = f64::from_bits(min_abs);
    let mut den = f64::from_bits(max_abs);

    // Check for exceptional cases, whether inputs are 0, inf, nan, or close to
    // overflow, or close to underflow.
    if max_exp > 0x7ffu64 - 128u64 || min_exp < 128u64 {
        if x.is_nan() || y.is_nan() {
            return f64::NAN;
        }
        let x_except = if x == 0.0 {
            0
        } else if x.is_infinite() {
            2
        } else {
            1
        };
        let y_except = if y == 0.0 {
            0
        } else if y.is_infinite() {
            2
        } else {
            1
        };

        if (x_except != 1) || (y_except != 1) {
            let r = EXCEPTS[y_except][x_except][x_sign];
            return f64::mul_add(IS_NEG[y_sign], r.hi, IS_NEG[y_sign] * r.lo);
        }
        let scale_up = min_exp < 128u64;
        let scale_down = max_exp > 0x7ffu64 - 128u64;
        // At least one input is denormal, multiply both numerator and denominator
        // by some large enough power of 2 to normalize denormal inputs.
        // if scale_up || scale_down {
        //     return atan2_hard(y, x).fast_as_f64();
        // }
        if scale_up {
            num *= f64::from_bits(0x43f0000000000000);
            if !scale_down {
                den *= f64::from_bits(0x43f0000000000000);
            }
        } else if scale_down {
            den *= f64::from_bits(0x3bf0000000000000);
            if !scale_up {
                num *= f64::from_bits(0x3bf0000000000000);
            }
        }

        min_abs = num.to_bits();
        max_abs = den.to_bits();
        min_exp = min_abs.wrapping_shr(52);
        max_exp = max_abs.wrapping_shr(52);
    }
    let final_sign = IS_NEG[((x_sign != y_sign) != recip) as usize];
    let const_term = CONST_ADJ[x_sign][y_sign][recip as usize];
    let exp_diff = max_exp - min_exp;
    // We have the following bound for normalized n and d:
    //   2^(-exp_diff - 1) < n/d < 2^(-exp_diff + 1).
    if exp_diff > 54 {
        return f64::mul_add(
            final_sign,
            const_term.hi,
            final_sign * (const_term.lo + num / den),
        );
    }

    let mut k = (64.0 * num / den).round();
    let idx = k as u64;
    // k = idx / 64
    k *= f64::from_bits(0x3f90000000000000);

    // Range reduction:
    // atan(n/d) - atan(k/64) = atan((n/d - k/64) / (1 + (n/d) * (k/64)))
    //                        = atan((n - d * k/64)) / (d + n * k/64))
    let num_k = DoubleDouble::from_exact_mult_fma(num, k);
    let den_k = DoubleDouble::from_exact_mult_fma(den, k);

    // num_dd = n - d * k
    let num_dd = DoubleDouble::from_exact_add(num - den_k.hi, -den_k.lo);
    // den_dd = d + n * k
    let mut den_dd = DoubleDouble::from_exact_add(den, num_k.hi);
    den_dd.lo += num_k.lo;

    // q = (n - d * k) / (d + n * k)
    let q = DoubleDouble::div_fma(num_dd, den_dd);
    // p ~ atan(q)
    let p = atan_eval_fma(q);

    let vl = ATAN_I[idx as usize];
    let vlo = DoubleDouble::from_bit_pair(vl);
    let mut r = DoubleDouble::add(const_term, DoubleDouble::add(vlo, p));

    let err = f64::mul_add(
        p.hi,
        f64::from_bits(0x3bd0000000000000),
        f64::from_bits(0x3c00000000000000),
    );

    let ub = r.hi + (r.lo + err);
    let lb = r.hi + (r.lo - err);

    if ub == lb {
        r.hi *= final_sign;
        r.lo *= final_sign;

        return r.to_f64();
    }
    atan2_hard(y, x).fast_as_f64()
}

/// Computes atan(x)
///
/// Max found ULP 0.5
pub fn f_atan2(y: f64, x: f64) -> f64 {
    #[cfg(not(any(target_arch = "x86", target_arch = "x86_64")))]
    {
        atan2_gen_fma(y, x)
    }
    #[cfg(any(target_arch = "x86", target_arch = "x86_64"))]
    {
        use std::sync::OnceLock;
        static EXECUTOR: OnceLock<unsafe fn(f64, f64) -> f64> = OnceLock::new();
        let q = EXECUTOR.get_or_init(|| {
            if std::arch::is_x86_feature_detected!("avx")
                && std::arch::is_x86_feature_detected!("fma")
            {
                atan2_fma_impl
            } else {
                fn def_atan2pi(y: f64, x: f64) -> f64 {
                    atan2_gen_fma(y, x)
                }
                def_atan2pi
            }
        });
        unsafe { q(y, x) }
    }
}

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

    #[test]
    fn test_atan2() {
        assert_eq!(
            f_atan2(0.05474853958030223, 0.9999995380640253),
            0.05469396182367716
        );
        assert_eq!(f_atan2(-5., 2.), -1.1902899496825317);
        assert_eq!(f_atan2(2., -5.), 2.761086276477428);
        assert_eq!(
            f_atan2(1.220342145227879E-321, 6.9806238698201653E-309),
            0.00000000000017481849301519772
        );
    }
}

Coverage Report

Created: 2025-12-11 07:11