pxfm 0.1.29

/*
 * // 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::bits::set_exponent_f64;
use crate::common::{dd_fmla, f_fmla};
use crate::double_double::DoubleDouble;
use crate::dyadic_float::{DyadicFloat128, DyadicSign};
use crate::rounding::CpuRound;
use crate::sincos_reduce_tables::ONE_TWENTY_EIGHT_OVER_PI;

#[derive(Debug)]
pub(crate) struct AngleReduced {
    pub(crate) angle: DoubleDouble,
}

#[derive(Default)]
pub(crate) struct LargeArgumentReduction {
    x_reduced: f64,
    idx: u64,
    y_hi: f64,
    y_lo: f64,
    // Low part of x * ONE_TWENTY_EIGHT_OVER_PI[idx][1].
    y_mid: DoubleDouble,
}

// For large range |x| >= 2^16, we perform the range reduction computations as:
//   u = x - k * pi/128 = (pi/128) * (x * (128/pi) - k).
// We use the exponent of x to find 4 double-chunks of 128/pi:
// c_hi, c_mid, c_lo, c_lo_2 such that:
//   1) ulp(round(x * c_hi, D, RN)) >= 2^8 = 256,
//   2) If x * c_hi = ph_hi + ph_lo and x * c_mid = pm_hi + pm_lo, then
//        min(ulp(ph_lo), ulp(pm_hi)) >= 2^-53.
// This will allow us to drop the high part ph_hi and the addition:
//   (ph_lo + pm_hi) mod 1
// can be exactly representable in a double precision.
// This will allow us to do split the computations as:
//   (x * 256/pi) ~ x * (c_hi + c_mid + c_lo + c_lo_2)    (mod 256)
//                ~ (ph_lo + pm_hi) + (pm_lo + x * c_lo) + x * c_lo_2.
// Then,
//   round(x * 128/pi) = round(ph_lo + pm_hi)    (mod 256)
// And the high part of fractional part of (x * 128/pi) can simply be:
//   {x * 128/pi}_hi = {ph_lo + pm_hi}.
// To prevent overflow when x is very large, we simply scale up
// (c_hi, c_mid, c_lo, c_lo_2) by a fixed power of 2 (based on the index) and
// scale down x by the same amount.
impl LargeArgumentReduction {
    #[cold]
    pub(crate) fn accurate(&self) -> DyadicFloat128 {
        // Sage math:
        // R = RealField(128)
        // π = R.pi()
        //
        // def format_hex(value):
        //     l = hex(value)[2:]
        //     n = 4
        //     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_hex(m)},")
        //     print("};")
        //
        // print_dyadic(π/128)
        const PI_OVER_128_F128: DyadicFloat128 = DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -133,
            mantissa: 0xc90f_daa2_2168_c234_c4c6_628b_80dc_1cd1_u128,
        };

        // y_lo = x * c_lo + pm.lo
        let one_pi_rot = ONE_TWENTY_EIGHT_OVER_PI[self.idx as usize];
        let y_lo_0 = DyadicFloat128::new_from_f64(self.x_reduced * f64::from_bits(one_pi_rot.3));
        let y_lo_1 = DyadicFloat128::new_from_f64(self.y_lo) + y_lo_0;
        let y_mid_f128 = DyadicFloat128::new_from_f64(self.y_mid.lo) + y_lo_1;
        let y_hi_f128 =
            DyadicFloat128::new_from_f64(self.y_hi) + DyadicFloat128::new_from_f64(self.y_mid.hi);
        let y = y_hi_f128 + y_mid_f128;

        y * PI_OVER_128_F128
    }

    pub(crate) fn reduce(&mut self, x: f64) -> (u64, DoubleDouble) {
        const E_BIAS: u64 = (1u64 << (11 - 1u64)) - 1u64;
        let mut xbits = x.to_bits();
        let x_e = ((x.to_bits() >> 52) & 0x7ff) as i64;
        let x_e_m62 = x_e.wrapping_sub(E_BIAS as i64 + 62);
        self.idx = (x_e_m62 >> 4).wrapping_add(3) as u64;
        // Scale x down by 2^(-(16 * (idx - 3))
        xbits = set_exponent_f64(
            xbits,
            (x_e_m62 & 15)
                .wrapping_add(E_BIAS as i64)
                .wrapping_add(62i64) as u64,
        );
        // 2^62 <= |x_reduced| < 2^(62 + 16) = 2^78
        self.x_reduced = f64::from_bits(xbits);
        // x * c_hi = ph.hi + ph.lo exactly.
        let one_pi = ONE_TWENTY_EIGHT_OVER_PI[self.idx as usize];
        let ph = DoubleDouble::from_exact_mult(self.x_reduced, f64::from_bits(one_pi.0));
        // x * c_mid = pm.hi + pm.lo exactly.
        let pm = DoubleDouble::from_exact_mult(self.x_reduced, f64::from_bits(one_pi.1));
        // x * c_lo = pl.hi + pl.lo exactly.
        let pl = DoubleDouble::from_exact_mult(self.x_reduced, f64::from_bits(one_pi.2));
        // Extract integral parts and fractional parts of (ph.lo + pm.hi).
        let sum_hi = ph.lo + pm.hi;
        let kd = sum_hi.cpu_round();

        // x * 128/pi mod 1 ~ y_hi + y_mid + y_lo
        self.y_hi = (ph.lo - kd) + pm.hi; // Exact
        self.y_mid = DoubleDouble::from_exact_add(pm.lo, pl.hi);
        self.y_lo = pl.lo;

        // y_l = x * c_lo_2 + pl.lo
        let y_l = dd_fmla(self.x_reduced, f64::from_bits(one_pi.3), self.y_lo);
        let mut y = DoubleDouble::from_exact_add(self.y_hi, self.y_mid.hi);
        y.lo += self.y_mid.lo + y_l;

        // Digits of pi/128, generated by SageMath with:
        // import struct
        // from sage.all import *
        //
        // def double_to_hex(f):
        //     return "0x" + format(struct.unpack('<Q', struct.pack('<d', f))[0], '016x')
        //
        // R = RealField(128)
        // π = R.pi()
        //
        // RN = RealField(53)
        //
        // hi = RN(π/128)
        // lo = RN(π/128 - R(hi))
        //
        // print("lo: " + double_to_hex(lo))
        // print("hi: " + double_to_hex(hi))
        const PI_OVER_128_DD: DoubleDouble = DoubleDouble::new(
            f64::from_bits(0x3c31a62633145c07),
            f64::from_bits(0x3f9921fb54442d18),
        );

        // Error bound: with {a} denote the fractional part of a, i.e.:
        //   {a} = a - round(a)
        // Then,
        //   | {x * 128/pi} - (y_hi + y_lo) | <=  ulp(ulp(y_hi)) <= 2^-105
        //   | {x mod pi/128} - (u.hi + u.lo) | < 2 * 2^-6 * 2^-105 = 2^-110
        let u = DoubleDouble::quick_mult(y, PI_OVER_128_DD);

        (kd as i64 as u64, u)
    }
}

// Generated by SageMath
// nwords = 20
// prec = nwords * 64 + 150
// R = RealField(prec)
// invpi = R(1) / (R(2) * R.pi())
//
// scale = R(2)**64
//
// words = []
// x = invpi
// for i in range(nwords):
//     y = floor(x * scale)
//     words.append(int(y))
//     x = x * scale - y
//
// for w in words:
//     print("0x{:016x},".format(w))
static INVPI_2_64: [u64; 20] = [
    0x28be60db9391054a,
    0x7f09d5f47d4d3770,
    0x36d8a5664f10e410,
    0x7f9458eaf7aef158,
    0x6dc91b8e909374b8,
    0x1924bba82746487,
    0x3f877ac72c4a69cf,
    0xba208d7d4baed121,
    0x3a671c09ad17df90,
    0x4e64758e60d4ce7d,
    0x272117e2ef7e4a0e,
    0xc7fe25fff7816603,
    0xfbcbc462d6829b47,
    0xdb4d9fb3c9f2c26d,
    0xd3d18fd9a797fa8b,
    0x5d49eeb1faf97c5e,
    0xcf41ce7de294a4ba,
    0x9afed7ec47e35742,
    0x1580cc11bf1edaea,
    0xfc33ef0826bd0d87,
];

#[inline]
fn create_dd(c1: u64, c0: u64) -> DoubleDouble {
    let mut c1 = c1;
    let mut c0 = c0;
    if c1 != 0 {
        let e = c1.leading_zeros();
        if e != 0 {
            c1 = (c1 << e) | (c0 >> (64 - e));
            c0 = c0.wrapping_shl(e);
        }
        let f = 0x3fe - e;
        let t_u = ((f as u64) << 52) | ((c1 << 1) >> 12);
        let hi = f64::from_bits(t_u);
        c0 = (c1 << 53) | (c0 >> 11);
        let l = if c0 != 0 {
            let g = c0.leading_zeros();
            if (g) != 0 {
                c0 = c0.wrapping_shl(g);
            }
            let t_u = (((f - 53 - g) as u64) << 52) | ((c0 << 1) >> 12);
            f64::from_bits(t_u)
        } else {
            0.
        };
        DoubleDouble::new(l, hi)
    } else if c0 != 0 {
        let e = c0.leading_zeros();
        let f = 0x3fe - 64 - e;
        c0 = c0.wrapping_shl(e + 1); // most significant bit shifted out

        /* put the upper 52 bits of c0 into h */
        let t_u = ((f as u64) << 52) | (c0 >> 12);
        let hi = f64::from_bits(t_u);
        /* put the lower 12 bits of c0 into l */
        c0 = c0.wrapping_shl(52);
        let l = if c0 != 0 {
            let g = c0.leading_zeros();
            c0 = c0.wrapping_shl(g + 1);
            let t_u = (((f - 64 - g) as u64) << 52) | (c0 >> 12);
            f64::from_bits(t_u)
        } else {
            0.
        };
        DoubleDouble::new(l, hi)
    } else {
        DoubleDouble::default()
    }
}

