zenith-float-num 1.0.0

//! Components used in MPFR integration tests

use gmp_mpfr_sys::mpfr::{self, rnd_t};
use rand::rngs::StdRng;
use rand::{Rng, SeedableRng};
use rug::Float;
use std::cell::RefCell;
use zenith_float_num::{
    Consts, ExactNum, Exponent, Radix, RoundingMode, Sign, Word, WORD_BIT_SIZE,
    WORD_SIGNIFICANT_BIT,
};

thread_local! {
    static TEST_RNG: RefCell<StdRng> = RefCell::new(StdRng::seed_from_u64(0x5EED_CAFE_BADC_0D00));
}

pub fn test_random<T>() -> T
where
    rand::distributions::Standard: rand::distributions::Distribution<T>,
{
    TEST_RNG.with(|rng| rng.borrow_mut().gen())
}

pub fn reset_test_rng() {
    TEST_RNG.with(|rng| {
        *rng.borrow_mut() = StdRng::seed_from_u64(0x5EED_CAFE_BADC_0D00);
    });
    zenith_float_num::reseed_random(zenith_float_num::DEFAULT_RANDOM_SEED);
}

macro_rules! test_zf_op {
    ($eq:literal, $n1:ident, $n2:ident, $zf_op:ident, $f1:ident, $f2:ident, $mpfr_op:ident, $p:ident, $rm:ident, $rnd:ident, $op_info:expr, $cc:ident) => {
        let n3 = ExactNum::$zf_op(&($n1), &($n2), $p, $rm, &mut $cc);

        let mut f3 = Float::with_val($p as u32, 1);

        unsafe { mpfr::$mpfr_op(f3.as_raw_mut(), ($f1).as_raw(), ($f2).as_raw(), $rnd) };

        //println!("\n{:b}\n{:b}", $n1, $n2);
        //println!("\n{:b}\n{}", $n1, $f1.to_string_radix(2, None));
        //println!("\n{:b}\n{}", n3, f3.to_string_radix(2, None));

        assert_float_close(n3, f3, $p, &format!("{:?}", $op_info), $eq, &mut $cc);
    };
    ($eq:literal, $n1:ident, $zf_op:ident, $f1:ident, $mpfr_op:ident, $p:ident, $rm:ident, $rnd:ident, $op_info:expr, $cc:ident) => {
        let n3 = ExactNum::$zf_op(&($n1), $p, $rm, &mut $cc);

        let mut f3 = Float::with_val($p as u32, 1);

        unsafe { mpfr::$mpfr_op(f3.as_raw_mut(), ($f1).as_raw(), $rnd) };

        // println!("\n{:b}\n{}", $n1, $f1.to_string_radix(2, None));
        // println!("\n{:b}\n{}", n3, f3.to_string_radix(2, None));

        assert_float_close(n3, f3, $p, &format!("{:?}", $op_info), $eq, &mut $cc);
    };
    ($eq:literal, $n1:ident, $zf_op:ident, $f1:ident, $mpfr_op:ident, $p:ident, $rm:ident, $rnd:ident, $op_info:expr, $cc:ident) => {
        let n3 = ExactNum::$zf_op(&($n1), $p, $rm, &mut $cc);

        let mut f3 = Float::with_val($p as u32, 1);

        unsafe { mpfr::$mpfr_op(f3.as_raw_mut(), ($f1).as_raw(), $rnd) };

        // println!("\n{:b}\n{}", $n1, $f1.to_string_radix(2, None));
        // println!("\n{:b}\n{}", n3, f3.to_string_radix(2, None));

        assert_float_close(n3, f3, $p, &format!("{:?}", $op_info), $eq, &mut $cc);
    };
}

