use crate::algos::exp::exp_generic as eg;
use crate::algos::support::fixed::Fixed;
use crate::algos::support::narrow_ziv::{self, WZiv};
use crate::int::types::Int;
use crate::support::rounding::{RoundingMode, is_nearest_mode};
pub(crate) const STRICT_GUARD: u32 = 30;
#[inline]
fn fixed_from_int256(raw: Int<4>) -> Fixed {
let words = raw.limbs_le();
Fixed {
negative: false,
mag: [
(words[0] as u128) | ((words[1] as u128) << 64),
(words[2] as u128) | ((words[3] as u128) << 64),
],
}
}
pub(crate) fn wide_ln2(w: u32) -> Fixed {
debug_assert!(w <= 75, "wide_ln2: working scale {w} exceeds Fixed capacity");
fixed_from_int256(crate::consts::ln2_const_n::<4>(
w,
crate::support::rounding::RoundingMode::HalfToEven,
))
}
pub(crate) fn wide_ln10(w: u32) -> Fixed {
debug_assert!(w <= 75, "wide_ln10: working scale {w} exceeds Fixed capacity");
fixed_from_int256(crate::consts::ln10_const_n::<4>(
w,
crate::support::rounding::RoundingMode::HalfToEven,
))
}
pub(crate) fn ln_fixed(v_w: Fixed, w: u32) -> Fixed {
let one_w = Fixed {
negative: false,
mag: Fixed::pow10(w),
};
let two_w = one_w.double();
let mut k: i32 = v_w.bit_length() as i32 - one_w.bit_length() as i32;
let m_w = loop {
let m = if k >= 0 {
v_w.shr(k as u32)
} else {
v_w.shl((-k) as u32)
};
if m.ge_mag(two_w) {
k += 1;
} else if !m.ge_mag(one_w) {
k -= 1;
} else {
break m;
}
};
let t = m_w.sub(one_w).div(m_w.add(one_w), w);
let t2 = t.mul(t, w);
let mut sum = t;
let mut term = t;
let mut j: u128 = 1;
loop {
term = term.mul(t2, w);
let contrib = term.div_small(2 * j + 1);
if contrib.is_zero() {
break;
}
sum = sum.add(contrib);
j += 1;
if j > 400 {
break;
}
}
let ln_m = sum.double();
let ln2 = wide_ln2(w);
let k_ln2 = if k >= 0 {
ln2.mul_u128(k as u128)
} else {
ln2.mul_u128((-k) as u128).neg()
};
k_ln2.add(ln_m)
}
#[inline]
#[must_use]
pub(crate) fn ln_with(
raw: Int<2>,
scale: u32,
working_digits: u32,
mode: RoundingMode,
) -> Option<Int<2>> {
ln_with_raw(raw.as_i128(), scale, working_digits, mode).map(Int::<2>::from_i128)
}
#[inline]
fn ln_with_raw(raw: i128, scale: u32, working_digits: u32, mode: RoundingMode) -> Option<i128> {
assert!(raw > 0, "ln kernel: argument must be positive");
let one_bits: i128 = 10_i128.pow(scale);
if raw == one_bits {
return Some(0);
}
let delta = raw - one_bits;
let ln1p_band: i128 = 10_i128.pow(scale.saturating_sub(1) / 2);
if delta.abs() <= ln1p_band && is_nearest_mode(mode) {
return Some(delta);
}
let w = scale + working_digits;
let v_w = Fixed::from_u128_mag(raw as u128, false).mul_u128(10u128.pow(working_digits));
ln_fixed(v_w, w).round_to_i128_with(w, scale, mode)
}
#[inline]
#[must_use]
pub(crate) fn ln_strict<const SCALE: u32>(raw: Int<2>, mode: RoundingMode) -> Option<Int<2>> {
ln_strict_raw::<SCALE>(raw.as_i128(), mode).map(Int::<2>::from_i128)
}
#[inline]
fn ln_strict_raw<const SCALE: u32>(raw: i128, mode: RoundingMode) -> Option<i128> {
assert!(raw > 0, "ln kernel: argument must be positive");
let one_bits: i128 = 10_i128.pow(SCALE);
if raw == one_bits {
return Some(0);
}
let delta = raw - one_bits;
let ln1p_band: i128 = 10_i128.pow(SCALE.saturating_sub(1) / 2);
if delta.abs() <= ln1p_band && is_nearest_mode(mode) {
return Some(delta);
}
let w = SCALE + STRICT_GUARD;
let v_w = Fixed::from_u128_mag(raw as u128, false).mul_u128(10u128.pow(STRICT_GUARD));
let lv = ln_fixed(v_w, w);
match lv.round_to_i128_clear_of_tie(w, SCALE, mode) {
Some(r) => r,
None => narrow_ziv::walk_checked(
lv.round_to_i128_with(w, SCALE, mode),
STRICT_GUARD,
SCALE,
mode,
|g| ln_ziv(raw, SCALE, g),
),
}
}
fn ln_ziv(raw: i128, scale: u32, g: u32) -> WZiv {
let w = scale + g;
eg::ln_fixed::<WZiv>(narrow_ziv::lift(raw, g), w, narrow_ziv::ln2_w(w))
}
fn log_ratio_ziv(x_raw: i128, b_sel: i128, scale: u32, g: u32) -> WZiv {
let w = scale + g;
let ln2 = narrow_ziv::ln2_w(w);
let lx = eg::ln_fixed::<WZiv>(narrow_ziv::lift(x_raw, g), w, ln2);
let lb = match b_sel {
-2 => ln2,
-10 => narrow_ziv::ln10_w(w),
b_raw => eg::ln_fixed::<WZiv>(narrow_ziv::lift(b_raw, g), w, ln2),
};
eg::div::<WZiv>(lx, lb, w)
}
pub(crate) fn gcd_u128(mut a: u128, mut b: u128) -> u128 {
while b != 0 {
let t = a % b;
a = b;
b = t;
}
a
}
pub(crate) fn pow_bounded(base: u128, e: u128) -> Option<WZiv> {
if base == 0 {
return None; }
let bl = (128 - base.leading_zeros()) as u128;
if bl.saturating_mul(e) > 1500 {
return None;
}
let mut acc = WZiv::from_i128(1);
let mut b = WZiv::from_u128(base);
let mut k = e;
while k > 0 {
if k & 1 == 1 {
acc = acc * b;
}
k >>= 1;
if k > 0 {
b = b * b;
}
}
Some(acc)
}
fn log_rational_pow_pin(
xn: u128,
xd: u128,
bn: u128,
bd: u128,
scale: u32,
num2: i128,
mode: RoundingMode,
) -> Option<i128> {
if num2 == 0 {
return None; }
let den: u128 = 2 * 10u128.pow(scale);
let neg = num2 < 0;
let num = num2.unsigned_abs();
let g = gcd_u128(num, den);
let p = num / g;
let q = den / g;
let (tn, td) = if neg { (bd, bn) } else { (bn, bd) };
let lx_n = pow_bounded(xn, q)?;
let lx_d = pow_bounded(xd, q)?;
let rb_n = pow_bounded(tn, p)?;
let rb_d = pow_bounded(td, p)?;
if lx_n != rb_n || lx_d != rb_d {
return None;
}
let q_mag = num / 2;
if num & 1 == 0 {
return Some(if neg { -(q_mag as i128) } else { q_mag as i128 });
}
let bump = crate::support::rounding::should_bump(
mode,
core::cmp::Ordering::Equal,
q_mag & 1 == 1,
!neg,
);
let m = (q_mag + u128::from(bump)) as i128;
Some(if neg { -m } else { m })
}
pub(crate) fn reduce_fraction(a: u128, b: u128) -> (u128, u128) {
let g = gcd_u128(a, b);
(a / g, b / g)
}
#[inline]
fn log_exact_int_pin(value_raw: i128, base_int: i128, scale: u32, k: i128) -> Option<i128> {
let one_s = 10i128.checked_pow(scale)?;
if k == 0 {
return (value_raw == one_s).then_some(0);
}
if base_int == 0 {
return None;
}
let kk = k.unsigned_abs();
let exact = if k > 0 {
if value_raw % one_s != 0 {
return None;
}
let value_int = value_raw / one_s;
let mut pow: i128 = 1;
let mut ok = true;
for _ in 0..kk {
match pow.checked_mul(base_int) {
Some(p) => pow = p,
None => {
ok = false;
break;
}
}
}
ok && pow == value_int
} else {
let mut cur = value_raw;
let mut ok = true;
for _ in 0..kk {
match cur.checked_mul(base_int) {
Some(p) => cur = p,
None => {
ok = false;
break;
}
}
}
ok && cur == one_s
};
if exact { k.checked_mul(one_s) } else { None }
}
#[inline]
#[must_use]
pub(crate) fn log_with(
raw: Int<2>,
base_raw: Int<2>,
scale: u32,
working_digits: u32,
mode: RoundingMode,
) -> Option<Int<2>> {
log_with_raw(
raw.as_i128(),
base_raw.as_i128(),
scale,
working_digits,
mode,
)
.map(Int::<2>::from_i128)
}
#[inline]
fn log_with_raw(
raw: i128,
base_raw: i128,
scale: u32,
working_digits: u32,
mode: RoundingMode,
) -> Option<i128> {
assert!(raw > 0, "D38::log: argument must be positive");
assert!(base_raw > 0, "D38::log: base must be positive");
let w = scale + working_digits;
let pow = 10u128.pow(working_digits);
let v_w = Fixed::from_u128_mag(raw as u128, false).mul_u128(pow);
let b_w = Fixed::from_u128_mag(base_raw as u128, false).mul_u128(pow);
let ln_b = ln_fixed(b_w, w);
assert!(
!ln_b.is_zero(),
"D38::log: base must not equal 1 (ln(1) is zero)"
);
let ratio = ln_fixed(v_w, w).div(ln_b, w);
let k = ratio.round_to_nearest_int(w);
let base_int = match 10i128.checked_pow(scale) {
Some(one_s) if base_raw % one_s == 0 => base_raw / one_s,
_ => 0,
};
if let Some(pinned) = log_exact_int_pin(raw, base_int, scale, k) {
return Some(pinned);
}
ratio.round_to_i128_with(w, scale, mode)
}
#[inline]
#[must_use]
pub(crate) fn log_strict<const SCALE: u32>(raw: Int<2>, base_raw: Int<2>, mode: RoundingMode) -> Option<Int<2>> {
log_strict_raw(raw.as_i128(), base_raw.as_i128(), SCALE, mode).map(Int::<2>::from_i128)
}
fn log_strict_raw(raw: i128, base_raw: i128, scale: u32, mode: RoundingMode) -> Option<i128> {
assert!(raw > 0, "D38::log: argument must be positive");
assert!(base_raw > 0, "D38::log: base must be positive");
let w = scale + STRICT_GUARD;
let pow = 10u128.pow(STRICT_GUARD);
let v_w = Fixed::from_u128_mag(raw as u128, false).mul_u128(pow);
let b_w = Fixed::from_u128_mag(base_raw as u128, false).mul_u128(pow);
let ln_b = ln_fixed(b_w, w);
assert!(
!ln_b.is_zero(),
"D38::log: base must not equal 1 (ln(1) is zero)"
);
let ratio = ln_fixed(v_w, w).div(ln_b, w);
let k = ratio.round_to_nearest_int(w);
let base_int = match 10i128.checked_pow(scale) {
Some(one_s) if base_raw % one_s == 0 => base_raw / one_s,
_ => 0,
};
if let Some(pinned) = log_exact_int_pin(raw, base_int, scale, k) {
return Some(pinned);
}
match ratio.round_to_i128_clear_of_tie(w, scale, mode) {
Some(r) => r,
None => {
if let Some(num2) = ratio.double().round_to_i128_with(w, scale, RoundingMode::HalfToEven)
{
let one_s = 10u128.pow(scale);
let (xn, xd) = reduce_fraction(raw as u128, one_s);
let (bn, bd) = reduce_fraction(base_raw as u128, one_s);
if let Some(pinned) =
log_rational_pow_pin(xn, xd, bn, bd, scale, num2, mode)
{
return Some(pinned);
}
}
narrow_ziv::walk_checked(
ratio.round_to_i128_with(w, scale, mode),
STRICT_GUARD,
scale,
mode,
|g| log_ratio_ziv(raw, base_raw, scale, g),
)
}
}
}
#[inline]
#[must_use]
pub(crate) fn log2_with(raw: Int<2>, scale: u32, working_digits: u32, mode: RoundingMode) -> Option<Int<2>> {
log2_with_raw(raw.as_i128(), scale, working_digits, mode).map(Int::<2>::from_i128)
}
#[inline]
fn log2_with_raw(raw: i128, scale: u32, working_digits: u32, mode: RoundingMode) -> Option<i128> {
assert!(raw > 0, "D38::log2: argument must be positive");
let w = scale + working_digits;
let v_w = Fixed::from_u128_mag(raw as u128, false).mul_u128(10u128.pow(working_digits));
let ratio = ln_fixed(v_w, w).div(wide_ln2(w), w);
let k = ratio.round_to_nearest_int(w);
if let Some(pinned) = log_exact_int_pin(raw, 2, scale, k) {
return Some(pinned);
}
ratio.round_to_i128_with(w, scale, mode)
}
#[inline]
#[must_use]
pub(crate) fn log2_strict<const SCALE: u32>(raw: Int<2>, mode: RoundingMode) -> Option<Int<2>> {
log2_strict_raw(raw.as_i128(), SCALE, mode).map(Int::<2>::from_i128)
}
fn log2_strict_raw(raw: i128, scale: u32, mode: RoundingMode) -> Option<i128> {
assert!(raw > 0, "D38::log2: argument must be positive");
let w = scale + STRICT_GUARD;
let v_w = Fixed::from_u128_mag(raw as u128, false).mul_u128(10u128.pow(STRICT_GUARD));
let ratio = ln_fixed(v_w, w).div(wide_ln2(w), w);
let k = ratio.round_to_nearest_int(w);
if let Some(pinned) = log_exact_int_pin(raw, 2, scale, k) {
return Some(pinned);
}
match ratio.round_to_i128_clear_of_tie(w, scale, mode) {
Some(r) => r,
None => narrow_ziv::walk_checked(
ratio.round_to_i128_with(w, scale, mode),
STRICT_GUARD,
scale,
mode,
|g| log_ratio_ziv(raw, -2, scale, g),
),
}
}
#[inline]
#[must_use]
pub(crate) fn log10_with(raw: Int<2>, scale: u32, working_digits: u32, mode: RoundingMode) -> Option<Int<2>> {
log10_with_raw(raw.as_i128(), scale, working_digits, mode).map(Int::<2>::from_i128)
}
#[inline]
fn log10_with_raw(raw: i128, scale: u32, working_digits: u32, mode: RoundingMode) -> Option<i128> {
assert!(raw > 0, "D38::log10: argument must be positive");
let w = scale + working_digits;
let v_w = Fixed::from_u128_mag(raw as u128, false).mul_u128(10u128.pow(working_digits));
let ratio = ln_fixed(v_w, w).div(wide_ln10(w), w);
let k = ratio.round_to_nearest_int(w);
if let Some(pinned) = log_exact_int_pin(raw, 10, scale, k) {
return Some(pinned);
}
ratio.round_to_i128_with(w, scale, mode)
}
#[inline]
#[must_use]
pub(crate) fn log10_strict<const SCALE: u32>(raw: Int<2>, mode: RoundingMode) -> Option<Int<2>> {
log10_strict_raw(raw.as_i128(), SCALE, mode).map(Int::<2>::from_i128)
}
fn log10_strict_raw(raw: i128, scale: u32, mode: RoundingMode) -> Option<i128> {
assert!(raw > 0, "D38::log10: argument must be positive");
let w = scale + STRICT_GUARD;
let v_w = Fixed::from_u128_mag(raw as u128, false).mul_u128(10u128.pow(STRICT_GUARD));
let ratio = ln_fixed(v_w, w).div(wide_ln10(w), w);
let k = ratio.round_to_nearest_int(w);
if let Some(pinned) = log_exact_int_pin(raw, 10, scale, k) {
return Some(pinned);
}
match ratio.round_to_i128_clear_of_tie(w, scale, mode) {
Some(r) => r,
None => narrow_ziv::walk_checked(
ratio.round_to_i128_with(w, scale, mode),
STRICT_GUARD,
scale,
mode,
|g| log_ratio_ziv(raw, -10, scale, g),
),
}
}
#[cfg(test)]
mod near_tie_pins {
use super::*;
#[test]
fn ln_directed_near_one_deep_d38_s38() {
let raw = Int::<2>::from_i128(10_i128.pow(38) + 1);
assert_eq!(
ln_strict::<38>(raw, RoundingMode::Floor).map(|v| v.as_i128()),
Some(0),
"ln Floor"
);
assert_eq!(
ln_strict::<38>(raw, RoundingMode::Trunc).map(|v| v.as_i128()),
Some(0),
"ln Trunc"
);
assert_eq!(
ln_strict::<38>(raw, RoundingMode::Ceiling).map(|v| v.as_i128()),
Some(1),
"ln Ceiling"
);
let raw_lo = Int::<2>::from_i128(10_i128.pow(38) - 1);
assert_eq!(
ln_strict::<38>(raw_lo, RoundingMode::Floor).map(|v| v.as_i128()),
Some(-2),
"ln(1−ulp) Floor"
);
assert_eq!(
ln_strict::<38>(raw_lo, RoundingMode::Trunc).map(|v| v.as_i128()),
Some(-1),
"ln(1−ulp) Trunc"
);
}
#[test]
fn ln_directed_deep_band_sweep() {
let mut bad = 0u32;
for s in 31u32..=38 {
let one: i128 = 10_i128.pow(s);
let max_d = ((2.0 * 10f64.powi(s as i32 - 30)).sqrt() as i128).max(1);
let mut d: i128 = 1;
while d <= max_d {
let raw = Int::<2>::from_i128(one + d);
let fl = ln_strict_dyn(s, raw, RoundingMode::Floor);
let ce = ln_strict_dyn(s, raw, RoundingMode::Ceiling);
if fl != Some(d - 1) || ce != Some(d) {
println!("LIVE ln(1+{d}ulp) S={s}: Floor={fl:?} (want {}), Ceil={ce:?} (want {d})", d-1);
bad += 1;
}
let raw_lo = Int::<2>::from_i128(one - d);
let fl2 = ln_strict_dyn(s, raw_lo, RoundingMode::Floor);
let tr2 = ln_strict_dyn(s, raw_lo, RoundingMode::Trunc);
if fl2 != Some(-d - 1) || tr2 != Some(-d) {
println!("LIVE ln(1-{d}ulp) S={s}: Floor={fl2:?} (want {}), Trunc={tr2:?} (want {})", -d-1, -d);
bad += 1;
}
d = if d < 20 { d + 1 } else { d + d / 3 };
}
}
assert_eq!(bad, 0, "ln deep directed band has live misrounds");
}
fn ln_strict_dyn(s: u32, raw: Int<2>, mode: RoundingMode) -> Option<i128> {
let v = match s {
31 => ln_strict::<31>(raw, mode),
32 => ln_strict::<32>(raw, mode),
33 => ln_strict::<33>(raw, mode),
34 => ln_strict::<34>(raw, mode),
35 => ln_strict::<35>(raw, mode),
36 => ln_strict::<36>(raw, mode),
37 => ln_strict::<37>(raw, mode),
38 => ln_strict::<38>(raw, mode),
_ => unreachable!(),
};
v.map(|v| v.as_i128())
}
#[test]
fn log_fractional_base_high_scale_budget_product_does_not_wrap() {
let xr: i128 = 1096614350149675660535769418_i128 * 1_000_000_000;
let br: i128 = 182017066872921546105935121_i128 * 1_000_000_000;
for mode in [
RoundingMode::HalfToEven,
RoundingMode::Floor,
RoundingMode::Ceiling,
] {
let r = log_strict_raw(xr, br, 37, mode);
assert!(r.is_some(), "log fractional-base s37 mode={mode:?}");
}
}
#[test]
fn log_rational_exact_d38() {
let one19 = 10_i128.pow(19);
let x = Int::<2>::from_i128(8 * one19);
let b = Int::<2>::from_i128(4 * one19);
for mode in [
RoundingMode::HalfToEven,
RoundingMode::HalfAwayFromZero,
RoundingMode::HalfTowardZero,
RoundingMode::Floor,
RoundingMode::Ceiling,
RoundingMode::Trunc,
] {
assert_eq!(
log_strict_raw(x.as_i128(), b.as_i128(), 19, mode),
Some(15 * one19 / 10),
"log_4(8) mode={mode:?}"
);
}
}
#[test]
fn log_rational_half_tie_s0() {
let x = Int::<2>::from_i128(32);
let b = Int::<2>::from_i128(4);
assert_eq!(
log_strict_raw(x.as_i128(), b.as_i128(), 0, RoundingMode::HalfToEven),
Some(2),
"log_4(32) HalfToEven"
);
assert_eq!(
log_strict_raw(x.as_i128(), b.as_i128(), 0, RoundingMode::HalfAwayFromZero),
Some(3),
"log_4(32) HalfAwayFromZero"
);
assert_eq!(
log_strict_raw(x.as_i128(), b.as_i128(), 0, RoundingMode::HalfTowardZero),
Some(2),
"log_4(32) HalfTowardZero"
);
let x8 = Int::<2>::from_i128(80);
let b16 = Int::<2>::from_i128(160);
assert_eq!(
log_strict_raw(x8.as_i128(), b16.as_i128(), 1, RoundingMode::HalfToEven),
Some(8),
"log_16(8) s1 HalfToEven"
);
assert_eq!(
log_strict_raw(x8.as_i128(), b16.as_i128(), 1, RoundingMode::HalfTowardZero),
Some(7),
"log_16(8) s1 HalfTowardZero"
);
}
}