#[inline]
fn frac_2pi(x: f64) -> DoubleDouble {
    if x <= f64::from_bits(0x401921fb54442d17)
    // x < 2*pi
    {
        /* | CH+CL - 1/(2pi) | < 2^-110.523 */
        const C: DoubleDouble = DoubleDouble::new(
            f64::from_bits(0xbc66b01ec5417056),
            f64::from_bits(0x3fc45f306dc9c883),
        );
        let mut z = DoubleDouble::quick_mult_f64(C, x);
        z.lo = f_fmla(C.lo, x, z.lo);
        z
    } else
    // x > 0x1.921fb54442d17p+2
    {
        let t = x.to_bits();
        let mut e = ((t >> 52) & 0x7ff) as i32; /* 1025 <= e <= 2046 */
        let m = (1u64 << 52) | (t & 0xfffffffffffffu64);
        let mut c0: u64;
        let mut c1: u64;
        let mut c2: u64;
        // x = m/2^53 * 2^(e-1022)
        if e <= 1074
        // 1025 <= e <= 1074: 2^2 <= x < 2^52
        {
            let mut u = m as u128 * INVPI_2_64[1] as u128;
            c0 = u as u64;
            c1 = (u >> 64) as u64;
            u = m as u128 * INVPI_2_64[0] as u128;
            c1 = c1.wrapping_add(u as u64);
            c2 = (u >> 64) as u64 + (c1 < (u as u64)) as u64;
            e = 1075 - e; // 1 <= e <= 50
        } else
        // 1075 <= e <= 2046, 2^52 <= x < 2^1024
        {
            let i = (e - 1138 + 63) / 64; // i = ceil((e-1138)/64), 0 <= i <= 15
            let mut u = m as u128 * INVPI_2_64[i as usize + 2] as u128;
            c0 = u as u64;
            c1 = (u >> 64) as u64;
            u = m as u128 * INVPI_2_64[i as usize + 1] as u128;
            c1 = c1.wrapping_add(u as u64);
            c2 = (u >> 64) as u64 + ((c1) < (u as u64)) as u64;
            u = m as u128 * INVPI_2_64[i as usize] as u128;
            c2 = c2.wrapping_add(u as u64);
            e = 1139 + (i << 6) - e; // 1 <= e <= 64
        }
        if e == 64 {
            c0 = c1;
            c1 = c2;
        } else {
            c0 = (c1 << (64 - e)) | c0 >> e;
            c1 = (c2 << (64 - e)) | c1 >> e;
        }
        create_dd(c1, c0)
    }
}