macro_rules! test_zf_op_no_cc {
    ($eq:literal, $n1:ident, $n2:ident, $zf_op:ident, $f1:ident, $f2:ident, $mpfr_op:ident, $p:ident, $rm:ident, $rnd:ident, $op_info:expr, $cc:ident) => {
        let n3 = ExactNum::$zf_op(&($n1), &($n2), $p, $rm);

        let mut f3 = Float::with_val($p as u32, 1);

        unsafe { mpfr::$mpfr_op(f3.as_raw_mut(), ($f1).as_raw(), ($f2).as_raw(), $rnd) };

        //println!("\n{:b}\n{:b}", $n1, $n2);
        //println!("\n{:b}\n{}", $n1, $f1.to_string_radix(2, None));
        //println!("\n{:b}\n{}", n3, f3.to_string_radix(2, None));

        assert_float_close(n3, f3, $p, &format!("{:?}", $op_info), $eq, &mut $cc);
    };
    ($eq:literal, $n1:ident, $zf_op:ident, $f1:ident, $mpfr_op:ident, $p:ident, $rm:ident, $rnd:ident, $op_info:expr, $cc:ident) => {
        let n3 = ExactNum::$zf_op(&($n1), $p, $rm);

        let mut f3 = Float::with_val($p as u32, 1);

        unsafe { mpfr::$mpfr_op(f3.as_raw_mut(), ($f1).as_raw(), $rnd) };

        // println!("\n{:b}\n{}", $n1, $f1.to_string_radix(2, None));
        // println!("\n{:b}\n{}", n3, f3.to_string_radix(2, None));

        assert_float_close(n3, f3, $p, &format!("{:?}", $op_info), $eq, &mut $cc);
    };
}

macro_rules! test_zf_fma {
    ($eq:literal, $n1:ident, $n2:ident, $nc:ident, $f1:ident, $f2:ident, $fc:ident, $p:ident, $rm:ident, $rnd:ident, $op_info:expr, $cc:ident) => {
        let n4 = ExactNum::fma(&($n1), &($n2), &($nc), $p, $rm);
        let mut f4 = Float::with_val($p as u32, 1);
        unsafe {
            mpfr::fma(
                f4.as_raw_mut(),
                ($f1).as_raw(),
                ($f2).as_raw(),
                ($fc).as_raw(),
                $rnd,
            )
        };
        assert_float_close(n4, f4, $p, &format!("{:?}", $op_info), $eq, &mut $cc);
    };
}

macro_rules! test_zf_rem_pi {
    ($n1:ident, $f1:ident, $p:ident, $rm:ident, $rnd:ident, $op_info:expr, $cc:ident) => {
        let mut n3 = ExactNum::rem_pi(&($n1), $p, $rm, &mut $cc);
        let _ = n3.set_precision($p, $rm);
        let _ = $rnd;

        if ($n1).exponent().unwrap_or(0) <= 2 {
            // Identity path: rem_pi leaves |x| < 4 unchanged aside from working precision.
            assert_float_close(
                n3,
                ($f1).clone(),
                $p,
                &format!("{:?}", $op_info),
                true,
                &mut $cc,
            );
        } else if !n3.is_nan() && !n3.is_inf() && !n3.is_zero() {
            // Not IEEE fmod; reduction is validated by sin/cos oracles on the unreduced argument.
            let e = n3.exponent().unwrap_or(0);
            assert!(
                e <= 3,
                "{}",
                format!("{:?} rem_pi exponent {e} (want <= 3)", $op_info)
            );
        }
    };
}

// test constant value match
macro_rules! test_zf_const {
    ($zf_const:ident, $mpfr_const:ident, $p:ident, $rm:ident, $rnd:ident, $op_info:expr, $cc:ident) => {
        let n1: ExactNum = $cc.$zf_const($p, $rm);

        let mut f1 = Float::with_val($p as u32, 1);

        unsafe {
            mpfr::$mpfr_const(f1.as_raw_mut(), $rnd);
        }

        assert_float_close(n1, f1, $p, &format!("{:?}", $op_info), true, &mut $cc);
    };
}

pub(crate) use test_zf_const;
pub(crate) use test_zf_fma;
pub(crate) use test_zf_op;
pub(crate) use test_zf_op_no_cc;
pub(crate) use test_zf_rem_pi;

pub const fn get_prec_rng() -> usize {
    157
}

pub fn get_float_pair(
    p: usize,
    emin: Exponent,
    emax: Exponent,
    cc: &mut Consts,
) -> (ExactNum, Float) {
    let n = ExactNum::random_normal(p, emin, emax);
    let f = conv_to_mpfr(p, &n, cc);
    (n, f)
}

