zenith_float_num/common/

util.rs

//! Auxiliary functions.

use crate::{
    defs::{Error, Word, WORD_BIT_SIZE, WORD_MAX, WORD_SIGNIFICANT_BIT},
    RoundingMode,
};

#[cfg(test)]
use crate::{num::ExactNumNumber, Sign, EXPONENT_MIN};

#[cfg(test)]
#[cfg(not(feature = "std"))]
use alloc::vec::Vec;

/// integer logarithm base 2 of a number.
pub fn log2_ceil(mut n: usize) -> usize {
    let mut ret = 0;
    let mut sticky = 0;
    while n > 1 {
        if n & 1 != 0 {
            sticky = 1;
        }
        ret += 1;
        n >>= 1;
    }
    ret + sticky
}

/// integer logarithm base 2 of a number.
pub fn log2_floor(mut n: usize) -> usize {
    let mut ret = 0;
    while n > 1 {
        ret += 1;
        n >>= 1;
    }
    ret
}

/// square root integer approximation.
pub fn sqrt_int(a: u32) -> u32 {
    let a = a as u64;
    let mut x = a;
    for _ in 0..20 {
        if x == 0 {
            break;
        }
        x = (a / x + x) >> 1;
    }
    x as u32
}

// cost of multiplication of two numbers with precision p.
pub fn calc_mul_cost(p: usize) -> usize {
    if p < 70 {
        p * p
    } else {
        // toom-3
        if p < 1625 {
            sqrt_int((p * p * p) as u32) as usize
        } else {
            let q = sqrt_int(p as u32) as usize;
            q * q * q
        }
    }
}

// cost of addition/subtraction of two numbers with precision p.
#[inline]
pub fn calc_add_cost(p: usize) -> usize {
    p
}

/// Maximum extra correct-rounding retries.
pub const MAX_PREC_RETRY: usize = 256;

#[cfg(test)]
pub const MAX_DEC_SCALE: usize = 1_000_000;

/// True when working precision has grown too far past the requested precision.
#[inline]
pub fn prec_retry_exhausted(p_wrk: usize, p: usize) -> bool {
    p_wrk > p.saturating_add(WORD_BIT_SIZE.saturating_mul(MAX_PREC_RETRY))
}

/// Bump working precision for a correct-rounding retry, or error instead of spinning.
#[inline]
pub fn bump_prec_retry(p_wrk: &mut usize, p_inc: &mut usize, p: usize) -> Result<(), Error> {
    if prec_retry_exhausted(*p_wrk, p) {
        return Err(Error::PrecisionRetryExhausted);
    }
    *p_wrk += *p_inc;
    *p_inc = round_p(*p_wrk / 5);
    Ok(())
}

#[cfg(test)]
pub const TEST_EXP_BOUND: crate::Exponent = 1024;

#[cfg(test)]
pub fn test_loop_count(full: usize) -> usize {
    full.min(200)
}

// Estimate of sqrt op cost.
#[inline]
pub fn calc_sqrt_cost(p: usize, cost_mul: usize, cost_add: usize) -> usize {
    let log3_estimate = (log2_floor(p) * 41349) >> 16;
    log3_estimate * (5 * cost_mul + 2 * cost_add) / 2
}

#[inline(always)]
pub fn add_carry(a: Word, b: Word, c: Word, r: &mut Word) -> Word {
    #[cfg(target_arch = "x86_64")]
    {
        // platform-specific operation
        core::arch::x86_64::_addcarry_u64(c as u8, a, b, r) as Word
    }

    #[cfg(target_arch = "x86")]
    {
        // platform-specific operation
        core::arch::x86::_addcarry_u32(c as u8, a, b, r) as Word
    }

    #[cfg(not(any(target_arch = "x86_64", target_arch = "x86")))]
    {
        use crate::defs::DoubleWord;
        use crate::defs::WORD_BASE;

        let mut s = c as DoubleWord + a as DoubleWord + b as DoubleWord;
        if s >= WORD_BASE {
            s -= WORD_BASE;
            *r = s as Word;
            1
        } else {
            *r = s as Word;
            0
        }
    }
}