/// Returns x mod 2*PI
#[inline]
pub(crate) fn rem2pi_any(x: f64) -> AngleReduced {
    const TWO_PI: DoubleDouble = DoubleDouble::new(
        f64::from_bits(0x3cb1a62633145c07),
        f64::from_bits(0x401921fb54442d18),
    );
    let a = frac_2pi(x);
    let z = DoubleDouble::quick_mult(a, TWO_PI);
    AngleReduced { angle: z }
}

/**
Generated by SageMath:
```python
def triple_double_split(x, n):
    VR = RealField(1000)
    R32 = RealField(n)
    hi = R32(x)
    rem1 = x - VR(hi)
    med = R32(rem1)
    rem2 = rem1 - VR(med)
    lo = R32(rem2)
    rem2 = rem1 - VR(med)
    lo = R32(rem2)
    return (hi, med, lo)

hi, med, lo = triple_double_split((RealField(600).pi() * RealField(600)(2)), 51)

print(double_to_hex(hi))
print(double_to_hex(med))
print(double_to_hex(lo))
```
 **/
#[inline]
fn rem2pif_small(x: f32) -> f64 {
    const ONE_OVER_PI_2: f64 = f64::from_bits(0x3fc45f306dc9c883);
    const TWO_PI: [u64; 3] = [0x401921fb54440000, 0x3da68c234c4c0000, 0x3b498a2e03707345];
    let dx = x as f64;
    let kd = (dx * ONE_OVER_PI_2).cpu_round();
    let mut y = dd_fmla(-kd, f64::from_bits(TWO_PI[0]), dx);
    y = dd_fmla(-kd, f64::from_bits(TWO_PI[1]), y);
    y = dd_fmla(-kd, f64::from_bits(TWO_PI[2]), y);
    y
}