pub fn conv_to_mpfr(p: usize, n: &ExactNum, cc: &mut Consts) -> Float {
    let s1 = conv_str_to_mpfr_compat(n.format(Radix::Bin, RoundingMode::ToEven, cc).unwrap());
    let f = Float::with_val(p as u32, Float::parse_radix(s1, 2).unwrap());
    let s2 = conv_str_from_mpfr_compat(f.to_string_radix(2, None));
    //println!("\n{}\n{}", s1, s2);
    assert_eq!(
        *n,
        ExactNum::parse(&s2, Radix::Bin, p, RoundingMode::None, cc)
    );
    f
}

// assert float values are equal or differ not more than 1 ulp.
pub fn assert_float_close(n: ExactNum, f: Float, p: usize, op: &str, eq: bool, cc: &mut Consts) {
    if n.is_inf() {
        // inf
        let ovf = unsafe { mpfr::overflow_p() };
        assert!(f.is_infinite() || ovf != 0, "{}", op);
        if eq {
            if n.is_positive() {
                assert!(f.is_sign_positive(), "{}", op);
            } else {
                assert!(f.is_sign_negative(), "{}", op);
            }
        }
    } else if n.is_nan() {
        // nan
        if eq {
            assert!(f.is_nan(), "{}", op);
        }
        return;
    } else if (n.is_subnormal() || n.is_zero()) && eq {
        // subnormal (bit-exact mode only — otherwise compare below)
        let unf = unsafe { mpfr::underflow_p() };
        if !f.is_zero() {
            assert!(unf != 0, "{}", op);
        }
    } else if eq {
        // n == f
        let s1 = f.to_string_radix(2, None);
        let s2 = conv_str_from_mpfr_compat(s1);
        let n2 = ExactNum::parse(&s2, Radix::Bin, p, RoundingMode::None, cc);
        assert_eq!(n, n2, "{}", op);
    } else {
        // at most 1 ulp difference
        let s1 = f.to_string_radix(2, None);
        let s2 = conv_str_from_mpfr_compat(s1);
        let n2 = ExactNum::parse(&s2, Radix::Bin, p, RoundingMode::None, cc);
        assert_eq!(
            n.mantissa_max_bit_len(),
            n2.mantissa_max_bit_len(),
            "{}",
            op
        );

        let d = n.sub(&n2, 1, RoundingMode::None).abs();

        let e1 = n.exponent().unwrap() as isize
            - (n.mantissa_max_bit_len().unwrap() - n.precision().unwrap()) as isize;
        let e2 = d.exponent().unwrap() as isize
            - (d.mantissa_max_bit_len().unwrap() - d.precision().unwrap()) as isize;

        //println!("\n{:?}\n{:?}\n{:?}\n{} {}", n, n2, d, e1, e2);

        if !d.is_zero() {
            assert!(
                e1 - e2 >= n.mantissa_max_bit_len().unwrap() as isize - 1,
                "{}",
                op
            );
        }
    }
}

pub fn conv_str_to_mpfr_compat(s: String) -> String {
    let (sig, exp) = if let Some(pos) = s.find('e') {
        s.split_at(pos + 1)
    } else {
        (s.as_str(), "0")
    };
    let expn = i64::from_str_radix(exp, 2).unwrap();
    sig.to_owned() + &expn.to_string()
}

pub fn conv_str_from_mpfr_compat(s: String) -> String {
    if let Some(epos) = s.find('e') {
        let (sig, exp) = s.split_at(epos + 1);
        let expn = exp.parse::<i64>().unwrap();
        if expn < 0 {
            sig.to_owned() + "-" + &format!("{:b}", -expn)
        } else {
            sig.to_owned() + &format!("{:b}", expn)
        }
    } else {
        s
    }
}

pub fn get_random_rnd_pair() -> (RoundingMode, rnd_t) {
    match test_random::<u8>() % 5 {
        0 => (RoundingMode::ToEven, rnd_t::RNDN),
        1 => (RoundingMode::Up, rnd_t::RNDU),
        2 => (RoundingMode::Down, rnd_t::RNDD),
        3 => (RoundingMode::FromZero, rnd_t::RNDA),
        4 => (RoundingMode::ToZero, rnd_t::RNDZ),
        _ => unreachable!(),
    }
}