#[inline(always)]
pub fn sub_borrow(a: Word, b: Word, c: Word, r: &mut Word) -> Word {
    #[cfg(target_arch = "x86_64")]
    {
        // platform-specific operation
        core::arch::x86_64::_subborrow_u64(c as u8, a, b, r) as Word
    }

    #[cfg(target_arch = "x86")]
    {
        // platform-specific operation
        core::arch::x86::_subborrow_u32(c as u8, a, b, r) as Word
    }

    #[cfg(not(any(target_arch = "x86_64", target_arch = "x86")))]
    {
        use crate::defs::DoubleWord;
        use crate::defs::WORD_BASE;

        let v1 = a as DoubleWord;
        let v2 = b as DoubleWord + c as DoubleWord;

        if v1 < v2 {
            *r = (v1 + WORD_BASE - v2) as Word;
            1
        } else {
            *r = (v1 - v2) as Word;
            0
        }
    }
}

// Shift m left by n digits.
pub fn shift_slice_left(m: &mut [Word], n: usize) {
    let idx = n / WORD_BIT_SIZE;
    let shift = n % WORD_BIT_SIZE;
    if idx >= m.len() {
        m.fill(0);
    } else if shift > 0 {
        let l = m.len() - 1;
        let end = m.as_mut_ptr();
        unsafe {
            let mut dst = end.add(l);
            let mut src = end.add(l - idx);
            loop {
                if src > end {
                    let mut d = *src << shift;
                    src = src.sub(1);
                    d |= *src >> (WORD_BIT_SIZE - shift);
                    *dst = d;
                    dst = dst.sub(1);
                } else {
                    break;
                }
            }
            *dst = *src << shift;
        }
        m[0..idx].fill(0);
    } else if idx > 0 {
        let r = m.len() - idx;
        m.copy_within(0..r, idx);
        m[..idx].fill(0);
    }
}

// Shift m left by n digits and put result in m2.
pub fn shift_slice_left_copy(m: &[Word], m2: &mut [Word], n: usize) {
    let idx = n / WORD_BIT_SIZE;
    let shift = n % WORD_BIT_SIZE;
    if idx >= m2.len() {
        m2.fill(0);
    } else if shift > 0 {
        m2[..idx].fill(0);
        let mut dst = m2[idx..].iter_mut();
        let src = m.iter();
        let mut prev = 0;
        for (a, b) in src.zip(dst.by_ref()) {
            *b = (prev >> (WORD_BIT_SIZE - shift)) | (*a << shift);
            prev = *a;
        }
        if let Some(b) = dst.next() {
            *b = prev >> (WORD_BIT_SIZE - shift);
        }
        for b in dst {
            *b = 0;
        }
    } else {
        m2[..idx].fill(0);
        let mut dst = m2[idx..].iter_mut();
        for (a, b) in m.iter().zip(dst.by_ref()) {
            *b = *a;
        }
        for b in dst {
            *b = 0;
        }
    }
}

// Shift m right by n digits.
pub fn shift_slice_right(m: &mut [Word], n: usize) {
    let idx = n / WORD_BIT_SIZE;
    let shift = n % WORD_BIT_SIZE;
    if idx >= m.len() {
        m.fill(0);
    } else if shift > 0 {
        let l = m.len();
        let mut dst = m.as_mut_ptr();
        unsafe {
            let end = dst.add(l - 1);
            let mut src = dst.add(idx);
            loop {
                if src < end {
                    let mut d = *src >> shift;
                    src = src.add(1);
                    d |= *src << (WORD_BIT_SIZE - shift);
                    *dst = d;
                    dst = dst.add(1);
                } else {
                    break;
                }
            }
            *dst = *src >> shift;
        }
        m[l - idx..].fill(0);
    } else if idx > 0 {
        let r = m.len() - idx;
        m.copy_within(idx.., 0);
        m[r..].fill(0);
    }
}