#[inline]
pub(crate) fn rem2pif_any(x: f32) -> f64 {
    let x_abs = x.abs();
    if x_abs.to_bits() < 0x5600_0000u32 {
        rem2pif_small(x_abs)
    } else {
        let p = rem2pi_any(x_abs as f64);
        p.angle.to_f64()
    }
}

pub(crate) fn frac2pi_d128(x: DyadicFloat128) -> DyadicFloat128 {
    let e = x.biased_exponent();

    let mut fe = x;

    if e <= 1
    // |X| < 2
    {
        /* multiply by T[0]/2^64 + T[1]/2^128, where
        |T[0]/2^64 + T[1]/2^128 - 1/(2pi)| < 2^-130.22 */
        let mut x_hi = (x.mantissa >> 64) as u64;
        let mut x_lo: u64;
        let mut u: u128 = x_hi as u128 * INVPI_2_64[1] as u128;
        let tiny: u64 = u as u64;
        x_lo = (u >> 64) as u64;
        u = x_hi as u128 * INVPI_2_64[0] as u128;
        x_lo = x_lo.wrapping_add(u as u64);
        x_hi = (u >> 64) as u64 + (x_lo < u as u64) as u64;
        /* hi + lo/2^64 + tiny/2^128 = hi_in * (T[0]/2^64 + T[1]/2^128) thus
        |hi + lo/2^64 + tiny/2^128 - hi_in/(2*pi)| < hi_in * 2^-130.22
        Since X is normalized at input, hi_in >= 2^63, and since T[0] >= 2^61,
        we have hi >= 2^(63+61-64) = 2^60, thus the normalize() below
        perform a left shift by at most 3 bits */
        let mut e = x.exponent;
        fe.mantissa = (x_hi as u128).wrapping_shl(64) | (x_lo as u128);
        fe.normalize();
        e -= fe.exponent;
        // put the upper e bits of tiny into X->lo
        if (e) != 0 {
            x_hi = (fe.mantissa >> 64) as u64;
            x_lo = (fe.mantissa & 0xffff_ffff_ffff_ffff) as u64;
            x_lo |= tiny >> (64 - e);
            fe.mantissa = (x_hi as u128).wrapping_shl(64) | (x_lo as u128);
        }
        /* The error is bounded by 2^-130.22 (relative) + ulp(lo) (absolute).
           Since now X->hi >= 2^63, the absolute error of ulp(lo) converts into
           a relative error of less than 2^-127.
           This yields a maximal relative error of:
           (1 + 2^-130.22) * (1 + 2^-127) - 1 < 2^-126.852.
        */
        return fe;
    }

    // now 2 <= e <= 1024

    /* The upper 64-bit word X->hi corresponds to hi/2^64*2^e, if multiplied by
    T[i]/2^((i+1)*64) it yields hi*T[i]/2^128 * 2^(e-i*64).
    If e-64i <= -128, it contributes to less than 2^-128;
    if e-64i >= 128, it yields an integer, which is 0 modulo 1.
    We thus only consider the values of i such that -127 <= e-64i <= 127,
    i.e., (-127+e)/64 <= i <= (127+e)/64.
    Up to 4 consecutive values of T[i] can contribute (only 3 when e is a
    multiple of 64). */
    let i = (if e < 127 { 0 } else { (e - 127 + 64 - 1) / 64 }) as usize; // ceil((e-127)/64)
    // 0 <= i <= 15
    let mut c: [u64; 5] = [0u64; 5];

    let mut x_hi = (x.mantissa >> 64) as u64;
    let mut x_lo: u64;

    let mut u: u128 = x_hi as u128 * INVPI_2_64[i + 3] as u128; // i+3 <= 18
    c[0] = u as u64;
    c[1] = (u >> 64) as u64;
    u = x_hi as u128 * INVPI_2_64[i + 2] as u128;
    c[1] = c[1].wrapping_add(u as u64);
    c[2] = (u >> 64) as u64 + (c[1] < u as u64) as u64;
    u = x_hi as u128 * INVPI_2_64[i + 1] as u128;
    c[2] = c[2].wrapping_add(u as u64);
    c[3] = (u >> 64) as u64 + (c[2] < u as u64) as u64;
    u = x_hi as u128 * INVPI_2_64[i] as u128;
    c[3] = c[3].wrapping_add(u as u64);
    c[4] = (u >> 64) as u64 + (c[3] < u as u64) as u64;

    let f = e as i32 - 64 * i as i32; // hi*T[i]/2^128 is multiplied by 2^f

    /* {c, 5} = hi*(T[i]+T[i+1]/2^64+T[i+2]/2^128+T[i+3]/2^192) */
    /* now shift c[0..4] by f bits to the left */
    let tiny;
    if f < 64 {
        x_hi = (c[4] << f) | (c[3] >> (64 - f));
        x_lo = (c[3] << f) | (c[2] >> (64 - f));
        tiny = (c[2] << f) | (c[1] >> (64 - f));
        /* the ignored part was less than 1 in c[1],
        thus less than 2^(f-64) <= 1/2 in tiny */
    } else if f == 64 {
        x_hi = c[3];
        x_lo = c[2];
        tiny = c[1];
        /* the ignored part was less than 1 in c[1],
        thus less than 1 in tiny */
    } else
    /* 65 <= f <= 127: this case can only occur when e >= 65 */
    {
        let g = f - 64; /* 1 <= g <= 63 */
        /* we compute an extra term */
        u = x_hi as u128 * INVPI_2_64[i + 4] as u128; // i+4 <= 19
        u >>= 64;
        c[0] = c[0].wrapping_add(u as u64);
        c[1] += (c[0] < u as u64) as u64;
        c[2] += ((c[0] < u as u64) && c[1] == 0) as u64;
        c[3] += ((c[0] < u as u64) && c[1] == 0 && c[2] == 0) as u64;
        c[4] += ((c[0] < u as u64) && c[1] == 0 && c[2] == 0 && c[3] == 0) as u64;
        x_hi = (c[3] << g) | (c[2] >> (64 - g));
        x_lo = (c[2] << g) | (c[1] >> (64 - g));
        tiny = (c[1] << g) | (c[0] >> (64 - g));
        /* the ignored part was less than 1 in c[0],
        thus less than 1/2 in tiny */
    }
    let mut fe = x;
    fe.exponent = -127;
    fe.mantissa = (x_hi as u128).wrapping_shl(64) | (x_lo as u128);
    fe.normalize();
    let ze = fe.biased_exponent();
    if ze < 0 {
        x_hi = (fe.mantissa >> 64) as u64;
        x_lo = (fe.mantissa & 0xffff_ffff_ffff_ffff) as u64;
        x_lo |= tiny >> (64 + ze);
        fe.mantissa = (x_hi as u128).wrapping_shl(64) | (x_lo as u128);
    }
    fe
}

pub(crate) fn rem2pi_f128(x: DyadicFloat128) -> DyadicFloat128 {
    /*
        Generated by SageMath:
        ```python
    def double_to_hex(f):
        # Converts Python float (f64) to hex string
        packed = struct.pack('>d', float(f))
        return '0x' + packed.hex()

    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("},")

    print_dyadic(RealField(300).pi() * 2)
        ```
         */
    const TWO_PI: DyadicFloat128 = DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -125,
        mantissa: 0xc90fdaa2_2168c234_c4c6628b_80dc1cd1_u128,
    };
    let frac2pi = frac2pi_d128(x);
    TWO_PI * frac2pi
}
//
// #[cfg(test)]
// mod tests {
//     use crate::dyadic_float::DyadicFloat128;
//     use crate::sincos_reduce::{frac_2pi, frac2pi_d128};
//
//     #[test]
//     fn test_reduce() {
//         let x = 1.8481363e36f32;
//         let reduced = frac_2pi(x as f64);
//         let to_reduced2 = DyadicFloat128::new_from_f64(x as f64);
//         let reduced2 = frac2pi_d128(to_reduced2);
//         println!("{:?}", reduced);
//         println!("{:?}", reduced2.fast_as_f64());
//     }
// }