// Generates a number with mantissa like 1111111..1110000..00000
pub fn get_oned_zeroed(p: usize, exp_from: Exponent, exp_to: Exponent) -> ExactNum {
    let zero_bits = test_random::<usize>() % (p - 1);
    let mut m1 = vec![Word::MAX; p / WORD_BIT_SIZE];
    let i = zero_bits / WORD_BIT_SIZE;
    m1.iter_mut().take(i).for_each(|v| *v = 0);
    m1[i] <<= zero_bits % WORD_BIT_SIZE;

    bf_from_mantissa_and_exp_rng(&m1, exp_from, exp_to)
}

// Generates a number with mantissa like 111..111000..000111..111
pub fn get_oned_sides(p: usize, exp_from: Exponent, exp_to: Exponent) -> ExactNum {
    let one_bits = test_random::<usize>() % (p / 2).min(WORD_BIT_SIZE * 2) + 1;
    let mut m1 = vec![0; p / WORD_BIT_SIZE];
    let i = one_bits / WORD_BIT_SIZE;
    let m1l = m1.len();
    m1.iter_mut().take(i + 1).for_each(|v| *v = Word::MAX);
    m1[i] >>= one_bits % WORD_BIT_SIZE;
    m1.iter_mut().rev().take(i + 1).for_each(|v| *v = Word::MAX);
    m1[m1l - i - 1] <<= one_bits % WORD_BIT_SIZE;

    bf_from_mantissa_and_exp_rng(&m1, exp_from, exp_to)
}

// Generates a number with periodic mantissa
pub fn get_periodic(p: usize, exp_from: Exponent, exp_to: Exponent) -> ExactNum {
    let nbits = 4;
    let bits = test_random::<Word>() % ((1 << nbits) - 1) + 1;
    let mut w: Word = bits;
    for _ in 1..WORD_BIT_SIZE / nbits {
        w <<= nbits;
        w |= bits;
    }
    let mut m1 = vec![w; p / WORD_BIT_SIZE];
    *(m1.last_mut().unwrap()) |= 1 << (WORD_BIT_SIZE - 1);

    bf_from_mantissa_and_exp_rng(&m1, exp_from, exp_to)
}

// Generates a number with oned mantissa, but last bit is zero
pub fn get_last_zero(p: usize, exp_from: Exponent, exp_to: Exponent) -> ExactNum {
    let mut m1 = vec![Word::MAX; p / WORD_BIT_SIZE];
    m1[0] = Word::MAX - 1;

    bf_from_mantissa_and_exp_rng(&m1, exp_from, exp_to)
}

// Generates a number near 1.
pub fn get_near_one(p: usize) -> ExactNum {
    let e = (test_random::<u8>() & 1) as Exponent;

    let random_bits = test_random::<usize>() % (p - 1);
    let i = random_bits / WORD_BIT_SIZE;

    let mut m1;
    if e == 0 {
        m1 = vec![Word::MAX; p / WORD_BIT_SIZE];
    } else {
        m1 = vec![0; p / WORD_BIT_SIZE];
    }
    m1[i] ^= test_random::<Word>() >> (random_bits % WORD_BIT_SIZE);

    m1.iter_mut().take(i).for_each(|v| *v = test_random());

    *m1.last_mut().unwrap() |= WORD_SIGNIFICANT_BIT;

    bf_from_mantissa_and_exp_rng(&m1, e, e)
}

pub fn bf_from_mantissa_and_exp_rng(m: &[Word], exp_from: Exponent, exp_to: Exponent) -> ExactNum {
    let e = if exp_from < exp_to {
        (test_random::<isize>().abs() % (exp_to as isize - exp_from as isize) + exp_from as isize)
            as Exponent
    } else {
        exp_from
    };

    let s = if test_random::<u8>() & 1 == 0 { Sign::Pos } else { Sign::Neg };

    ExactNum::from_words(m, s, e)
}