use crate::common::buf::WordBuf;
use crate::common::int::SliceWithSign;
use crate::common::util::find_one_from;
use crate::common::util::shift_slice_left;
use crate::common::util::shift_slice_right;
use crate::defs::Error;
use crate::defs::WORD_BIT_SIZE;
use crate::mantissa::util::root_estimate;
use crate::mantissa::util::RightShiftedSlice;
use crate::mantissa::Mantissa;
use crate::Word;
impl Mantissa {
fn crbt_normalize_div(
m1: &mut WordBuf,
m2: &mut [Word],
m_shift: usize,
) -> Result<usize, Error> {
if let Some(shift) = find_one_from(m2, 0) {
shift_slice_left(m2, shift);
if m_shift > shift {
shift_slice_right(m1 as &mut [Word], m_shift - shift);
} else if m_shift < shift {
shift_slice_left(m1 as &mut [Word], shift - m_shift);
}
Ok(shift)
} else {
Ok(m_shift)
}
}
pub fn cbrt_rem(mut m: WordBuf) -> Result<(WordBuf, WordBuf), Error> {
let mut sbuf = root_estimate(&m, 3)?;
let mut tmpbuf = WordBuf::new(sbuf.len() * 5)?;
let mut wb = WordBuf::new(sbuf.len() + 1)?;
let mut m_shift = 0;
m.try_extend_2((m.len() + 2) * WORD_BIT_SIZE)?;
loop {
let sqbuf = &mut tmpbuf[..sbuf.len() * 2];
Self::mul_unbalanced(&sbuf, &sbuf, sqbuf)?;
m_shift = Self::crbt_normalize_div(&mut m, sqbuf, m_shift)?;
let (mut qbuf, _rbuf) = Self::div_unbalanced(&m, sqbuf)?;
wb.iter_mut()
.zip(RightShiftedSlice::new(&sbuf, WORD_BIT_SIZE - 1, 0, 1))
.for_each(|(u, v)| *u = v);
qbuf.try_extend_2((qbuf.len().max(wb.len()) + 1) * WORD_BIT_SIZE)?;
let s = SliceWithSign::new(&sbuf, 1);
let w = SliceWithSign::new(&wb, 1);
let mut q = SliceWithSign::new_mut(&mut qbuf, 1);
q.add_assign(&w);
q.div_by_word(3);
if q.cmp(&s) >= 0 {
let mut rbuf = WordBuf::new(m.len() - 2)?;
shift_slice_right(&mut m, m_shift);
rbuf.copy_from_slice(&m[..m.len() - 2]);
let mut r = SliceWithSign::new_mut(&mut rbuf, 1);
Self::cbrt_remainder(&mut tmpbuf, &mut sbuf, &mut r)?;
debug_assert!(r.sign() >= 0);
break Ok((sbuf, rbuf));
}
qbuf.trunc_leading_zeroes();
sbuf = qbuf;
}
}
fn cbrt_remainder(
tmpbuf: &mut WordBuf,
sbuf: &mut WordBuf,
r: &mut SliceWithSign,
) -> Result<(), Error> {
let (sqbuf, rest) = tmpbuf.split_at_mut(sbuf.len() * 2);
let scbuf = &mut rest[..sbuf.len() * 3];
Self::mul_unbalanced(sbuf, sbuf, sqbuf)?;
Self::mul_unbalanced(sqbuf, sbuf, scbuf)?;
let sc = SliceWithSign::new(scbuf, 1);
r.sub_assign(&sc);
Ok(())
}
}
#[cfg(test)]
mod tests {
use crate::{
defs::{WORD_MAX, WORD_SIGNIFICANT_BIT},
Word,
};
use super::*;
use rand::random;
macro_rules! assert_sqrt {
($s1:expr, $qb:expr, $rb:expr, $MAX_BUF:ident, $op:literal) => {
let mut wb = [0; MAX_BUF * 3];
let mut buf = [0; MAX_BUF * 3];
let q = SliceWithSign::new($qb, 1);
let r = SliceWithSign::new($rb, 1);
buf[q.len()..].fill(0);
buf[..q.len()].copy_from_slice(&q);
let mut qq = SliceWithSign::new_mut(&mut buf, 1);
qq.mul_assign(&q, &mut wb);
qq.mul_assign(&q, &mut wb);
qq.add_assign(&r);
assert_eq!(&qq[..$s1.len()], $s1, "{}", $op);
};
}
fn wordbuf_from_words(s: &[Word]) -> WordBuf {
let mut ret = WordBuf::new(s.len()).unwrap();
ret.copy_from_slice(s);
ret
}
#[test]
fn test_cbrt_rem() {
const MAX_BUF: usize = 100;
for s1 in [
&[WORD_MAX] as &[Word],
&[WORD_MAX, WORD_MAX],
&[WORD_MAX, WORD_MAX, WORD_MAX],
&[WORD_MAX, WORD_MAX, WORD_MAX, WORD_MAX],
&[WORD_MAX, WORD_MAX, WORD_MAX, WORD_MAX, WORD_MAX],
&[WORD_MAX, WORD_MAX, WORD_MAX, WORD_MAX, WORD_MAX, WORD_MAX],
&[WORD_MAX, WORD_MAX, WORD_MAX, WORD_MAX, WORD_MAX, WORD_MAX, WORD_MAX],
&[WORD_MAX, WORD_MAX, WORD_MAX, WORD_MAX, WORD_MAX, WORD_MAX, WORD_MAX, WORD_MAX],
] {
let (qb, rb) = Mantissa::cbrt_rem(wordbuf_from_words(s1)).unwrap();
assert_sqrt!(s1 as &[Word], &qb, &rb, MAX_BUF, "max");
}
for s1 in [
&[WORD_SIGNIFICANT_BIT] as &[Word],
&[0, WORD_SIGNIFICANT_BIT],
&[0, 0, WORD_SIGNIFICANT_BIT],
&[0, 0, 0, WORD_SIGNIFICANT_BIT],
&[0, 0, 0, 0, WORD_SIGNIFICANT_BIT],
&[0, 0, 0, 0, 0, WORD_SIGNIFICANT_BIT],
&[0, 0, 0, 0, 0, 0, WORD_SIGNIFICANT_BIT],
&[0, 0, 0, 0, 0, 0, 0, WORD_SIGNIFICANT_BIT],
] {
let (qb, rb) = Mantissa::cbrt_rem(wordbuf_from_words(s1)).unwrap();
assert_sqrt!(s1 as &[Word], &qb, &rb, MAX_BUF, "significant");
}
for s1 in [
&[WORD_SIGNIFICANT_BIT + 1] as &[Word],
&[1, WORD_SIGNIFICANT_BIT],
&[1, 0, WORD_SIGNIFICANT_BIT],
&[1, 0, 0, WORD_SIGNIFICANT_BIT],
&[1, 0, 0, 0, WORD_SIGNIFICANT_BIT],
&[1, 0, 0, 0, 0, WORD_SIGNIFICANT_BIT],
&[1, 0, 0, 0, 0, 0, WORD_SIGNIFICANT_BIT],
&[1, 0, 0, 0, 0, 0, 0, WORD_SIGNIFICANT_BIT],
] {
let (qb, rb) = Mantissa::cbrt_rem(wordbuf_from_words(s1)).unwrap();
assert_sqrt!(s1 as &[Word], &qb, &rb, MAX_BUF, "lsb + msb");
}
for s1 in [
&[1 as Word] as &[Word],
&[1, 2],
&[1, 0, 3],
&[1, 0, 0, 4],
&[1, 0, 0, 0, 5],
&[1, 0, 0, 0, 0, 6],
&[1, 0, 0, 0, 0, 0, 7],
&[1, 0, 0, 0, 0, 0, 0, 8],
] {
let (qb, rb) = Mantissa::cbrt_rem(wordbuf_from_words(s1)).unwrap();
assert_sqrt!(s1 as &[Word], &qb, &rb, MAX_BUF, "subnormal");
}
for _ in 0..1000 {
let s1 = random_normalized_slice(1, MAX_BUF);
let (qb, rb) = Mantissa::cbrt_rem(wordbuf_from_words(&s1)).unwrap();
assert_sqrt!(&s1 as &[Word], &qb, &rb, MAX_BUF, "rand");
}
}
fn random_normalized_slice(min_len: usize, max_len: usize) -> WordBuf {
let l = if max_len > min_len {
random::<usize>() % (max_len - min_len) + min_len
} else {
min_len
};
let mut s1 = WordBuf::new(l).unwrap();
for v in s1.iter_mut() {
*v = random();
}
let l = s1.len();
s1[l - 1] |= WORD_SIGNIFICANT_BIT;
s1
}
#[ignore]
#[test]
#[cfg(feature = "std")]
fn cbrt_perf() {
for _ in 0..5 {
let sz = 34;
let mut n = vec![];
let l = 10000;
for _ in 0..l {
let v = random_normalized_slice(sz, sz);
n.push(v);
}
let start_time = std::time::Instant::now();
for ni in n.drain(..) {
Mantissa::cbrt_rem(ni).unwrap();
}
let time = start_time.elapsed();
println!("cbrt {}", time.as_millis());
}
}
}