pub fn count_leading_zeroes_skip_first(m: &[Word]) -> usize {
    let mut iter = m.iter().rev();
    let mut w;
    let mut ret = 0;

    if let Some(v) = iter.next() {
        w = *v & (WORD_MAX >> 1);

        while w == 0 {
            ret += WORD_BIT_SIZE;

            w = match iter.next() {
                Some(v) => *v,
                None => break,
            }
        }

        if w != 0 {
            while w & WORD_SIGNIFICANT_BIT == 0 {
                w <<= 1;
                ret += 1;
            }
        }
    }

    ret
}

pub fn count_leading_ones(m: &[Word]) -> usize {
    let mut ret = 0;

    for &v in m.iter().rev() {
        if v == WORD_MAX {
            ret += WORD_BIT_SIZE;
        } else {
            let mut v = v;

            while v & WORD_SIGNIFICANT_BIT != 0 {
                v <<= 1;
                ret += 1;
            }

            break;
        }
    }

    ret
}

/// Round precision to word bounday.
pub fn round_p(p: usize) -> usize {
    ((p.saturating_add(WORD_BIT_SIZE - 1)) / WORD_BIT_SIZE) * WORD_BIT_SIZE
}

// Convert rounding mode for an opposite sign.
pub fn invert_rm_for_sign(rm: RoundingMode) -> RoundingMode {
    if rm == RoundingMode::Up {
        RoundingMode::Down
    } else if rm == RoundingMode::Down {
        RoundingMode::Up
    } else {
        rm
    }
}

pub fn find_one_from(slice: &[Word], start_pos: usize) -> Option<usize> {
    let start_idx = start_pos / WORD_BIT_SIZE;
    if start_idx >= slice.len() {
        None
    } else {
        let mut iter = slice.iter().rev().skip(start_idx);
        if let Some(v) = iter.next() {
            let mut d = *v;

            let start_bit = start_pos % WORD_BIT_SIZE;
            let mut shift = start_pos;

            d <<= start_bit;

            if d != 0 {
                while d & WORD_SIGNIFICANT_BIT == 0 {
                    d <<= 1;
                    shift += 1;
                }

                return Some(shift);
            }

            shift += WORD_BIT_SIZE - start_bit;

            for v in iter {
                d = *v;

                if d != 0 {
                    while d & WORD_SIGNIFICANT_BIT == 0 {
                        d <<= 1;
                        shift += 1;
                    }

                    return Some(shift);
                }

                shift += WORD_BIT_SIZE;
            }
        }

        None
    }
}

/// Returns random subnormal number.
#[cfg(test)]
pub(crate) fn random_subnormal(p: usize) -> ExactNumNumber {
    let p = round_p(if p < 3 * WORD_BIT_SIZE { 3 * WORD_BIT_SIZE } else { p });
    let n = p - crate::common::test_rng::random::<usize>() % (2 * WORD_BIT_SIZE) - 1;
    let mut m = Vec::with_capacity(p / WORD_BIT_SIZE);

    for _ in 0..n / WORD_BIT_SIZE {
        m.push(crate::common::test_rng::random::<Word>());
    }

    if n % WORD_BIT_SIZE > 0 {
        let w = (crate::common::test_rng::random::<Word>() | WORD_SIGNIFICANT_BIT)
            >> (WORD_BIT_SIZE - n % WORD_BIT_SIZE);
        m.push(w);
    } else {
        *m.last_mut().unwrap() |= WORD_SIGNIFICANT_BIT;
    }

    m.resize(p / WORD_BIT_SIZE, 0);

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

    ExactNumNumber::from_raw_parts(&m, n, s, EXPONENT_MIN, false).unwrap()
}

#[cfg(test)]
#[inline]
pub fn rand_p() -> usize {
    crate::common::test_rng::random::<usize>() % 1000 + crate::defs::DEFAULT_P
}