use num_bigint::BigInt;
use num_rational::Ratio;
use num_traits::Zero;
use rustc_hash::FxHashMap;
use astro_float::{BigFloat, Consts, Radix, RoundingMode, Sign, WORD_BIT_SIZE};
use num_complex::Complex64;
use crate::base::arena::Arena;
use crate::base::errors::SymplexError;
use crate::base::libfn::LibFn;
use crate::base::node::{ExprId, ExprNode};
use crate::base::walk;
use tracing::debug;
use crate::base::bigcomplex::Complex;
pub(crate) const QUADRATURE_MAX_DIGITS: u32 = 16;
pub(crate) fn evalf(arena: &Arena, expr: ExprId, digits: u32) -> Result<String, SymplexError> {
let binary_prec = (digits as usize) * 34 / 10 + 64;
let prec = binary_prec.max(128);
debug!(
digits = digits,
binary_precision = prec,
"evalf: starting numerical evaluation"
);
let max_prec = arena.config.max_evalf_precision as usize;
if prec > max_prec {
return Err(SymplexError::PrecisionExhausted {
requested: digits,
achieved: (max_prec * 10 / 34).saturating_sub(6) as u32,
});
}
let free = walk::free_symbols(arena, expr);
if !free.is_empty() {
let mut names: Vec<&str> = free
.iter()
.filter_map(|&id| match arena.node(id) {
ExprNode::Symbol(sid) => Some(arena.symbol_name(*sid)),
_ => None,
})
.collect();
names.sort_unstable();
if let Some(name) = names.first() {
debug!(symbol = name, "evalf: expression has free symbols");
return Err(SymplexError::FreeSymbol {
name: (*name).to_owned(),
});
}
}
if digits > QUADRATURE_MAX_DIGITS && contains_definite_integral(arena, expr) {
return Err(SymplexError::NotImplemented(format!(
"arbitrary-precision evaluation of a definite integral is not implemented \
({digits} digits requested, at most {QUADRATURE_MAX_DIGITS} available from \
f64 Gauss–Kronrod quadrature); use eval_f64 or integrate_definite"
)));
}
let rm = RoundingMode::ToEven;
let mut cc = Consts::new().map_err(|e| {
SymplexError::NotImplemented(format!("astro-float constants init failed: {e:?}"))
})?;
let post_order = walk::post_order_ids(arena, expr);
let mut cache: FxHashMap<ExprId, Complex> = FxHashMap::default();
for &id in &post_order {
match eval_node(arena, id, &cache, prec, rm, &mut cc) {
Ok(value) => {
cache.insert(id, value);
}
Err(e) => {
if id == expr {
return Err(e);
}
}
}
}
let result = cache.get(&expr).ok_or_else(|| {
SymplexError::NotImplemented("evalf: expression not found in cache".into())
})?;
if result.0.is_nan() || result.1.is_nan() {
return Err(SymplexError::PrecisionExhausted {
requested: digits,
achieved: 0,
});
}
format_complex(result, digits, prec, rm, &mut cc)
}
pub(crate) const F64_DIGITS: u32 = 16;
thread_local! {
static F64_CONSTS: std::cell::RefCell<Option<Consts>> = const { std::cell::RefCell::new(None) };
}
fn with_f64_consts<T>(
f: impl FnOnce(&mut Consts) -> Result<T, SymplexError>,
) -> Result<T, SymplexError> {
let mut cc = match F64_CONSTS.with(|slot| slot.borrow_mut().take()) {
Some(cc) => cc,
None => Consts::new().map_err(|e| {
SymplexError::NotImplemented(format!("astro-float constants init failed: {e:?}"))
})?,
};
let out = f(&mut cc);
F64_CONSTS.with(|slot| {
let mut slot = slot.borrow_mut();
if slot.is_none() {
*slot = Some(cc);
}
});
out
}
fn evalf_value(arena: &Arena, expr: ExprId, digits: u32) -> Result<Complex, SymplexError> {
let prec = ((digits as usize) * 34 / 10 + 64).max(128);
let max_prec = arena.config.max_evalf_precision as usize;
if prec > max_prec {
return Err(SymplexError::PrecisionExhausted {
requested: digits,
achieved: (max_prec * 10 / 34).saturating_sub(6) as u32,
});
}
let free = walk::free_symbols(arena, expr);
if !free.is_empty() {
let mut names: Vec<&str> = free
.iter()
.filter_map(|&id| match arena.node(id) {
ExprNode::Symbol(sid) => Some(arena.symbol_name(*sid)),
_ => None,
})
.collect();
names.sort_unstable();
if let Some(name) = names.first() {
debug!(symbol = name, "evalf_value: expression has free symbols");
return Err(SymplexError::FreeSymbol {
name: (*name).to_owned(),
});
}
}
let rm = RoundingMode::ToEven;
let post_order = walk::post_order_ids(arena, expr);
let result = with_f64_consts(|cc| {
let mut cache: FxHashMap<ExprId, Complex> = FxHashMap::default();
for &id in &post_order {
match eval_node(arena, id, &cache, prec, rm, cc) {
Ok(value) => {
cache.insert(id, value);
}
Err(e) if id == expr => return Err(e),
Err(_) => {}
}
}
cache.remove(&expr).ok_or_else(|| {
SymplexError::NotImplemented("evalf: expression not found in cache".into())
})
})?;
if result.0.is_nan() || result.1.is_nan() {
return Err(SymplexError::PrecisionExhausted {
requested: digits,
achieved: 0,
});
}
Ok(result)
}
fn finite_part_to_f64(part: &BigFloat, other: &BigFloat) -> Result<f64, SymplexError> {
if is_negligible_part(part, other, F64_DIGITS) {
return Ok(0.0);
}
if part.is_inf() {
return Err(SymplexError::Unevaluable {
reason: "the value overflowed the arbitrary-precision exponent range".into(),
});
}
bigfloat_to_f64_rounded(part, RoundingMode::ToEven)
}
pub(crate) fn evalf_complex64(arena: &Arena, expr: ExprId) -> Result<Complex64, SymplexError> {
let z = evalf_value(arena, expr, F64_DIGITS)?;
let re = finite_part_to_f64(&z.0, &z.1)?;
let im = finite_part_to_f64(&z.1, &z.0)?;
Ok(Complex64::new(re, im))
}
pub(crate) fn evalf_f64(arena: &Arena, expr: ExprId) -> Result<f64, SymplexError> {
let z = evalf_complex64(arena, expr)?;
if z.im.abs() > 1e-15 {
return Err(SymplexError::ComputationFailed {
operation: "eval_f64",
reason: format!(
"expression has nonzero imaginary part (im={}); use eval_complex64() for complex results",
z.im
),
});
}
Ok(z.re)
}
pub(crate) fn eval_const_f64(arena: &mut Arena, expr: ExprId) -> Option<f64> {
let evaled = crate::transforms::eval::eval(arena, expr);
if let Some(r) = arena.as_num(evaled) {
let result = crate::base::numeric::ratio_to_f64(r)?;
tracing::trace!(result, "eval_const_f64: rational fast path");
return Some(result);
}
match evalf(arena, evaled, 16) {
Ok(s) => {
let result = s.parse::<f64>().ok();
tracing::trace!(?result, decimal_str = %s, "eval_const_f64: evalf path");
result
}
Err(e) => {
tracing::debug!(?e, "eval_const_f64: evalf failed");
None
}
}
}
fn eval_node(
arena: &Arena,
id: ExprId,
cache: &FxHashMap<ExprId, Complex>,
prec: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<Complex, SymplexError> {
let node = arena.node(id);
match node {
ExprNode::Num(nid) => {
let r = arena.num(*nid);
Ok((ratio_to_bigfloat(r, prec, rm), BigFloat::new(prec)))
}
ExprNode::Pi => Ok((cc.pi(prec, rm).clone(), BigFloat::new(prec))),
ExprNode::E => Ok((cc.e(prec, rm).clone(), BigFloat::new(prec))),
ExprNode::ImaginaryUnit => Ok((BigFloat::new(prec), BigFloat::from_i32(1, prec))),
ExprNode::EulerGamma => {
debug!(prec, "evalf: EulerGamma via Brent–McMillan");
let g = arb_euler_gamma(prec, rm, cc)?;
Ok((g, BigFloat::new(prec)))
}
ExprNode::Catalan => {
debug!(prec, "evalf: Catalan via Ramanujan series");
let g = arb_catalan(prec, rm, cc)?;
Ok((g, BigFloat::new(prec)))
}
ExprNode::GoldenRatio => {
let wp = prec + 16;
let five = BigFloat::from_i32(5, wp);
let sqrt5 = five.sqrt(wp, rm);
let one = BigFloat::from_i32(1, wp);
let two = BigFloat::from_i32(2, wp);
let phi = one.add(&sqrt5, wp, rm).div(&two, prec, rm);
Ok((phi, BigFloat::new(prec)))
}
ExprNode::Re(inner) => {
let val = get_cached(cache, *inner)?;
Ok((val.0.clone(), BigFloat::new(prec)))
}
ExprNode::Im(inner) => {
let val = get_cached(cache, *inner)?;
Ok((val.1.clone(), BigFloat::new(prec)))
}
ExprNode::Conjugate(inner) => {
let val = get_cached(cache, *inner)?;
Ok((val.0.clone(), val.1.neg()))
}
ExprNode::Arg(inner) => {
let val = get_cached(cache, *inner)?;
if val.0.is_zero() && val.1.is_zero() {
return Err(SymplexError::Unevaluable {
reason: "arg(0) is undefined".into(),
});
}
Ok((atan2_bf(&val.1, &val.0, prec, rm, cc), BigFloat::new(prec)))
}
ExprNode::Si(inner) => {
let val = get_cached(cache, *inner)?;
if !val.1.is_zero() {
return Err(SymplexError::Unevaluable {
reason: "Si of complex argument not yet supported in evalf".into(),
});
}
debug!(prec, "evalf: Si via series/asymptotic");
let (si, _ci) = arb_si_ci(&val.0, true, prec, rm, cc)?;
Ok((si, BigFloat::new(prec)))
}
ExprNode::Ci(inner) => {
let val = get_cached(cache, *inner)?;
if !val.1.is_zero() {
return Err(SymplexError::Unevaluable {
reason: "Ci of complex argument not yet supported in evalf".into(),
});
}
if val.0.is_zero() {
return Err(SymplexError::Unevaluable {
reason: "Ci(0) is -∞".into(),
});
}
debug!(prec, "evalf: Ci via series/asymptotic");
let x_abs = val.0.abs();
let (_si, ci) = arb_si_ci(&x_abs, false, prec, rm, cc)?;
if val.0.is_negative() {
Ok((ci, cc.pi(prec, rm).clone()))
} else {
Ok((ci, BigFloat::new(prec)))
}
}
ExprNode::Ei(inner) => {
let val = get_cached(cache, *inner)?;
if !val.1.is_zero() {
return Err(SymplexError::Unevaluable {
reason: "Ei of complex argument not yet supported in evalf".into(),
});
}
if val.0.is_zero() {
return Err(SymplexError::Unevaluable {
reason: "Ei(0) is -∞".into(),
});
}
debug!(prec, "evalf: Ei via series/asymptotic");
let r = arb_ei(&val.0, prec, rm, cc)?;
Ok((r, BigFloat::new(prec)))
}
ExprNode::Li(inner) => {
let val = get_cached(cache, *inner)?;
if !val.1.is_zero() || val.0.is_negative() {
return Err(SymplexError::Unevaluable {
reason: "li of complex or negative argument not yet supported in evalf".into(),
});
}
if val.0.is_zero() {
return Ok(c_zero(prec));
}
let wp = prec + 32;
let ln_x = val.0.ln(wp, rm, cc);
if ln_x.is_zero() {
return Err(SymplexError::Unevaluable {
reason: "li(1) is -∞".into(),
});
}
debug!(prec, "evalf: li via Ei(ln x)");
let r = arb_ei(&ln_x, prec, rm, cc)?;
Ok((r, BigFloat::new(prec)))
}
ExprNode::Zeta(inner) => {
let val = get_cached(cache, *inner)?;
if !val.1.is_zero() {
return Err(SymplexError::Unevaluable {
reason: "zeta of complex argument not yet supported in evalf".into(),
});
}
debug!(prec, "evalf: zeta via Borwein / functional equation");
let r = arb_zeta(&val.0, prec, rm, cc)?;
Ok((r, BigFloat::new(prec)))
}
ExprNode::Polygamma(n_id, x_id) => {
let n_val = get_cached(cache, *n_id)?;
let x_val = get_cached(cache, *x_id)?;
if !n_val.1.is_zero() || !x_val.1.is_zero() {
return Err(SymplexError::Unevaluable {
reason: "polygamma of complex arguments not yet supported in evalf".into(),
});
}
let n_f = bigfloat_to_f64(&n_val.0, rm, cc)?;
let n_round = n_f.round();
if (n_f - n_round).abs() > 1e-12 || !(0.0..=10_000.0).contains(&n_round) {
return Err(SymplexError::Unevaluable {
reason: "polygamma order must be a non-negative integer".into(),
});
}
let n = n_round as u32;
if n == 0 {
let r = arb_digamma(&x_val.0, prec, rm, cc)?;
return Ok((r, BigFloat::new(prec)));
}
debug!(
prec,
n, "evalf: polygamma via recurrence + asymptotic series"
);
let r = arb_polygamma(n, &x_val.0, prec, rm, cc)?;
Ok((r, BigFloat::new(prec)))
}
ExprNode::KroneckerDelta(i_id, j_id) => {
let iv = get_cached(cache, *i_id)?;
let jv = get_cached(cache, *j_id)?;
let d_re = iv.0.sub(&jv.0, prec, rm);
let d_im = iv.1.sub(&jv.1, prec, rm);
if d_re.is_zero() && d_im.is_zero() {
Ok(c_one(prec))
} else {
Ok(c_zero(prec))
}
}
ExprNode::PhysicalConstant(_, value_id) => {
if let Some(v) = cache.get(value_id) {
return Ok(v.clone());
}
let mut sub_cache: FxHashMap<ExprId, Complex> = FxHashMap::default();
for id in walk::post_order_ids(arena, *value_id) {
if !sub_cache.contains_key(&id) {
let v = eval_node(arena, id, &sub_cache, prec, rm, cc)?;
sub_cache.insert(id, v);
}
}
sub_cache
.remove(value_id)
.ok_or_else(|| SymplexError::Unevaluable {
reason: "physical constant value could not be evaluated".into(),
})
}
ExprNode::Infinity | ExprNode::NegInfinity | ExprNode::ComplexInfinity => {
Err(SymplexError::Unevaluable {
reason: "cannot evaluate infinity to finite precision".into(),
})
}
ExprNode::Factorial(x_id) => {
let x_val = get_cached(cache, *x_id)?;
if !x_val.1.is_zero() {
return Err(SymplexError::Unevaluable {
reason: "factorial of a complex argument not yet supported in evalf".into(),
});
}
debug!(
prec,
"evalf: factorial via arbitrary-precision Gamma(x + 1)"
);
let x_plus_1 = x_val.0.add(&BigFloat::from_i32(1, prec), prec, rm);
let r = arb_gamma_real(&x_plus_1, prec, rm, cc)?;
Ok((r, BigFloat::new(prec)))
}
ExprNode::Binomial(n_id, k_id) => {
debug!(prec, "evalf: Binomial via arbitrary-precision Gamma");
let n_val = get_cached(cache, *n_id)?;
let k_val = get_cached(cache, *k_id)?;
if !n_val.1.is_zero() || !k_val.1.is_zero() {
return Err(SymplexError::Unevaluable {
reason: "Binomial of complex arguments not yet supported in evalf".into(),
});
}
if let (Some(n_i), Some(k_i)) =
(bf_as_int(&n_val.0, rm, cc)?, bf_as_int(&k_val.0, rm, cc)?)
&& n_i >= 0
&& k_i > n_i
{
return Ok((BigFloat::new(prec), BigFloat::new(prec)));
}
let one_bf = BigFloat::from_i32(1, prec);
let n_plus_1 = n_val.0.add(&one_bf, prec, rm);
let k_plus_1 = k_val.0.add(&one_bf, prec, rm);
let n_minus_k = n_val.0.sub(&k_val.0, prec, rm);
let n_minus_k_plus_1 = n_minus_k.add(&one_bf, prec, rm);
let g_numer = arb_gamma_real(&n_plus_1, prec, rm, cc)?;
let g_k = arb_gamma_real(&k_plus_1, prec, rm, cc)?;
let g_nmk = arb_gamma_real(&n_minus_k_plus_1, prec, rm, cc)?;
let g_denom = g_k.mul(&g_nmk, prec, rm);
if g_denom.is_zero() {
return Err(SymplexError::Unevaluable {
reason: "Binomial denominator is zero (pole in Gamma)".into(),
});
}
let result = g_numer.div(&g_denom, prec, rm);
Ok((result, BigFloat::new(prec)))
}
ExprNode::Gamma(inner) => {
let val = get_cached(cache, *inner)?;
if !val.1.is_zero() {
return Err(SymplexError::Unevaluable {
reason: "Gamma of complex argument not yet supported in evalf".into(),
});
}
debug!(prec, "evalf: Gamma via Stirling series");
let result = arb_gamma_real(&val.0, prec, rm, cc)?;
Ok((result, BigFloat::new(prec)))
}
ExprNode::LogGamma(inner) => {
let val = get_cached(cache, *inner)?;
if !val.1.is_zero() {
return Err(SymplexError::Unevaluable {
reason: "LogGamma of complex argument not yet supported in evalf".into(),
});
}
let x_approx = bigfloat_to_f64(&val.0, rm, cc)?;
if x_approx <= 0.0 {
let rounded = x_approx.round();
if (rounded - x_approx).abs() < 1e-12 && rounded <= 0.0 {
return Err(SymplexError::Unevaluable {
reason: "LogGamma at non-positive integer pole".into(),
});
}
let one = BigFloat::from_i32(1, prec);
let one_minus_x = one.sub(&val.0, prec, rm);
let log_gamma_1mx = stirling_log_gamma(&one_minus_x, prec, rm, cc)?;
let pi_val = cc.pi(prec, rm).clone();
let pi_x = pi_val.mul(&val.0, prec, rm);
let sin_pi_x = pi_x.sin(prec, rm, cc);
let abs_sin = sin_pi_x.abs();
if abs_sin.is_zero() {
return Err(SymplexError::Unevaluable {
reason: "LogGamma at non-positive integer pole".into(),
});
}
let ln_pi = cc.pi(prec, rm).clone().ln(prec, rm, cc);
let ln_abs_sin = abs_sin.ln(prec, rm, cc);
let result = ln_pi
.sub(&ln_abs_sin, prec, rm)
.sub(&log_gamma_1mx, prec, rm);
Ok((result, BigFloat::new(prec)))
} else {
debug!(prec, "evalf: LogGamma via Stirling series");
let result = stirling_log_gamma(&val.0, prec, rm, cc)?;
Ok((result, BigFloat::new(prec)))
}
}
ExprNode::Digamma(inner) => {
let val = get_cached(cache, *inner)?;
if !val.1.is_zero() {
return Err(SymplexError::Unevaluable {
reason: "Digamma of complex argument not yet supported in evalf".into(),
});
}
let result = arb_digamma(&val.0, prec, rm, cc)?;
Ok((result, BigFloat::new(prec)))
}
ExprNode::Erf(inner) => {
let val = get_cached(cache, *inner)?;
if !val.1.is_zero() {
return Err(SymplexError::Unevaluable {
reason: "erf of complex argument not yet supported in evalf".into(),
});
}
debug!(prec, "evalf: erf via arbitrary-precision series");
let result = arb_erf(&val.0, prec, rm, cc)?;
Ok((result, BigFloat::new(prec)))
}
ExprNode::Erfc(inner) => {
let val = get_cached(cache, *inner)?;
if !val.1.is_zero() {
return Err(SymplexError::Unevaluable {
reason: "erfc of complex argument not yet supported in evalf".into(),
});
}
debug!(prec, "evalf: erfc via arbitrary-precision series");
let result = arb_erfc(&val.0, prec, rm, cc)?;
Ok((result, BigFloat::new(prec)))
}
ExprNode::LambertW(inner) => {
let val = get_cached(cache, *inner)?;
if !val.1.is_zero() {
return Err(SymplexError::Unevaluable {
reason: "LambertW of complex argument not yet supported in evalf".into(),
});
}
debug!(prec, "evalf: LambertW via Halley iteration");
let result = arb_lambert_w(&val.0, prec, rm, cc)?;
Ok((result, BigFloat::new(prec)))
}
ExprNode::Beta(a_id, b_id) => {
let a_val = get_cached(cache, *a_id)?;
let b_val = get_cached(cache, *b_id)?;
if !a_val.1.is_zero() || !b_val.1.is_zero() {
return Err(SymplexError::Unevaluable {
reason: "Beta of complex arguments not yet supported in evalf".into(),
});
}
debug!(prec, "evalf: Beta via arbitrary-precision Gamma");
let ga = arb_gamma_real(&a_val.0, prec, rm, cc)?;
let gb = arb_gamma_real(&b_val.0, prec, rm, cc)?;
let a_plus_b = a_val.0.add(&b_val.0, prec, rm);
let gab = arb_gamma_real(&a_plus_b, prec, rm, cc)?;
if gab.is_zero() {
return Err(SymplexError::Unevaluable {
reason: "Beta function: Gamma(a+b) is zero".into(),
});
}
let numer = ga.mul(&gb, prec, rm);
let result = numer.div(&gab, prec, rm);
Ok((result, BigFloat::new(prec)))
}
ExprNode::NaN => Err(SymplexError::Unevaluable {
reason: "cannot evaluate NaN".into(),
}),
ExprNode::Symbol(sid) => {
let name = arena.symbol_name(*sid);
Err(SymplexError::FreeSymbol {
name: name.to_owned(),
})
}
ExprNode::Add(children) => {
let mut sum = c_zero(prec);
for &child in children.iter() {
let val = get_cached(cache, child)?;
sum = c_add(&sum, val, prec, rm);
}
Ok(sum)
}
ExprNode::Mul(children) => {
let mut product = c_one(prec);
for &child in children.iter() {
let val = get_cached(cache, child)?;
product = c_mul(&product, val, prec, rm);
}
Ok(product)
}
ExprNode::Pow(base, exp) => {
let b = get_cached(cache, *base)?;
let e = get_cached(cache, *exp)?;
let b_is_real = b.1.is_zero();
let e_is_real = e.1.is_zero();
if b_is_real
&& !b.0.is_negative()
&& let ExprNode::Num(nid) = arena.node(*exp)
{
let r = arena.num(*nid);
if *r == Ratio::new(BigInt::from(1), BigInt::from(2)) {
return Ok((b.0.sqrt(prec, rm), BigFloat::new(prec)));
}
}
if let Some(n) = try_as_small_int(arena, *exp) {
if b_is_real {
if n >= 0 {
return Ok((b.0.powi(n as usize, prec, rm), BigFloat::new(prec)));
} else {
let pow_pos = b.0.powi((-n) as usize, prec, rm);
let one_bf = BigFloat::from_i32(1, prec);
return Ok((one_bf.div(&pow_pos, prec, rm), BigFloat::new(prec)));
}
} else {
if n >= 0 {
return Ok(c_powi(b, n as usize, prec, rm));
} else {
let pow_pos = c_powi(b, (-n) as usize, prec, rm);
let one_c = c_one(prec);
return Ok(c_div(&one_c, &pow_pos, prec, rm));
}
}
}
if b_is_real && e_is_real && b.0.is_positive() {
return Ok((bf_pow(&b.0, &e.0, prec, rm, cc), BigFloat::new(prec)));
}
Ok(c_pow(b, e, prec, rm, cc))
}
ExprNode::Neg(inner) => {
let val = get_cached(cache, *inner)?;
Ok(c_neg(val))
}
ExprNode::Sin(inner) => {
let val = get_cached(cache, *inner)?;
if val.1.is_zero() {
Ok((val.0.sin(prec, rm, cc), BigFloat::new(prec)))
} else {
Ok(c_sin(val, prec, rm, cc))
}
}
ExprNode::Cos(inner) => {
let val = get_cached(cache, *inner)?;
if val.1.is_zero() {
Ok((val.0.cos(prec, rm, cc), BigFloat::new(prec)))
} else {
Ok(c_cos(val, prec, rm, cc))
}
}
ExprNode::Tan(inner) => {
let val = get_cached(cache, *inner)?;
if val.1.is_zero() {
Ok((val.0.tan(prec, rm, cc), BigFloat::new(prec)))
} else {
let s = c_sin(val, prec, rm, cc);
let c = c_cos(val, prec, rm, cc);
Ok(c_div(&s, &c, prec, rm))
}
}
ExprNode::Exp(inner) => {
let val = get_cached(cache, *inner)?;
if val.1.is_zero() {
Ok((val.0.exp(prec, rm, cc), BigFloat::new(prec)))
} else {
Ok(c_exp(val, prec, rm, cc))
}
}
ExprNode::Ln(inner) => {
let val = get_cached(cache, *inner)?;
if val.1.is_zero() && val.0.is_positive() {
Ok((val.0.ln(prec, rm, cc), BigFloat::new(prec)))
} else {
Ok(c_ln(val, prec, rm, cc))
}
}
ExprNode::Abs(inner) => {
let val = get_cached(cache, *inner)?;
if val.1.is_zero() {
Ok((val.0.abs(), BigFloat::new(prec)))
} else {
Ok((c_abs(val, prec, rm), BigFloat::new(prec)))
}
}
ExprNode::Asin(inner) => {
let val = get_cached(cache, *inner)?;
if val.1.is_zero() {
Ok((val.0.asin(prec, rm, cc), BigFloat::new(prec)))
} else {
Ok(c_asin(val, prec, rm, cc))
}
}
ExprNode::Acos(inner) => {
let val = get_cached(cache, *inner)?;
if val.1.is_zero() {
Ok((val.0.acos(prec, rm, cc), BigFloat::new(prec)))
} else {
Ok(c_acos(val, prec, rm, cc))
}
}
ExprNode::Atan(inner) => {
let val = get_cached(cache, *inner)?;
if val.1.is_zero() {
Ok((val.0.atan(prec, rm, cc), BigFloat::new(prec)))
} else {
Ok(c_atan(val, prec, rm, cc))
}
}
ExprNode::Atan2(y_id, x_id) => {
let yv = get_cached(cache, *y_id)?;
let xv = get_cached(cache, *x_id)?;
if yv.1.is_zero() && xv.1.is_zero() {
Ok((atan2_bf(&yv.0, &xv.0, prec, rm, cc), BigFloat::new(prec)))
} else {
Err(SymplexError::Unevaluable {
reason: "atan2 with complex arguments is not supported".into(),
})
}
}
ExprNode::Sinh(inner) => {
let val = get_cached(cache, *inner)?;
if val.1.is_zero() {
Ok((val.0.sinh(prec, rm, cc), BigFloat::new(prec)))
} else {
Ok(c_sinh(val, prec, rm, cc))
}
}
ExprNode::Cosh(inner) => {
let val = get_cached(cache, *inner)?;
if val.1.is_zero() {
Ok((val.0.cosh(prec, rm, cc), BigFloat::new(prec)))
} else {
Ok(c_cosh(val, prec, rm, cc))
}
}
ExprNode::Tanh(inner) => {
let val = get_cached(cache, *inner)?;
if val.1.is_zero() {
Ok((val.0.tanh(prec, rm, cc), BigFloat::new(prec)))
} else {
let s = c_sinh(val, prec, rm, cc);
let c = c_cosh(val, prec, rm, cc);
Ok(c_div(&s, &c, prec, rm))
}
}
ExprNode::Asinh(inner) => {
let val = get_cached(cache, *inner)?;
if val.1.is_zero() {
Ok((val.0.asinh(prec, rm, cc), BigFloat::new(prec)))
} else {
Ok(c_asinh(val, prec, rm, cc))
}
}
ExprNode::Acosh(inner) => {
let val = get_cached(cache, *inner)?;
if val.1.is_zero() {
Ok((val.0.acosh(prec, rm, cc), BigFloat::new(prec)))
} else {
Ok(c_acosh(val, prec, rm, cc))
}
}
ExprNode::Atanh(inner) => {
let val = get_cached(cache, *inner)?;
if val.1.is_zero() {
Ok((val.0.atanh(prec, rm, cc), BigFloat::new(prec)))
} else {
Ok(c_atanh(val, prec, rm, cc))
}
}
ExprNode::Sign(inner) => {
let val = get_cached(cache, *inner)?;
if val.1.is_zero() {
if val.0.is_zero() {
Ok(c_zero(prec))
} else if val.0.is_negative() {
Ok((BigFloat::from_i32(-1, prec), BigFloat::new(prec)))
} else {
Ok((BigFloat::from_i32(1, prec), BigFloat::new(prec)))
}
} else {
let modulus = c_abs(val, prec, rm);
if modulus.is_zero() {
Ok(c_zero(prec))
} else {
let denom = (modulus, BigFloat::new(prec));
Ok(c_div(val, &denom, prec, rm))
}
}
}
ExprNode::Floor(inner) => {
let val = get_cached(cache, *inner)?;
if !val.1.is_zero() {
return Err(SymplexError::Unevaluable {
reason: "floor of complex number is not defined".into(),
});
}
let truncated = val.0.int();
let result = if val.0.is_negative() && truncated != val.0 {
truncated.sub(&BigFloat::from_i32(1, prec), prec, rm)
} else {
truncated
};
Ok((result, BigFloat::new(prec)))
}
ExprNode::Ceiling(inner) => {
let val = get_cached(cache, *inner)?;
if !val.1.is_zero() {
return Err(SymplexError::Unevaluable {
reason: "ceiling of complex number is not defined".into(),
});
}
let truncated = val.0.int();
let result = if val.0.is_positive() && truncated != val.0 {
truncated.add(&BigFloat::from_i32(1, prec), prec, rm)
} else {
truncated
};
Ok((result, BigFloat::new(prec)))
}
ExprNode::Min(children) => {
if children.is_empty() {
return Err(SymplexError::Unevaluable {
reason: "min of empty set".into(),
});
}
let mut best = get_cached(cache, children[0])?.clone();
if !best.1.is_zero() {
return Err(SymplexError::Unevaluable {
reason: "min requires real arguments".into(),
});
}
for &child in &children[1..] {
let val = get_cached(cache, child)?;
if !val.1.is_zero() {
return Err(SymplexError::Unevaluable {
reason: "min requires real arguments".into(),
});
}
if val.0.sub(&best.0, prec, rm).is_negative() {
best = val.clone();
}
}
Ok(best)
}
ExprNode::Max(children) => {
if children.is_empty() {
return Err(SymplexError::Unevaluable {
reason: "max of empty set".into(),
});
}
let mut best = get_cached(cache, children[0])?.clone();
if !best.1.is_zero() {
return Err(SymplexError::Unevaluable {
reason: "max requires real arguments".into(),
});
}
for &child in &children[1..] {
let val = get_cached(cache, child)?;
if !val.1.is_zero() {
return Err(SymplexError::Unevaluable {
reason: "max requires real arguments".into(),
});
}
if val.0.sub(&best.0, prec, rm).is_positive() {
best = val.clone();
}
}
Ok(best)
}
ExprNode::Heaviside(inner) => {
let val = get_cached(cache, *inner)?;
if !val.1.is_zero() {
return Err(SymplexError::Unevaluable {
reason: "Heaviside of complex argument not supported in evalf".into(),
});
}
if val.0.is_zero() {
Ok((BigFloat::from_f64(0.5, prec), BigFloat::new(prec)))
} else if val.0.is_negative() {
Ok(c_zero(prec))
} else {
Ok((BigFloat::from_i32(1, prec), BigFloat::new(prec)))
}
}
ExprNode::DiracDelta(inner) => {
let _val = get_cached(cache, *inner)?;
Ok(c_zero(prec))
}
ExprNode::Apply(sid, args) => {
let name = arena.symbol_name(*sid);
match arena.lib_fn(*sid) {
Some(f) if f.arity().accepts(args.len()) => {
eval_lib_fn(f, args, cache, prec, rm, cc)
}
Some(_) => Err(unevaluable(format!(
"{name} called with {} argument(s)",
args.len()
))),
None => Err(unevaluable(format!("cannot evaluate function '{name}'"))),
}
}
ExprNode::Derivative(_, _) => Err(SymplexError::Unevaluable {
reason: "cannot evaluate unevaluated derivative".into(),
}),
ExprNode::Integral(_, _) => Err(SymplexError::Unevaluable {
reason: "cannot evaluate unevaluated integral".into(),
}),
ExprNode::DefiniteIntegral(body_id, var_id, lo_id, hi_id) => {
debug!("evalf: DefiniteIntegral — f64 Gauss–Kronrod quadrature");
let var_name = match arena.node(*var_id) {
ExprNode::Symbol(sid) => arena.symbol_name(*sid).to_string(),
_ => {
return Err(SymplexError::Unevaluable {
reason: "integration variable of a definite integral must be a symbol"
.into(),
});
}
};
let a = definite_bound_f64(arena, cache, *lo_id, rm, cc)?;
let b = definite_bound_f64(arena, cache, *hi_id, rm, cc)?;
let func = crate::output::lambdify::compile_raw(arena, *body_id, &[&var_name])
.map_err(|e| SymplexError::Unevaluable {
reason: format!(
"definite integral body cannot be compiled for quadrature: {e}"
),
})?;
let f = |t: f64| func(&[t]);
let opts = crate::calculus::definite::QuadOpts::default();
let crate::calculus::definite::QuadResult { value, error: err } =
crate::calculus::definite::quadrature(&f, a, b, &opts)?;
let tol = opts.abs_tol.max(opts.rel_tol * value.abs());
if err > 1e3 * tol {
return Err(SymplexError::ComputationFailed {
operation: "evalf",
reason: format!(
"quadrature of the definite integral did not converge: estimate {value} \
with error {err:e} (integral may diverge)"
),
});
}
debug!(value, err, "evalf: DefiniteIntegral evaluated");
Ok((BigFloat::from_f64(value, prec), BigFloat::new(prec)))
}
ExprNode::Sum(body_id, var_id, lo_id, hi_id) => {
debug!("evalf: Sum — attempting finite evaluation");
let lo_val = get_cached(cache, *lo_id)?;
let hi_val = get_cached(cache, *hi_id)?;
if !lo_val.1.is_zero() || !hi_val.1.is_zero() {
return Err(SymplexError::Unevaluable {
reason: "Sum bounds must be real".into(),
});
}
let lo_f = bigfloat_to_f64(&lo_val.0, rm, cc)?;
let hi_f = bigfloat_to_f64(&hi_val.0, rm, cc)?;
let lo_i = lo_f.round() as i64;
let hi_i = hi_f.round() as i64;
if (lo_f - lo_i as f64).abs() > 1e-9 || (hi_f - hi_i as f64).abs() > 1e-9 {
return Err(SymplexError::Unevaluable {
reason: "Sum bounds are not integers".into(),
});
}
if hi_i - lo_i > 10_000 {
return Err(SymplexError::Unevaluable {
reason: "Sum range too large for numerical evaluation (> 10000 terms)".into(),
});
}
let mut acc = c_zero(prec);
for i in lo_i..=hi_i {
let sub_val = (BigFloat::from_f64(i as f64, prec), BigFloat::new(prec));
let term =
evalf_subtree_with_sub(arena, *body_id, *var_id, &sub_val, prec, rm, cc)?;
acc = c_add(&acc, &term, prec, rm);
}
debug!(lo = lo_i, hi = hi_i, "evalf: Sum evaluated");
Ok(acc)
}
ExprNode::Product_(body_id, var_id, lo_id, hi_id) => {
debug!("evalf: Product — attempting finite evaluation");
let lo_val = get_cached(cache, *lo_id)?;
let hi_val = get_cached(cache, *hi_id)?;
if !lo_val.1.is_zero() || !hi_val.1.is_zero() {
return Err(SymplexError::Unevaluable {
reason: "Product bounds must be real".into(),
});
}
let lo_f = bigfloat_to_f64(&lo_val.0, rm, cc)?;
let hi_f = bigfloat_to_f64(&hi_val.0, rm, cc)?;
let lo_i = lo_f.round() as i64;
let hi_i = hi_f.round() as i64;
if (lo_f - lo_i as f64).abs() > 1e-9 || (hi_f - hi_i as f64).abs() > 1e-9 {
return Err(SymplexError::Unevaluable {
reason: "Product bounds are not integers".into(),
});
}
if hi_i - lo_i > 10_000 {
return Err(SymplexError::Unevaluable {
reason: "Product range too large for numerical evaluation (> 10000 terms)"
.into(),
});
}
let mut acc = c_one(prec);
for i in lo_i..=hi_i {
let sub_val = (BigFloat::from_f64(i as f64, prec), BigFloat::new(prec));
let term =
evalf_subtree_with_sub(arena, *body_id, *var_id, &sub_val, prec, rm, cc)?;
acc = c_mul(&acc, &term, prec, rm);
}
debug!(lo = lo_i, hi = hi_i, "evalf: Product evaluated");
Ok(acc)
}
ExprNode::Piecewise(pairs) => {
debug!(
branches = pairs.len(),
"evalf: Piecewise — evaluating conditions"
);
for &(value_id, cond_id) in pairs.iter() {
match decide_condition(arena, cond_id, cache, prec, rm) {
Some(true) => {
return eval_node_or_subtree(arena, value_id, cache, prec, rm, cc);
}
Some(false) => continue,
None => {
return Err(SymplexError::Unevaluable {
reason: format!(
"cannot evaluate piecewise: condition `{}` is undecided",
arena.display(cond_id)
),
});
}
}
}
Err(SymplexError::Unevaluable {
reason: "cannot evaluate piecewise: every condition is false".into(),
})
}
ExprNode::BoolTrue
| ExprNode::BoolFalse
| ExprNode::Gt(_, _)
| ExprNode::Ge(_, _)
| ExprNode::Eq_(_, _)
| ExprNode::Ne(_, _)
| ExprNode::And(_)
| ExprNode::Or(_)
| ExprNode::Not(_) => Err(SymplexError::Unevaluable {
reason: "boolean expression".into(),
}),
ExprNode::EmptySet
| ExprNode::UniversalSet
| ExprNode::Interval(_, _, _)
| ExprNode::FiniteSet(_)
| ExprNode::SetUnion(_)
| ExprNode::SetIntersection(_)
| ExprNode::SetComplement(_, _) => Err(SymplexError::Unevaluable {
reason: "set-valued expressions cannot be numerically evaluated".into(),
}),
ExprNode::Limit(_, _, _) => Err(SymplexError::Unevaluable {
reason: "cannot numerically evaluate unevaluated Limit".into(),
}),
ExprNode::Series(_, _, _, _) => Err(SymplexError::Unevaluable {
reason: "cannot numerically evaluate unevaluated Series".into(),
}),
ExprNode::LaplaceTransform(_, _, _) => Err(SymplexError::Unevaluable {
reason: "cannot numerically evaluate unevaluated LaplaceTransform".into(),
}),
ExprNode::InverseLaplaceTransform(_, _, _) => Err(SymplexError::Unevaluable {
reason: "cannot numerically evaluate unevaluated InverseLaplaceTransform".into(),
}),
ExprNode::Residue(_, _, _) => Err(SymplexError::Unevaluable {
reason: "cannot numerically evaluate unevaluated Residue".into(),
}),
ExprNode::RootOf(poly_id, idx_id) => {
let idx: usize = match arena.node(*idx_id) {
ExprNode::Num(nid) => {
let r = arena.num(*nid);
if r.is_integer() {
let n = r.to_integer();
let n_i64: i64 = n.try_into().map_err(|_| SymplexError::Unevaluable {
reason: "RootOf index out of range".into(),
})?;
if n_i64 < 0 {
return Err(SymplexError::Unevaluable {
reason: "RootOf index must be non-negative".into(),
});
}
n_i64 as usize
} else {
return Err(SymplexError::Unevaluable {
reason: "RootOf index must be an integer".into(),
});
}
}
_ => {
return Err(SymplexError::Unevaluable {
reason: "RootOf index must be numeric".into(),
});
}
};
let syms = walk::free_symbols(arena, *poly_id);
if syms.is_empty() {
return Err(SymplexError::Unevaluable {
reason: "RootOf polynomial has no variables".into(),
});
}
let var_id = syms[0];
let poly = match crate::poly::polybridge::expr_to_poly(arena, *poly_id, var_id) {
Some(p) => p,
None => {
return Err(SymplexError::Unevaluable {
reason: "could not convert RootOf expression to polynomial".into(),
});
}
};
let roots = crate::poly::roots::rootof_roots(&poly, prec);
if idx >= roots.len() {
return Err(SymplexError::Unevaluable {
reason: format!(
"RootOf index {} exceeds the {} root(s) found",
idx,
roots.len()
),
});
}
let (re, im) = &roots[idx];
Ok((re.clone(), im.clone()))
}
ExprNode::RootSum(poly, body, sumvar) => {
let poly_id = *poly;
let body_id = *body;
let sumvar_id = *sumvar;
tracing::debug!("evalf: RootSum — attempting numerical evaluation via Aberth roots");
let poly_obj = crate::poly::polybridge::expr_to_poly(arena, poly_id, sumvar_id)
.ok_or_else(|| SymplexError::Unevaluable {
reason: "RootSum: cannot convert polynomial expression to Poly \
(may contain free symbols)"
.into(),
})?;
let poly_deg = poly_obj.degree().unwrap_or(0);
tracing::debug!(
degree = poly_deg,
"evalf: RootSum — polynomial extracted, finding roots"
);
let roots = crate::poly::roots::aberth_roots(&poly_obj, prec, 200);
if roots.len() != poly_deg {
tracing::debug!(
expected = poly_deg,
found = roots.len(),
"evalf: RootSum — Aberth returned fewer roots than expected"
);
}
let mut sum = c_zero(prec);
for (k, root) in roots.iter().enumerate() {
let term = evalf_subtree_with_sub(arena, body_id, sumvar_id, root, prec, rm, cc)?;
tracing::trace!(root_idx = k, "evalf: RootSum — evaluated body at root");
sum = c_add(&sum, &term, prec, rm);
}
tracing::debug!(
n_roots = roots.len(),
"evalf: RootSum — numerical evaluation complete"
);
Ok(sum)
}
ExprNode::DSolve(_, _, _) => Err(SymplexError::Unevaluable {
reason: "cannot numerically evaluate unevaluated DSolve".into(),
}),
ExprNode::ConditionSet(_, _) => Err(SymplexError::Unevaluable {
reason: "cannot numerically evaluate ConditionSet".into(),
}),
}
}
fn contains_definite_integral(arena: &Arena, root: ExprId) -> bool {
let mut stack = vec![root];
let mut seen = rustc_hash::FxHashSet::default();
while let Some(id) = stack.pop() {
if !seen.insert(id) {
continue;
}
let node = arena.node(id);
if matches!(node, ExprNode::DefiniteIntegral(..)) {
return true;
}
node.for_each_child(|c| stack.push(c));
}
false
}
fn definite_bound_f64(
arena: &Arena,
cache: &FxHashMap<ExprId, Complex>,
id: ExprId,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<f64, SymplexError> {
match arena.node(id) {
ExprNode::Infinity => Ok(f64::INFINITY),
ExprNode::NegInfinity => Ok(f64::NEG_INFINITY),
_ => {
let v = get_cached(cache, id)?;
if !v.1.is_zero() {
return Err(SymplexError::Unevaluable {
reason: "integration bounds of a definite integral must be real".into(),
});
}
bigfloat_to_f64(&v.0, rm, cc)
}
}
}
fn evalf_subtree_with_sub(
arena: &Arena,
expr: ExprId,
var_id: ExprId,
var_value: &Complex,
prec: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<Complex, SymplexError> {
let post_order = walk::post_order_ids(arena, expr);
let mut sub_cache: FxHashMap<ExprId, Complex> = FxHashMap::default();
sub_cache.insert(var_id, var_value.clone());
for &id in &post_order {
if sub_cache.contains_key(&id) {
continue; }
let value = eval_node(arena, id, &sub_cache, prec, rm, cc)?;
sub_cache.insert(id, value);
}
sub_cache
.get(&expr)
.cloned()
.ok_or_else(|| SymplexError::Unevaluable {
reason: "subtree evaluation with substitution failed".into(),
})
}
fn decide_condition(
arena: &Arena,
cond: ExprId,
cache: &FxHashMap<ExprId, Complex>,
prec: usize,
rm: RoundingMode,
) -> Option<bool> {
let real_diff = |a: &ExprId, b: &ExprId| -> Option<BigFloat> {
let (av, bv) = (cache.get(a)?, cache.get(b)?);
if av.1.is_zero() && bv.1.is_zero() {
Some(av.0.sub(&bv.0, prec, rm))
} else {
None
}
};
let complex_equal = |a: &ExprId, b: &ExprId| -> Option<bool> {
let (av, bv) = (cache.get(a)?, cache.get(b)?);
let d_re = av.0.sub(&bv.0, prec, rm);
let d_im = av.1.sub(&bv.1, prec, rm);
Some(d_re.is_zero() && d_im.is_zero())
};
let order = walk::post_order_ids(arena, cond);
let mut truth: FxHashMap<ExprId, Option<bool>> = FxHashMap::default();
for &id in &order {
let v: Option<bool> = match arena.node(id) {
ExprNode::BoolTrue => Some(true),
ExprNode::BoolFalse => Some(false),
ExprNode::Gt(a, b) => real_diff(a, b).map(|d| d.is_positive()),
ExprNode::Ge(a, b) => real_diff(a, b).map(|d| d.is_positive() || d.is_zero()),
ExprNode::Eq_(a, b) => complex_equal(a, b),
ExprNode::Ne(a, b) => complex_equal(a, b).map(|e| !e),
ExprNode::Not(inner) => truth.get(inner).copied().flatten().map(|b| !b),
ExprNode::And(kids) => {
let vals: Vec<Option<bool>> = kids
.iter()
.map(|k| truth.get(k).copied().flatten())
.collect();
if vals.contains(&Some(false)) {
Some(false)
} else if vals.iter().all(|v| *v == Some(true)) {
Some(true)
} else {
None
}
}
ExprNode::Or(kids) => {
let vals: Vec<Option<bool>> = kids
.iter()
.map(|k| truth.get(k).copied().flatten())
.collect();
if vals.contains(&Some(true)) {
Some(true)
} else if vals.iter().all(|v| *v == Some(false)) {
Some(false)
} else {
None
}
}
_ => None,
};
truth.insert(id, v);
}
truth.get(&cond).copied().flatten()
}
fn eval_node_or_subtree(
arena: &Arena,
id: ExprId,
cache: &FxHashMap<ExprId, Complex>,
prec: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<Complex, SymplexError> {
if let Some(val) = cache.get(&id) {
return Ok(val.clone());
}
eval_node(arena, id, cache, prec, rm, cc)
}
pub(crate) use crate::base::bigcomplex::{
c_abs, c_add, c_div, c_from_real, c_i, c_mul, c_neg, c_one, c_powi, c_sub, c_zero,
};
fn c_exp(z: &Complex, prec: usize, rm: RoundingMode, cc: &mut Consts) -> Complex {
let exp_a = z.0.exp(prec, rm, cc);
let cos_b = z.1.cos(prec, rm, cc);
let sin_b = z.1.sin(prec, rm, cc);
(exp_a.mul(&cos_b, prec, rm), exp_a.mul(&sin_b, prec, rm))
}
fn c_ln(z: &Complex, prec: usize, rm: RoundingMode, cc: &mut Consts) -> Complex {
let r = c_abs(z, prec, rm);
let theta = atan2_bf(&z.1, &z.0, prec, rm, cc);
(r.ln(prec, rm, cc), theta)
}
fn c_sin(z: &Complex, prec: usize, rm: RoundingMode, cc: &mut Consts) -> Complex {
let sin_a = z.0.sin(prec, rm, cc);
let cos_a = z.0.cos(prec, rm, cc);
let cosh_b = z.1.cosh(prec, rm, cc);
let sinh_b = z.1.sinh(prec, rm, cc);
(sin_a.mul(&cosh_b, prec, rm), cos_a.mul(&sinh_b, prec, rm))
}
fn c_cos(z: &Complex, prec: usize, rm: RoundingMode, cc: &mut Consts) -> Complex {
let cos_a = z.0.cos(prec, rm, cc);
let sin_a = z.0.sin(prec, rm, cc);
let cosh_b = z.1.cosh(prec, rm, cc);
let sinh_b = z.1.sinh(prec, rm, cc);
(
cos_a.mul(&cosh_b, prec, rm),
sin_a.mul(&sinh_b, prec, rm).neg(),
)
}
fn c_sinh(z: &Complex, prec: usize, rm: RoundingMode, cc: &mut Consts) -> Complex {
let sinh_a = z.0.sinh(prec, rm, cc);
let cosh_a = z.0.cosh(prec, rm, cc);
let cos_b = z.1.cos(prec, rm, cc);
let sin_b = z.1.sin(prec, rm, cc);
(sinh_a.mul(&cos_b, prec, rm), cosh_a.mul(&sin_b, prec, rm))
}
fn c_cosh(z: &Complex, prec: usize, rm: RoundingMode, cc: &mut Consts) -> Complex {
let cosh_a = z.0.cosh(prec, rm, cc);
let sinh_a = z.0.sinh(prec, rm, cc);
let cos_b = z.1.cos(prec, rm, cc);
let sin_b = z.1.sin(prec, rm, cc);
(cosh_a.mul(&cos_b, prec, rm), sinh_a.mul(&sin_b, prec, rm))
}
fn c_pow(base: &Complex, exp: &Complex, prec: usize, rm: RoundingMode, cc: &mut Consts) -> Complex {
if base.0.is_zero() && base.1.is_zero() {
if exp.0.is_zero() && exp.1.is_zero() {
return c_one(prec);
}
return c_zero(prec);
}
let ln_base = c_ln(base, prec, rm, cc);
let product = c_mul(exp, &ln_base, prec, rm);
c_exp(&product, prec, rm, cc)
}
fn c_sqrt(z: &Complex, prec: usize, rm: RoundingMode, cc: &mut Consts) -> Complex {
if z.0.is_zero() && z.1.is_zero() {
return c_zero(prec);
}
let half = c_from_real(
BigFloat::from_i32(1, prec).div(&BigFloat::from_i32(2, prec), prec, rm),
prec,
);
c_pow(z, &half, prec, rm, cc)
}
fn c_asin(z: &Complex, prec: usize, rm: RoundingMode, cc: &mut Consts) -> Complex {
let i_unit = c_i(prec);
let one = c_one(prec);
let z_sq = c_mul(z, z, prec, rm);
let one_minus_z_sq = c_sub(&one, &z_sq, prec, rm);
let sqrt_term = c_sqrt(&one_minus_z_sq, prec, rm, cc);
let iz = c_mul(&i_unit, z, prec, rm);
let sum = c_add(&iz, &sqrt_term, prec, rm);
let ln_sum = c_ln(&sum, prec, rm, cc);
let neg_i = c_neg(&i_unit);
c_mul(&neg_i, &ln_sum, prec, rm)
}
fn c_acos(z: &Complex, prec: usize, rm: RoundingMode, cc: &mut Consts) -> Complex {
let half_pi_val = cc.pi(prec, rm).div(&BigFloat::from_i32(2, prec), prec, rm);
let half_pi = c_from_real(half_pi_val, prec);
let asin_z = c_asin(z, prec, rm, cc);
c_sub(&half_pi, &asin_z, prec, rm)
}
fn c_atan(z: &Complex, prec: usize, rm: RoundingMode, cc: &mut Consts) -> Complex {
let i_unit = c_i(prec);
let two = c_from_real(BigFloat::from_i32(2, prec), prec);
let i_over_2 = c_div(&i_unit, &two, prec, rm);
let i_plus_z = c_add(&i_unit, z, prec, rm);
let i_minus_z = c_sub(&i_unit, z, prec, rm);
let ratio = c_div(&i_plus_z, &i_minus_z, prec, rm);
let ln_ratio = c_ln(&ratio, prec, rm, cc);
c_mul(&i_over_2, &ln_ratio, prec, rm)
}
fn c_asinh(z: &Complex, prec: usize, rm: RoundingMode, cc: &mut Consts) -> Complex {
let one = c_one(prec);
let z_sq = c_mul(z, z, prec, rm);
let z_sq_plus_one = c_add(&z_sq, &one, prec, rm);
let sqrt_term = c_sqrt(&z_sq_plus_one, prec, rm, cc);
let sum = c_add(z, &sqrt_term, prec, rm);
c_ln(&sum, prec, rm, cc)
}
fn c_acosh(z: &Complex, prec: usize, rm: RoundingMode, cc: &mut Consts) -> Complex {
let one = c_one(prec);
let z_sq = c_mul(z, z, prec, rm);
let z_sq_minus_one = c_sub(&z_sq, &one, prec, rm);
let sqrt_term = c_sqrt(&z_sq_minus_one, prec, rm, cc);
let sum = c_add(z, &sqrt_term, prec, rm);
c_ln(&sum, prec, rm, cc)
}
fn c_atanh(z: &Complex, prec: usize, rm: RoundingMode, cc: &mut Consts) -> Complex {
let one = c_one(prec);
let two = c_from_real(BigFloat::from_i32(2, prec), prec);
let half = c_div(&one, &two, prec, rm);
let one_plus_z = c_add(&one, z, prec, rm);
let one_minus_z = c_sub(&one, z, prec, rm);
let ratio = c_div(&one_plus_z, &one_minus_z, prec, rm);
let ln_ratio = c_ln(&ratio, prec, rm, cc);
c_mul(&half, &ln_ratio, prec, rm)
}
fn atan2_bf(
y: &BigFloat,
x: &BigFloat,
prec: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> BigFloat {
if x.is_zero() && y.is_zero() {
return BigFloat::new(prec);
}
let pi_val = cc.pi(prec, rm);
let two = BigFloat::from_i32(2, prec);
let half_pi = pi_val.div(&two, prec, rm);
if x.is_zero() {
return if y.is_positive() {
half_pi
} else {
half_pi.neg()
};
}
let ratio = y.div(x, prec, rm);
let atan_val = ratio.atan(prec, rm, cc);
if x.is_positive() {
atan_val
} else {
let pi_val2 = cc.pi(prec, rm);
if y.is_negative() {
atan_val.sub(&pi_val2, prec, rm)
} else {
atan_val.add(&pi_val2, prec, rm)
}
}
}
fn bigfloat_to_f64(bf: &BigFloat, rm: RoundingMode, _cc: &mut Consts) -> Result<f64, SymplexError> {
bigfloat_to_f64_rounded(bf, rm)
}
fn bigfloat_to_f64_rounded(bf: &BigFloat, rm: RoundingMode) -> Result<f64, SymplexError> {
if bf.is_zero() {
return Ok(0.0);
}
if bf.is_inf_pos() {
return Ok(f64::INFINITY);
}
if bf.is_inf_neg() {
return Ok(f64::NEG_INFINITY);
}
let Some((words, _, sign, exponent, _)) = bf.as_raw_parts() else {
return Err(SymplexError::Unevaluable {
reason: "NaN in BigFloat to f64 conversion".into(),
});
};
let negative = sign == Sign::Neg;
let signed = |v: f64| if negative { -v } else { v };
let mut acc: u128 = 0;
let mut filled = 0usize;
let mut sticky = false;
for &w in words.iter().rev() {
if filled < 64 {
acc = (acc << WORD_BIT_SIZE) | w as u128;
filled += WORD_BIT_SIZE;
} else if w != 0 {
sticky = true;
}
}
let top: u64 = if filled >= 64 {
let extra = filled - 64;
if extra > 0 && acc & ((1u128 << extra) - 1) != 0 {
sticky = true;
}
(acc >> extra) as u64
} else {
(acc << (64 - filled)) as u64
};
if top & (1 << 63) == 0 {
return Ok(signed(0.0));
}
let mut unbiased = exponent as i64 - 1;
if unbiased > 1023 {
return Ok(signed(f64::INFINITY));
}
let keep = if unbiased >= -1022 {
53
} else {
unbiased + 1075
};
let (mant, above_half, exactly_half) = if keep <= 0 {
let at_half = top == 1 << 63 && !sticky;
(0u64, keep == 0 && !at_half, keep == 0 && at_half)
} else {
let shift = (64 - keep) as u32;
let mant = top >> shift;
let rem = top & ((1u64 << shift) - 1);
let half = 1u64 << (shift - 1);
(
mant,
rem > half || (rem == half && sticky),
rem == half && !sticky,
)
};
let round_up = match rm {
RoundingMode::None => false,
_ if above_half => true,
_ if !exactly_half => false,
RoundingMode::Up => !negative,
RoundingMode::Down => negative,
RoundingMode::ToZero => false,
RoundingMode::FromZero => true,
RoundingMode::ToEven => mant & 1 == 1,
RoundingMode::ToOdd => mant & 1 == 0,
};
let mut mant = mant + u64::from(round_up);
let bits = if keep == 53 {
if mant == 1 << 53 {
mant >>= 1;
unbiased += 1;
if unbiased > 1023 {
return Ok(signed(f64::INFINITY));
}
}
(((unbiased + 1023) as u64) << 52) | (mant & ((1 << 52) - 1))
} else {
mant
};
Ok(signed(f64::from_bits(bits)))
}
#[allow(dead_code)]
#[allow(clippy::excessive_precision)]
fn lanczos_gamma_f64(x: f64) -> Result<f64, SymplexError> {
if x.is_nan() || x.is_infinite() {
return Err(SymplexError::Unevaluable {
reason: "Gamma of special float value".into(),
});
}
if x < 0.5 {
let reflected = lanczos_gamma_f64(1.0 - x)?;
let sin_pi_x = (std::f64::consts::PI * x).sin();
if sin_pi_x.abs() < 1e-300 {
return Err(SymplexError::Unevaluable {
reason: "Gamma at non-positive integer pole".into(),
});
}
let result = std::f64::consts::PI / (sin_pi_x * reflected);
if result.is_finite() {
return Ok(result);
}
return Err(SymplexError::Unevaluable {
reason: "Gamma at pole".into(),
});
}
const COEFFICIENTS: [f64; 9] = [
0.99999999999980993,
676.5203681218851,
-1259.1392167224028,
771.32342877765313,
-176.61502916214059,
12.507343278686905,
-0.13857109526572012,
9.9843695780195716e-6,
1.5056327351493116e-7,
];
let x = x - 1.0; let t = x + 7.0 + 0.5;
let mut sum = COEFFICIENTS[0];
for (i, &c) in COEFFICIENTS[1..].iter().enumerate() {
sum += c / (x + i as f64 + 1.0);
}
let result = (2.0 * std::f64::consts::PI).sqrt() * t.powf(x + 0.5) * (-t).exp() * sum;
Ok(result)
}
#[allow(dead_code)]
fn erf_f64(x: f64) -> f64 {
if x.abs() < 4.0 {
let mut sum = 0.0;
let mut term = x; sum += term;
for n in 1..50 {
term *= -x * x / n as f64;
sum += term / (2 * n + 1) as f64;
}
sum * 2.0 / std::f64::consts::PI.sqrt()
} else {
let sign = x.signum();
let ax = x.abs();
let mut sum = 1.0;
let mut term = 1.0;
for n in 1..20 {
term *= -(2 * n - 1) as f64 / (2.0 * ax * ax);
if term.abs() < 1e-16 {
break;
}
sum += term;
}
let erfc = (-ax * ax).exp() / (ax * std::f64::consts::PI.sqrt()) * sum;
sign * (1.0 - erfc)
}
}
fn arb_digamma(
x: &BigFloat,
prec: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
let guard = 20;
let wp = prec + guard;
if x.is_nan() || x.is_inf_pos() || x.is_inf_neg() {
return Err(SymplexError::Unevaluable {
reason: "Digamma of special float value".into(),
});
}
let one = BigFloat::from_i32(1, wp);
if x.is_negative() {
let pi_val = cc.pi(wp, rm).clone();
let pi_x = pi_val.mul(x, wp, rm);
let sin_val = pi_x.sin(wp, rm, cc);
if sin_val.is_zero() {
return Err(SymplexError::Unevaluable {
reason: "Digamma at non-positive integer pole".into(),
});
}
if let Some(s_exp) = sin_val.exponent()
&& (s_exp as i64) < -(wp as i64 / 2)
{
return Err(SymplexError::Unevaluable {
reason: "Digamma at non-positive integer pole".into(),
});
}
let cos_val = pi_x.cos(wp, rm, cc);
let cot_val = cos_val.div(&sin_val, wp, rm);
let one_minus_x = one.sub(x, wp, rm);
let psi_1mx = arb_digamma(&one_minus_x, prec, rm, cc)?;
let pi_cot = cc.pi(wp, rm).clone().mul(&cot_val, wp, rm);
return Ok(psi_1mx.sub(&pi_cot, wp, rm));
}
let threshold_val = (wp as i32 / 3).max(10);
let threshold = BigFloat::from_i32(threshold_val, wp);
let mut result = BigFloat::new(wp); let mut x = x.clone();
while x.sub(&threshold, wp, rm).is_negative() {
if x.is_zero() {
return Err(SymplexError::Unevaluable {
reason: "Digamma at non-positive integer pole".into(),
});
}
if let Some(x_exp) = x.exponent()
&& (x_exp as i64) < -(wp as i64 / 2)
{
return Err(SymplexError::Unevaluable {
reason: "Digamma at non-positive integer pole".into(),
});
}
let inv_x = one.div(&x, wp, rm);
result = result.sub(&inv_x, wp, rm);
x = x.add(&one, wp, rm);
}
let ln_x = x.ln(wp, rm, cc);
result = result.add(&ln_x, wp, rm);
let two = BigFloat::from_i32(2, wp);
let half_inv_x = one.div(&x.mul(&two, wp, rm), wp, rm);
result = result.sub(&half_inv_x, wp, rm);
let x2 = x.mul(&x, wp, rm);
let mut x_pow = x2.clone();
let bernoulli_nums: &[(i128, i128)] = &[
(1, 6),
(-1, 30),
(1, 42),
(-1, 30),
(5, 66),
(-691, 2730),
(7, 6),
(-3617, 510),
(43867, 798),
(-174611, 330),
(854513, 138),
(-236364091, 2730),
];
let n_terms = (prec / 6 + 2).min(bernoulli_nums.len());
for (k_idx, &(bn, bd)) in bernoulli_nums[..n_terms].iter().enumerate() {
let k = (k_idx + 1) as i128;
let two_k = 2 * k;
let coeff_n = BigFloat::from_i128(bn, wp);
let coeff_d = BigFloat::from_i128(bd * two_k, wp);
let coeff = coeff_n.div(&coeff_d, wp, rm);
let inv_xpow = one.div(&x_pow, wp, rm);
let term = coeff.mul(&inv_xpow, wp, rm);
result = result.sub(&term, wp, rm);
if let (Some(t_exp), Some(r_exp)) = (term.exponent(), result.exponent())
&& (r_exp as i64 - t_exp as i64) > wp as i64
{
break;
}
x_pow = x_pow.mul(&x2, wp, rm);
}
Ok(result)
}
fn stirling_log_gamma(
z: &BigFloat,
prec: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
tracing::debug!(prec, "stirling_log_gamma: entry");
let guard = 20; let wp = prec + guard;
let z_approx = bigfloat_to_f64(z, rm, cc)?;
let threshold = ((wp as f64) * 0.2).ceil() as i64;
let shift = if z_approx < threshold as f64 {
(threshold - z_approx.floor() as i64).max(0) as usize
} else {
0
};
tracing::trace!(
z_approx,
threshold,
shift,
"stirling_log_gamma: argument reduction"
);
let shift_bf = BigFloat::from_i128(shift as i128, wp);
let z_shifted = z.add(&shift_bf, wp, rm);
let half = BigFloat::from_f64(0.5, wp);
let z_minus_half = z_shifted.sub(&half, wp, rm);
let log_z = z_shifted.ln(wp, rm, cc);
let mut result = z_minus_half.mul(&log_z, wp, rm);
result = result.sub(&z_shifted, wp, rm);
let two = BigFloat::from_i32(2, wp);
let pi_val = cc.pi(wp, rm).clone();
let two_pi = two.mul(&pi_val, wp, rm);
let half_log_2pi = two_pi.ln(wp, rm, cc).mul(&half, wp, rm);
result = result.add(&half_log_2pi, wp, rm);
let z_sq = z_shifted.mul(&z_shifted, wp, rm);
let mut z_pow = z_shifted.clone();
let max_terms = ((wp as f64) * 0.35) as usize + 10;
let mut prev_term_exp: Option<i32> = None;
for k in 1..=max_terms {
let b2k = crate::base::bernoulli::bernoulli(2 * k);
if b2k.is_zero() {
continue;
}
if k > 1 {
z_pow = z_pow.mul(&z_sq, wp, rm);
}
let b2k_bf = ratio_to_bigfloat(&b2k, wp, rm);
let denom_int = (2 * k * (2 * k - 1)) as i128;
let denom = BigFloat::from_i128(denom_int, wp);
let denom_full = denom.mul(&z_pow, wp, rm);
let term = b2k_bf.div(&denom_full, wp, rm);
if let Some(t_exp) = term.exponent() {
if let Some(prev) = prev_term_exp
&& t_exp > prev + 10
{
tracing::trace!(k, "stirling series: diverging, stopping");
break;
}
prev_term_exp = Some(t_exp);
if let Some(r_exp) = result.exponent()
&& (r_exp as i64 - t_exp as i64) > wp as i64
{
tracing::trace!(k, "stirling series: converged");
break;
}
}
result = result.add(&term, wp, rm);
}
for i in 0..shift {
let z_plus_i = z.add(&BigFloat::from_i128(i as i128, wp), wp, rm);
let log_zi = z_plus_i.ln(wp, rm, cc);
result = result.sub(&log_zi, wp, rm);
}
Ok(result)
}
fn arb_gamma_real(
x: &BigFloat,
prec: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
if x.is_zero() {
return Err(SymplexError::Unevaluable {
reason: "Gamma at non-positive integer pole".into(),
});
}
let x_approx = bigfloat_to_f64(x, rm, cc)?;
if x_approx <= 0.0 {
let rounded = x_approx.round();
if (rounded - x_approx).abs() < 1e-12 && rounded <= 0.0 {
return Err(SymplexError::Unevaluable {
reason: "Gamma at non-positive integer pole".into(),
});
}
}
if x_approx > 0.0 {
tracing::debug!(x_approx, "arb_gamma_real: positive argument");
let log_gamma = stirling_log_gamma(x, prec, rm, cc)?;
Ok(log_gamma.exp(prec, rm, cc))
} else {
tracing::debug!(x_approx, "arb_gamma_real: reflection formula");
let one = BigFloat::from_i32(1, prec);
let one_minus_x = one.sub(x, prec, rm);
let log_gamma_1mx = stirling_log_gamma(&one_minus_x, prec, rm, cc)?;
let gamma_1mx = log_gamma_1mx.exp(prec, rm, cc);
let pi_val = cc.pi(prec, rm).clone();
let pi_x = pi_val.mul(x, prec, rm);
let sin_pi_x = pi_x.sin(prec, rm, cc);
if sin_pi_x.is_zero() {
return Err(SymplexError::Unevaluable {
reason: "Gamma at non-positive integer pole".into(),
});
}
if let Some(s_exp) = sin_pi_x.exponent()
&& (s_exp as i64) < -(prec as i64 / 2)
{
return Err(SymplexError::Unevaluable {
reason: "Gamma at non-positive integer pole".into(),
});
}
let denom = sin_pi_x.mul(&gamma_1mx, prec, rm);
let pi_val2 = cc.pi(prec, rm).clone();
Ok(pi_val2.div(&denom, prec, rm))
}
}
fn erfc_asymptotic(ax: &BigFloat, wp: usize, rm: RoundingMode, cc: &mut Consts) -> BigFloat {
let x_sq = ax.mul(ax, wp, rm);
let two_x_sq = x_sq.mul(&BigFloat::from_i32(2, wp), wp, rm);
let mut sum = BigFloat::from_i32(1, wp);
let mut term = BigFloat::from_i32(1, wp);
let mut prev_abs: Option<BigFloat> = None;
let max_terms = wp + 40;
for n in 1..=max_terms {
let factor = BigFloat::from_i128(2 * n as i128 - 1, wp);
term = term.mul(&factor, wp, rm).div(&two_x_sq, wp, rm).neg();
let t_abs = term.abs();
if let Some(ref p) = prev_abs
&& bf_gt(&t_abs, p)
{
break; }
sum = sum.add(&term, wp, rm);
if negligible(&term, &sum, wp) {
break;
}
prev_abs = Some(t_abs);
}
let exp_neg_x_sq = x_sq.neg().exp(wp, rm, cc);
let sqrt_pi = cc.pi(wp, rm).clone().sqrt(wp, rm);
let x_sqrt_pi = ax.mul(&sqrt_pi, wp, rm);
exp_neg_x_sq.mul(&sum, wp, rm).div(&x_sqrt_pi, wp, rm)
}
fn arb_erf(
x: &BigFloat,
prec: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
if x.is_zero() {
return Ok(BigFloat::new(prec));
}
let guard = 32;
let wp = prec + guard;
let x_approx = bigfloat_to_f64(x, rm, cc)?;
let threshold = ((wp as f64) * 0.35).sqrt() + 2.0;
if x_approx.abs() < threshold {
let x_sq_f = x_approx * x_approx;
let cancel = (x_sq_f * std::f64::consts::LOG2_E).ceil() as usize + 8;
let wpt = wp + cancel;
let mut xw = x.clone();
let _ = xw.set_precision(wpt, rm);
let neg_x_sq = xw.mul(&xw, wpt, rm).neg();
let mut prod = xw.clone(); let mut sum = xw.clone();
let max_terms = (6.0 * x_sq_f) as usize + wpt + 60;
let mut converged = false;
for n in 1..=max_terms {
prod = prod.mul(&neg_x_sq, wpt, rm);
prod = prod.div(&BigFloat::from_i128(n as i128, wpt), wpt, rm);
let divisor = BigFloat::from_i128(2 * n as i128 + 1, wpt);
let term = prod.div(&divisor, wpt, rm);
if negligible(&term, &sum, wpt) {
converged = true;
break;
}
sum = sum.add(&term, wpt, rm);
}
if !converged {
return Err(series_did_not_converge("erf", max_terms));
}
let two = BigFloat::from_i32(2, wpt);
let pi_val = cc.pi(wpt, rm).clone();
let sqrt_pi = pi_val.sqrt(wpt, rm);
let two_over_sqrt_pi = two.div(&sqrt_pi, wpt, rm);
Ok(round_to(sum.mul(&two_over_sqrt_pi, wpt, rm), wp, rm))
} else {
let mut ax = x.abs();
let _ = ax.set_precision(wp, rm);
let erfc_val = erfc_asymptotic(&ax, wp, rm, cc);
let one = BigFloat::from_i32(1, wp);
let erf_val = one.sub(&erfc_val, wp, rm);
if x.is_negative() {
Ok(erf_val.neg())
} else {
Ok(erf_val)
}
}
}
fn arb_lambert_w(
x: &BigFloat,
prec: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
let wp = prec + 32;
if x.is_zero() {
return Ok(BigFloat::new(prec));
}
let neg_inv_e = {
let e_val = cc.e(wp, rm).clone();
let one = BigFloat::from_i32(1, wp);
one.div(&e_val, wp, rm).neg()
};
let diff = x.sub(&neg_inv_e, wp, rm);
if diff.is_negative() {
return Err(SymplexError::Unevaluable {
reason: "LambertW: argument < -1/e, outside principal branch domain".into(),
});
}
let mut w = if x.is_negative() {
let e_val = cc.e(wp, rm).clone();
let ex = e_val.mul(x, wp, rm);
let one = BigFloat::from_i32(1, wp);
let two = BigFloat::from_i32(2, wp);
let inner = one.add(&ex, wp, rm).mul(&two, wp, rm);
if inner.is_negative() || inner.is_zero() {
BigFloat::from_i32(-1, wp)
} else {
let sq = inner.sqrt(wp, rm);
BigFloat::from_i32(-1, wp).add(&sq, wp, rm)
}
} else {
let threshold = BigFloat::from_f64(2.5, wp);
if x.sub(&threshold, wp, rm).is_negative() {
x.clone()
} else {
let ln_x = x.ln(wp, rm, cc);
let ln_ln_x = ln_x.ln(wp, rm, cc);
ln_x.sub(&ln_ln_x, wp, rm)
}
};
let max_iter = 100;
for _ in 0..max_iter {
let ew = w.exp(wp, rm, cc);
let wew = w.mul(&ew, wp, rm);
let residual = wew.sub(x, wp, rm);
let w_plus_1 = w.add(&BigFloat::from_i32(1, wp), wp, rm);
let denom_base = w_plus_1.mul(&ew, wp, rm);
let w_plus_2 = w.add(&BigFloat::from_i32(2, wp), wp, rm);
let two_w_plus_1 = w_plus_1.mul(&BigFloat::from_i32(2, wp), wp, rm);
let correction_numer = w_plus_2.mul(&residual, wp, rm);
let correction = if two_w_plus_1.is_zero() {
BigFloat::new(wp)
} else {
correction_numer.div(&two_w_plus_1, wp, rm)
};
let denom = denom_base.sub(&correction, wp, rm);
if denom.is_zero() {
break;
}
let delta = residual.div(&denom, wp, rm);
w = w.sub(&delta, wp, rm);
if delta.is_zero() {
break;
}
if let (Some(d_exp), Some(w_e)) = (delta.exponent(), w.exponent())
&& (d_exp as i64) < (w_e as i64) - (wp as i64)
{
break;
}
}
Ok(w)
}
fn arb_euler_gamma(
prec: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
let extra = 30;
let wp = prec + extra;
let p = ((wp as f64 / 4.0) * std::f64::consts::LN_2).log2().ceil() as u64 + 1;
let n: i128 = 1i128 << p;
let n_sq = n * n;
let two_bf = BigFloat::from_i32(2, wp);
let ln2 = two_bf.ln(wp, rm, cc);
let p_bf = BigFloat::from_i128(p as i128, wp);
let neg_p_ln2 = p_bf.mul(&ln2, wp, rm).neg();
let mut a = neg_p_ln2.clone();
let mut u = neg_p_ln2;
let one = BigFloat::from_i32(1, wp);
let mut b = one.clone();
let mut v = one.clone();
let n_sq_bf = BigFloat::from_i128(n_sq, wp);
let mut k: i128 = 1;
loop {
let k_bf = BigFloat::from_i128(k, wp);
let k_sq_bf = k_bf.mul(&k_bf, wp, rm);
b = b.mul(&n_sq_bf, wp, rm).div(&k_sq_bf, wp, rm);
a = a.mul(&n_sq_bf, wp, rm).div(&k_bf, wp, rm);
a = a.add(&b, wp, rm).div(&k_bf, wp, rm);
u = u.add(&a, wp, rm);
v = v.add(&b, wp, rm);
let a_small = a.exponent().is_none_or(|e| e < -(wp as i32) + 10);
let b_small = b.exponent().is_none_or(|e| e < -(wp as i32) + 10);
if a_small && b_small {
tracing::trace!(k, "arb_euler_gamma: converged");
break;
}
k += 1;
if k > 10 * (wp as i128) {
return Err(SymplexError::ComputationFailed {
operation: "euler_gamma",
reason: "Brent-McMillan did not converge".into(),
});
}
}
Ok(u.div(&v, prec, rm))
}
fn hankel_pq(
order: &BigFloat,
z: &BigFloat,
prec: usize,
rm: RoundingMode,
_cc: &mut Consts,
) -> (BigFloat, BigFloat, bool) {
let one = BigFloat::from_i32(1, prec);
let four = BigFloat::from_i32(4, prec);
let eight = BigFloat::from_i32(8, prec);
let four_nu_sq = four.mul(&order.mul(order, prec, rm), prec, rm);
let z_inv = one.div(z, prec, rm);
let mut p_sum = BigFloat::new(prec);
let mut q_sum = BigFloat::new(prec);
let mut a_k = one.clone(); let mut z_power = one.clone();
let max_terms = (3 * prec + 100).min(100_000);
let mut prev_a_abs: Option<BigFloat> = None;
let mut converged = false;
for k in 0..max_terms {
if k > 0 {
let two_km1 = BigFloat::from_i32((2 * k as i32) - 1, prec);
let two_km1_sq = two_km1.mul(&two_km1, prec, rm);
let numer = four_nu_sq.sub(&two_km1_sq, prec, rm);
let k_bf = BigFloat::from_i32(k as i32, prec);
let denom = eight.mul(&k_bf, prec, rm);
a_k = a_k.mul(&numer, prec, rm).div(&denom, prec, rm);
z_power = z_power.mul(&z_inv, prec, rm);
}
let term = a_k.mul(&z_power, prec, rm);
let term_abs = term.abs();
if let Some(ref prev) = prev_a_abs
&& term_abs.partial_cmp(prev) == Some(std::cmp::Ordering::Greater)
&& k > 2
{
tracing::trace!(k, "hankel_pq: optimal truncation (terms diverging)");
break;
}
let p_exp = p_sum.exponent().unwrap_or(0);
let q_exp = q_sum.exponent().unwrap_or(0);
let ref_exp = p_exp.max(q_exp);
let t_exp = term.exponent().unwrap_or(i32::MIN / 2);
if k > 0 && (ref_exp as i64 - t_exp as i64) > prec as i64 {
tracing::trace!(k, "hankel_pq: converged (term negligible)");
converged = true;
break;
}
prev_a_abs = Some(term_abs);
if k % 2 == 0 {
if (k / 2) % 2 == 0 {
p_sum = p_sum.add(&term, prec, rm);
} else {
p_sum = p_sum.sub(&term, prec, rm);
}
} else {
if ((k - 1) / 2) % 2 == 0 {
q_sum = q_sum.add(&term, prec, rm);
} else {
q_sum = q_sum.sub(&term, prec, rm);
}
}
}
(p_sum, q_sum, converged)
}
fn hankel_threshold(wp: usize, order_abs: f64) -> f64 {
(wp as f64 + 3.0 * order_abs) * std::f64::consts::LN_2 / 2.0 + 4.0
}
fn bessel_jy_hankel(
order: &BigFloat,
x: &BigFloat,
want_j: bool,
wp: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Option<BigFloat> {
let wp = wp + x.exponent().unwrap_or(0).max(0) as usize;
let mut xw = x.clone();
let _ = xw.set_precision(wp, rm);
let mut nu = order.clone();
let _ = nu.set_precision(wp, rm);
let (p_val, q_val, converged) = hankel_pq(&nu, &xw, wp, rm, cc);
if !converged {
tracing::trace!("bessel_jy_hankel: expansion did not reach target precision");
return None;
}
let pi = cc.pi(wp, rm).clone();
let two = BigFloat::from_i32(2, wp);
let four = BigFloat::from_i32(4, wp);
let two_over_pi_x = two.div(&pi.mul(&xw, wp, rm), wp, rm);
let amplitude = two_over_pi_x.sqrt(wp, rm);
let nu_pi_half = nu.mul(&pi, wp, rm).div(&two, wp, rm);
let pi_quarter = pi.div(&four, wp, rm);
let phase = xw.sub(&nu_pi_half, wp, rm).sub(&pi_quarter, wp, rm);
let cos_phase = phase.cos(wp, rm, cc);
let sin_phase = phase.sin(wp, rm, cc);
let combo = if want_j {
cos_phase
.mul(&p_val, wp, rm)
.sub(&sin_phase.mul(&q_val, wp, rm), wp, rm)
} else {
sin_phase
.mul(&p_val, wp, rm)
.add(&cos_phase.mul(&q_val, wp, rm), wp, rm)
};
Some(amplitude.mul(&combo, wp, rm))
}
fn arb_bessel_j(
order: &BigFloat,
x: &BigFloat,
prec: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
let guard = 32;
let wp = prec + guard;
if x.is_zero() {
let order_f64 = bigfloat_to_f64(order, rm, cc)?;
return if order_f64.abs() < 1e-15 {
Ok(BigFloat::from_i32(1, prec))
} else if order_f64 > 0.0 {
Ok(BigFloat::new(prec))
} else {
Err(SymplexError::Unevaluable {
reason: "BesselJ at x=0 with negative order".into(),
})
};
}
let x_f64 = bigfloat_to_f64(x, rm, cc)?;
let order_f64 = bigfloat_to_f64(order, rm, cc)?;
let order_int = order_f64.round() as i64;
let is_int_order = (order_f64 - order_int as f64).abs() < 1e-12;
if is_int_order && order_int < 0 {
let pos_order = BigFloat::from_i128((-order_int) as i128, wp);
let j = arb_bessel_j(&pos_order, x, prec, rm, cc)?;
return Ok(if order_int % 2 == 0 { j } else { j.neg() });
}
if x_f64.abs() >= hankel_threshold(wp, order_f64.abs()) {
if x.is_negative() {
if !is_int_order {
return Err(SymplexError::Unevaluable {
reason: "BesselJ of negative argument with non-integer order is complex".into(),
});
}
let j = arb_bessel_j(order, &x.abs(), prec, rm, cc)?;
return Ok(if order_int % 2 == 0 { j } else { j.neg() });
}
if let Some(j) = bessel_jy_hankel(order, x, true, wp, rm, cc) {
return Ok(round_to(j, prec, rm));
}
}
{
let wp = wp + cancellation_guard_bits(x_f64.abs());
let mut xw = x.clone();
let _ = xw.set_precision(wp, rm);
let mut order_w = order.clone();
let _ = order_w.set_precision(wp, rm);
let order = &order_w;
let two = BigFloat::from_i32(2, wp);
let x_half = xw.div(&two, wp, rm);
let x_half_sq = x_half.mul(&x_half, wp, rm);
let neg_x_half_sq = x_half_sq.neg();
let prefix = if is_int_order && order_int >= 0 {
let mut p = BigFloat::from_i32(1, wp);
for _ in 0..order_int {
p = p.mul(&x_half, wp, rm);
}
p
} else {
let ln_xh = x_half.abs().ln(wp, rm, cc);
let nu_ln = order.mul(&ln_xh, wp, rm);
nu_ln.exp(wp, rm, cc)
};
let mut sum = BigFloat::new(wp);
let gamma_nu1 = arb_gamma_real(&order.add(&BigFloat::from_i32(1, wp), wp, rm), wp, rm, cc)?;
let mut term = BigFloat::from_i32(1, wp).div(&gamma_nu1, wp, rm);
sum = sum.add(&term, wp, rm);
let max_terms = (x_f64.abs() * 1.5) as usize + wp + 40;
for k in 1..=max_terms {
let k_bf = BigFloat::from_i32(k as i32, wp);
let nu_plus_k = order.add(&k_bf, wp, rm);
let denom = k_bf.mul(&nu_plus_k, wp, rm);
term = term.mul(&neg_x_half_sq, wp, rm);
term = term.div(&denom, wp, rm);
if (k as f64) > x_f64.abs() / 2.0
&& let (Some(t_exp), Some(s_exp)) = (term.exponent(), sum.exponent())
&& (s_exp as i64 - t_exp as i64) > wp as i64
{
tracing::trace!(k, "arb_bessel_j: series converged");
break;
}
sum = sum.add(&term, wp, rm);
}
Ok(round_to(prefix.mul(&sum, wp, rm), prec, rm))
}
}
fn arb_bessel_y(
order: &BigFloat,
x: &BigFloat,
prec: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
let guard = 32;
let wp = prec + guard;
if x.is_zero() || x.is_negative() {
return Err(SymplexError::Unevaluable {
reason: "BesselY undefined at x ≤ 0".into(),
});
}
let x_f64 = bigfloat_to_f64(x, rm, cc)?;
let order_f64 = bigfloat_to_f64(order, rm, cc)?;
let order_int = order_f64.round() as i64;
let is_int_order = (order_f64 - order_int as f64).abs() < 1e-12;
if is_int_order && order_int < 0 {
let pos_order = BigFloat::from_i128((-order_int) as i128, wp);
let y = arb_bessel_y(&pos_order, x, prec, rm, cc)?;
return Ok(if order_int % 2 == 0 { y } else { y.neg() });
}
if x_f64 >= hankel_threshold(wp, order_f64.abs())
&& let Some(y) = bessel_jy_hankel(order, x, false, wp, rm, cc)
{
tracing::trace!("arb_bessel_y: used Hankel P/Q expansion");
return Ok(round_to(y, prec, rm));
}
if !is_int_order {
return Err(SymplexError::Unevaluable {
reason: format!(
"BesselY for non-integer order {order_f64} not yet implemented in small-|x| regime"
),
});
}
let n = order_int as usize;
let wp = wp + cancellation_guard_bits(x_f64);
let mut xw = x.clone();
let _ = xw.set_precision(wp, rm);
let pi = cc.pi(wp, rm).clone();
let one = BigFloat::from_i32(1, wp);
let two = BigFloat::from_i32(2, wp);
let inv_pi = one.div(&pi, wp, rm);
let x_half = xw.div(&two, wp, rm);
let x_half_sq = x_half.mul(&x_half, wp, rm);
let ln_x_half = x_half.ln(wp, rm, cc);
let euler_gamma = arb_euler_gamma(wp, rm, cc)?;
let n_bf = BigFloat::from_i128(n as i128, wp);
let j_n = arb_bessel_j(&n_bf, &xw, wp, rm, cc)?;
let mut finite = BigFloat::new(wp);
if n > 0 {
let mut f = Ratio::<BigInt>::from_integer(BigInt::from(1)); for i in 2..n {
f *= Ratio::from_integer(BigInt::from(i as u64));
}
let mut pow = one.clone();
for k in 0..n {
if k > 0 {
f /= Ratio::from_integer(BigInt::from(((n - k) * k) as u64));
pow = pow.mul(&x_half_sq, wp, rm);
}
let coeff = ratio_to_bigfloat(&f, wp, rm);
finite = finite.add(&coeff.mul(&pow, wp, rm), wp, rm);
}
let x_half_pow_neg_n = one.div(&x_half.powi(n, wp, rm), wp, rm);
finite = inv_pi
.mul(&x_half_pow_neg_n, wp, rm)
.mul(&finite, wp, rm)
.neg();
}
let log_term = two
.mul(&inv_pi, wp, rm)
.mul(&ln_x_half, wp, rm)
.mul(&j_n, wp, rm);
let mut h_k = BigFloat::new(wp);
let mut h_nk = BigFloat::new(wp);
for m in 1..=n {
h_nk = h_nk.add(
&one.div(&BigFloat::from_i128(m as i128, wp), wp, rm),
wp,
rm,
);
}
let mut n_fact = Ratio::<BigInt>::from_integer(BigInt::from(1));
for i in 2..=n {
n_fact *= Ratio::from_integer(BigInt::from(i as u64));
}
let mut inv_fact = one.div(&ratio_to_bigfloat(&n_fact, wp, rm), wp, rm);
let neg_x_half_sq = x_half_sq.neg();
let mut pow = one.clone();
let mut series = BigFloat::new(wp);
let two_gamma = euler_gamma.mul(&two, wp, rm);
let max_terms = (x_f64.abs() * 2.0) as usize + wp + 40;
for k in 0..=max_terms {
if k > 0 {
let k_bf = BigFloat::from_i128(k as i128, wp);
let nk_bf = BigFloat::from_i128((n + k) as i128, wp);
h_k = h_k.add(&one.div(&k_bf, wp, rm), wp, rm);
h_nk = h_nk.add(&one.div(&nk_bf, wp, rm), wp, rm);
inv_fact = inv_fact.div(&k_bf.mul(&nk_bf, wp, rm), wp, rm);
pow = pow.mul(&neg_x_half_sq, wp, rm);
}
let psi_sum = h_k.add(&h_nk, wp, rm).sub(&two_gamma, wp, rm);
let term = psi_sum.mul(&pow, wp, rm).mul(&inv_fact, wp, rm);
series = series.add(&term, wp, rm);
if k > 2
&& let (Some(t_exp), Some(s_exp)) = (term.exponent(), series.exponent())
&& (s_exp as i64 - t_exp as i64) > wp as i64
{
tracing::trace!(k, "arb_bessel_y: Neumann series converged");
break;
}
}
let series_term = inv_pi
.mul(&x_half.powi(n, wp, rm), wp, rm)
.mul(&series, wp, rm)
.neg();
let y = finite.add(&log_term, wp, rm).add(&series_term, wp, rm);
Ok(round_to(y, prec, rm))
}
fn arb_catalan(prec: usize, rm: RoundingMode, cc: &mut Consts) -> Result<BigFloat, SymplexError> {
let wp = prec + 32;
let one = BigFloat::from_i32(1, wp);
let two = BigFloat::from_i32(2, wp);
let three = BigFloat::from_i32(3, wp);
let eight = BigFloat::from_i32(8, wp);
let sqrt3 = three.sqrt(wp, rm);
let ln_term = two.add(&sqrt3, wp, rm).ln(wp, rm, cc);
let pi = cc.pi(wp, rm).clone();
let first = pi.div(&eight, wp, rm).mul(&ln_term, wp, rm);
let mut b = one.clone();
let mut sum = BigFloat::new(wp);
let max_terms = wp / 2 + 64;
for n in 0..max_terms {
let two_n_plus_1 = BigFloat::from_i128(2 * n as i128 + 1, wp);
let denom = two_n_plus_1.mul(&two_n_plus_1, wp, rm);
let term = b.div(&denom, wp, rm);
sum = sum.add(&term, wp, rm);
if let (Some(t_exp), Some(s_exp)) = (term.exponent(), sum.exponent())
&& (s_exp as i64 - t_exp as i64) > wp as i64
{
break;
}
let n_plus_1 = BigFloat::from_i128(n as i128 + 1, wp);
let two_2n_plus_1 = two.mul(&two_n_plus_1, wp, rm);
b = b.mul(&n_plus_1, wp, rm).div(&two_2n_plus_1, wp, rm);
}
let second = three.div(&eight, wp, rm).mul(&sum, wp, rm);
let mut g = first.add(&second, wp, rm);
g.set_precision(prec, rm)
.map_err(|e| SymplexError::ComputationFailed {
operation: "catalan",
reason: format!("precision adjustment failed: {e:?}"),
})?;
Ok(g)
}
fn cancellation_guard_bits(x_abs: f64) -> usize {
(x_abs * std::f64::consts::LOG2_E).ceil() as usize + 32
}
fn si_ci_asymptotic_fg(x: &BigFloat, wp: usize, rm: RoundingMode) -> (BigFloat, BigFloat) {
let one = BigFloat::from_i32(1, wp);
let inv_x = one.div(x, wp, rm);
let inv_x2 = inv_x.mul(&inv_x, wp, rm);
let mut f_sum = BigFloat::new(wp);
let mut g_sum = BigFloat::new(wp);
let mut t_f = one.clone();
let mut t_g = one.clone();
let mut prev_f: Option<BigFloat> = None;
let max_terms = wp * 2 + 100;
for k in 0..max_terms {
if k > 0 {
let a = BigFloat::from_i128((2 * k as i128 - 1) * (2 * k as i128), wp);
let b = BigFloat::from_i128((2 * k as i128) * (2 * k as i128 + 1), wp);
t_f = t_f.mul(&a, wp, rm).mul(&inv_x2, wp, rm);
t_g = t_g.mul(&b, wp, rm).mul(&inv_x2, wp, rm);
}
if let Some(ref p) = prev_f
&& t_f.abs().cmp(p).is_some_and(|c| c > 0)
{
break;
}
if k % 2 == 0 {
f_sum = f_sum.add(&t_f, wp, rm);
g_sum = g_sum.add(&t_g, wp, rm);
} else {
f_sum = f_sum.sub(&t_f, wp, rm);
g_sum = g_sum.sub(&t_g, wp, rm);
}
if let (Some(t_exp), Some(s_exp)) = (t_f.exponent(), f_sum.exponent())
&& (s_exp as i64 - t_exp as i64) > wp as i64
{
break;
}
prev_f = Some(t_f.abs());
}
(f_sum.mul(&inv_x, wp, rm), g_sum.mul(&inv_x2, wp, rm))
}
fn arb_si_ci(
x: &BigFloat,
want_si: bool,
prec: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<(BigFloat, BigFloat), SymplexError> {
if x.is_nan() || x.is_inf_pos() || x.is_inf_neg() {
return Err(SymplexError::Unevaluable {
reason: "Si/Ci of special float value".into(),
});
}
let x_f64 = bigfloat_to_f64(x, rm, cc)?;
let x_abs = x_f64.abs();
let base_wp = prec + 32;
if x.is_zero() {
return if want_si {
Ok((BigFloat::new(prec), BigFloat::new(prec)))
} else {
Err(SymplexError::Unevaluable {
reason: "Ci(0) is -∞".into(),
})
};
}
if !want_si && x.is_negative() {
return Err(SymplexError::Unevaluable {
reason: "Ci requires a positive argument here".into(),
});
}
let switch = (base_wp as f64) * std::f64::consts::LN_2;
if x_abs > switch {
let wp = base_wp;
let ax = x.abs();
let (f, g) = si_ci_asymptotic_fg(&ax, wp, rm);
let sin_x = ax.sin(wp, rm, cc);
let cos_x = ax.cos(wp, rm, cc);
let pi = cc.pi(wp, rm).clone();
let two = BigFloat::from_i32(2, wp);
let half_pi = pi.div(&two, wp, rm);
let mut si =
half_pi
.sub(&f.mul(&cos_x, wp, rm), wp, rm)
.sub(&g.mul(&sin_x, wp, rm), wp, rm);
if x.is_negative() {
si = si.neg();
}
let ci = f.mul(&sin_x, wp, rm).sub(&g.mul(&cos_x, wp, rm), wp, rm);
return Ok((round_to(si, prec, rm), round_to(ci, prec, rm)));
}
let wp = base_wp + cancellation_guard_bits(x_abs);
let mut xw = x.clone();
let _ = xw.set_precision(wp, rm);
let x2 = xw.mul(&xw, wp, rm);
let max_terms = (x_abs * 1.5) as usize + wp + 40;
let si = if want_si {
let mut t = xw.clone();
let mut sum = BigFloat::new(wp);
for k in 0..max_terms {
if k > 0 {
let d = BigFloat::from_i128((2 * k as i128) * (2 * k as i128 + 1), wp);
t = t.mul(&x2, wp, rm).div(&d, wp, rm);
}
let term = t.div(&BigFloat::from_i128(2 * k as i128 + 1, wp), wp, rm);
if k % 2 == 0 {
sum = sum.add(&term, wp, rm);
} else {
sum = sum.sub(&term, wp, rm);
}
if let (Some(t_exp), Some(s_exp)) = (term.exponent(), sum.exponent())
&& (s_exp as i64 - t_exp as i64) > wp as i64
{
break;
}
}
round_to(sum, prec, rm)
} else {
BigFloat::new(prec)
};
let ci = if want_si {
BigFloat::new(prec)
} else {
let gamma = arb_euler_gamma(wp, rm, cc)?;
let ln_x = xw.ln(wp, rm, cc);
let mut sum = BigFloat::new(wp);
let mut u = BigFloat::from_i32(1, wp);
for k in 1..max_terms {
let d = BigFloat::from_i128((2 * k as i128 - 1) * (2 * k as i128), wp);
u = u.mul(&x2, wp, rm).div(&d, wp, rm);
let term = u.div(&BigFloat::from_i128(2 * k as i128, wp), wp, rm);
if k % 2 == 0 {
sum = sum.add(&term, wp, rm);
} else {
sum = sum.sub(&term, wp, rm);
}
if let (Some(t_exp), Some(s_exp)) = (term.exponent(), sum.exponent())
&& (s_exp as i64 - t_exp as i64) > wp as i64
{
break;
}
}
let ci = gamma.add(&ln_x, wp, rm).add(&sum, wp, rm);
round_to(ci, prec, rm)
};
Ok((si, ci))
}
fn arb_ei(
x: &BigFloat,
prec: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
if x.is_nan() || x.is_inf_pos() || x.is_inf_neg() || x.is_zero() {
return Err(SymplexError::Unevaluable {
reason: "Ei of zero or special float value".into(),
});
}
let x_f64 = bigfloat_to_f64(x, rm, cc)?;
let x_abs = x_f64.abs();
let base_wp = prec + 32;
let switch = (base_wp as f64) * std::f64::consts::LN_2;
if x_abs > switch {
let wp = base_wp;
let one = BigFloat::from_i32(1, wp);
let inv_x = one.div(x, wp, rm);
let mut term = one.clone();
let mut sum = BigFloat::new(wp);
let mut prev: Option<BigFloat> = None;
let max_terms = wp * 2 + 100;
for k in 0..max_terms {
if k > 0 {
term = term
.mul(&BigFloat::from_i128(k as i128, wp), wp, rm)
.mul(&inv_x, wp, rm);
}
if let Some(ref p) = prev
&& term.abs().cmp(p).is_some_and(|c| c > 0)
{
break;
}
sum = sum.add(&term, wp, rm);
if let (Some(t_exp), Some(s_exp)) = (term.exponent(), sum.exponent())
&& (s_exp as i64 - t_exp as i64) > wp as i64
{
break;
}
prev = Some(term.abs());
}
let ex = x.exp(wp, rm, cc);
let r = ex.mul(&inv_x, wp, rm).mul(&sum, wp, rm);
return Ok(round_to(r, prec, rm));
}
let wp = if x.is_negative() {
base_wp + cancellation_guard_bits(x_abs)
} else {
base_wp
};
let mut xw = x.clone();
let _ = xw.set_precision(wp, rm);
let gamma = arb_euler_gamma(wp, rm, cc)?;
let ln_abs_x = xw.abs().ln(wp, rm, cc);
let mut v = BigFloat::from_i32(1, wp); let mut sum = BigFloat::new(wp);
let max_terms = (x_abs * 3.0) as usize + wp + 40;
for k in 1..max_terms {
let k_bf = BigFloat::from_i128(k as i128, wp);
v = v.mul(&xw, wp, rm).div(&k_bf, wp, rm);
let term = v.div(&k_bf, wp, rm);
sum = sum.add(&term, wp, rm);
if let (Some(t_exp), Some(s_exp)) = (term.exponent(), sum.exponent())
&& (s_exp as i64 - t_exp as i64) > wp as i64
&& (k as f64) > x_abs
{
break;
}
}
let r = gamma.add(&ln_abs_x, wp, rm).add(&sum, wp, rm);
Ok(round_to(r, prec, rm))
}
fn round_to(mut v: BigFloat, prec: usize, rm: RoundingMode) -> BigFloat {
let _ = v.set_precision(prec, rm);
v
}
fn bf_pow(b: &BigFloat, e: &BigFloat, wp: usize, rm: RoundingMode, cc: &mut Consts) -> BigFloat {
if e.is_zero() {
return BigFloat::from_i32(1, wp);
}
if e.is_int()
&& let Ok(ef) = bigfloat_to_f64(e, rm, cc)
&& ef.abs() < 1.0e9
{
let n = ef.abs() as usize;
let p = b.powi(n, wp, rm);
return if ef < 0.0 {
BigFloat::from_i32(1, wp).div(&p, wp, rm)
} else {
p
};
}
let ln_b = b.ln(wp + 32, rm, cc);
ln_b.mul(e, wp + 32, rm).exp(wp, rm, cc)
}
fn arb_zeta(
s: &BigFloat,
prec: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
if s.is_nan() || s.is_inf_neg() {
return Err(SymplexError::Unevaluable {
reason: "zeta of special float value".into(),
});
}
if s.is_inf_pos() {
return Ok(BigFloat::from_i32(1, prec));
}
if s.is_zero() {
return Ok(BigFloat::from_f64(-0.5, prec));
}
let s_f64 = bigfloat_to_f64(s, rm, cc)?;
if (s_f64 - 1.0).abs() < 1e-300 {
return Err(SymplexError::Unevaluable {
reason: "zeta(1) is a pole".into(),
});
}
let wp = prec + 32;
if s.is_negative() {
let s_round = s_f64.round();
if (s_f64 - s_round).abs() < 1e-14 && (s_round as i64) % 2 == 0 {
let mut sw = s.clone();
let _ = sw.set_precision(wp, rm);
let diff = sw.sub(&BigFloat::from_f64(s_round, wp), wp, rm);
if diff.is_zero() {
return Ok(BigFloat::new(prec));
}
}
let mut sw = s.clone();
let _ = sw.set_precision(wp, rm);
let one = BigFloat::from_i32(1, wp);
let two = BigFloat::from_i32(2, wp);
let pi = cc.pi(wp, rm).clone();
let one_minus_s = one.sub(&sw, wp, rm);
let two_pow_s = bf_pow(&two, &sw, wp, rm, cc);
let s_minus_1 = sw.sub(&one, wp, rm);
let pi_pow = bf_pow(&pi, &s_minus_1, wp, rm, cc);
let half_pi_s = pi.mul(&sw, wp, rm).div(&two, wp, rm);
let sin_term = half_pi_s.sin(wp, rm, cc);
let gamma_term = arb_gamma_real(&one_minus_s, wp, rm, cc)?;
let zeta_term = arb_zeta_borwein(&one_minus_s, wp, rm, cc)?;
let r = two_pow_s
.mul(&pi_pow, wp, rm)
.mul(&sin_term, wp, rm)
.mul(&gamma_term, wp, rm)
.mul(&zeta_term, wp, rm);
return Ok(round_to(r, prec, rm));
}
let r = arb_zeta_borwein(s, wp, rm, cc)?;
Ok(round_to(r, prec, rm))
}
fn arb_zeta_borwein(
s: &BigFloat,
wp: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
let eta = arb_eta_borwein(s, wp, rm, cc)?;
let one = BigFloat::from_i32(1, wp);
let two = BigFloat::from_i32(2, wp);
let one_minus_s = one.sub(s, wp, rm);
let two_pow = bf_pow(&two, &one_minus_s, wp, rm, cc);
let denom_factor = one.sub(&two_pow, wp, rm);
if denom_factor.is_zero() {
return Err(SymplexError::Unevaluable {
reason: "zeta(1) is a pole".into(),
});
}
Ok(eta.div(&denom_factor, wp, rm))
}
fn arb_eta_borwein(
s: &BigFloat,
wp: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
let n = ((wp as f64) * std::f64::consts::LN_2 / (3.0 + 8f64.sqrt()).ln()).ceil() as usize + 4;
let mut d: Vec<Ratio<BigInt>> = Vec::with_capacity(n + 1);
let mut acc = Ratio::<BigInt>::zero();
let mut term = Ratio::from_integer(BigInt::from(1));
for i in 0..=n {
if i > 0 {
let numer = BigInt::from((n + i - 1) as u64)
* BigInt::from(4u64)
* BigInt::from((n - i + 1) as u64);
let denom = BigInt::from((2 * i - 1) as u64) * BigInt::from((2 * i) as u64);
term *= Ratio::new(numer, denom);
}
acc += &term;
d.push(acc.clone());
}
let d_n = ratio_to_bigfloat(&d[n], wp, rm);
let mut sum = BigFloat::new(wp);
let s_is_int = s.is_int();
let s_int = if s_is_int {
bigfloat_to_f64(s, rm, cc)?.round() as usize
} else {
0
};
for k in 0..n {
let coeff = ratio_to_bigfloat(&(&d[k] - &d[n]), wp, rm);
let kp1 = BigFloat::from_i128(k as i128 + 1, wp);
let kp1_pow_s = if s_is_int {
kp1.powi(s_int, wp, rm)
} else {
bf_pow(&kp1, s, wp, rm, cc)
};
let term = coeff.div(&kp1_pow_s, wp, rm);
if k % 2 == 0 {
sum = sum.add(&term, wp, rm);
} else {
sum = sum.sub(&term, wp, rm);
}
}
Ok(sum.div(&d_n, wp, rm).neg())
}
fn arb_polygamma(
n: u32,
x: &BigFloat,
prec: usize,
rm: RoundingMode,
_cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
if x.is_nan() || x.is_inf_pos() || x.is_inf_neg() {
return Err(SymplexError::Unevaluable {
reason: "polygamma of special float value".into(),
});
}
let guard = 32;
let wp = prec + guard;
let n_us = n as usize;
let one = BigFloat::from_i32(1, wp);
let mut n_fact = Ratio::<BigInt>::from_integer(BigInt::from(1));
for i in 2..=n_us {
n_fact *= Ratio::from_integer(BigInt::from(i as u64));
}
let nm1_fact = &n_fact / Ratio::from_integer(BigInt::from(n_us as u64));
let n_fact_bf = ratio_to_bigfloat(&n_fact, wp, rm);
let nm1_fact_bf = ratio_to_bigfloat(&nm1_fact, wp, rm);
let outer_sign_neg = n_us.is_multiple_of(2);
let threshold_val = ((wp / 3).max(n_us) + 10) as i128;
let threshold = BigFloat::from_i128(threshold_val, wp);
let mut xw = x.clone();
let _ = xw.set_precision(wp, rm);
let mut shift_sum = BigFloat::new(wp); let mut steps: usize = 0;
while xw.sub(&threshold, wp, rm).is_negative() {
if xw.is_zero() {
return Err(SymplexError::Unevaluable {
reason: "polygamma at non-positive integer pole".into(),
});
}
if let Some(x_exp) = xw.exponent()
&& (x_exp as i64) < -(wp as i64 / 2)
{
return Err(SymplexError::Unevaluable {
reason: "polygamma at non-positive integer pole".into(),
});
}
let p = xw.powi(n_us + 1, wp, rm);
let inv = one.div(&p, wp, rm);
shift_sum = shift_sum.add(&inv, wp, rm);
xw = xw.add(&one, wp, rm);
steps += 1;
if steps > 10 * wp + 1_000_000 {
return Err(SymplexError::ComputationFailed {
operation: "polygamma",
reason: "argument shift did not terminate".into(),
});
}
}
let shift_term = n_fact_bf.mul(&shift_sum, wp, rm);
let x_pow_n = xw.powi(n_us, wp, rm);
let x_pow_np1 = x_pow_n.mul(&xw, wp, rm);
let two = BigFloat::from_i32(2, wp);
let mut bracket = nm1_fact_bf.div(&x_pow_n, wp, rm);
bracket = bracket.add(&n_fact_bf.div(&two.mul(&x_pow_np1, wp, rm), wp, rm), wp, rm);
let x2 = xw.mul(&xw, wp, rm);
let mut x_pow = x_pow_n.mul(&x2, wp, rm); let mut ratio_fact = &n_fact * Ratio::from_integer(BigInt::from(n_us as u64 + 1))
/ Ratio::from_integer(BigInt::from(2));
let mut prev_term_exp: Option<i32> = None;
let max_terms = wp / 2 + 20;
for k in 1..=max_terms {
if k > 1 {
let a =
BigInt::from((2 * k + n_us - 2) as u64) * BigInt::from((2 * k + n_us - 1) as u64);
let b = BigInt::from((2 * k - 1) as u64) * BigInt::from((2 * k) as u64);
ratio_fact *= Ratio::new(a, b);
x_pow = x_pow.mul(&x2, wp, rm);
}
let b2k = crate::base::bernoulli::bernoulli(2 * k);
let coeff = ratio_to_bigfloat(&(&b2k * &ratio_fact), wp, rm);
let term = coeff.div(&x_pow, wp, rm);
if let Some(t_exp) = term.exponent() {
if let Some(prev) = prev_term_exp
&& t_exp > prev + 10
{
break; }
prev_term_exp = Some(t_exp);
if let Some(b_exp) = bracket.exponent()
&& (b_exp as i64 - t_exp as i64) > wp as i64
{
bracket = bracket.add(&term, wp, rm);
break;
}
}
bracket = bracket.add(&term, wp, rm);
}
let total = bracket.add(&shift_term, wp, rm);
let result = if outer_sign_neg { total.neg() } else { total };
Ok(round_to(result, prec, rm))
}
fn bessel_i_series(
order: &BigFloat,
x: &BigFloat,
wp: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
let two = BigFloat::from_i32(2, wp);
let x_half = x.div(&two, wp, rm);
let x_half_sq = x_half.mul(&x_half, wp, rm);
let order_f64 = bigfloat_to_f64(order, rm, cc)?;
let order_int = order_f64.round() as i64;
let is_int_order = (order_f64 - order_int as f64).abs() < 1e-12;
if x.is_zero() {
return if is_int_order && order_int == 0 {
Ok(BigFloat::from_i32(1, wp))
} else if order_f64 > 0.0 {
Ok(BigFloat::new(wp))
} else {
Err(SymplexError::Unevaluable {
reason: "BesselI at x=0 with negative order".into(),
})
};
}
let prefix = if is_int_order && order_int >= 0 {
x_half.powi(order_int as usize, wp, rm)
} else if is_int_order {
x_half.powi((-order_int) as usize, wp, rm)
} else {
if x.is_negative() {
return Err(SymplexError::Unevaluable {
reason: "BesselI of negative argument with non-integer order is complex".into(),
});
}
let ln_xh = x_half.ln(wp, rm, cc);
order.mul(&ln_xh, wp, rm).exp(wp, rm, cc)
};
let eff_order = if is_int_order && order_int < 0 {
BigFloat::from_i128((-order_int) as i128, wp)
} else {
order.clone()
};
let one = BigFloat::from_i32(1, wp);
let gamma_nu1 = arb_gamma_real(&eff_order.add(&one, wp, rm), wp, rm, cc)?;
let mut term = one.div(&gamma_nu1, wp, rm);
let mut sum = term.clone();
let x_f64 = bigfloat_to_f64(x, rm, cc)?.abs();
let max_terms = (x_f64 * 2.0) as usize + wp + 40;
for k in 1..=max_terms {
let k_bf = BigFloat::from_i128(k as i128, wp);
let nu_plus_k = eff_order.add(&k_bf, wp, rm);
let denom = k_bf.mul(&nu_plus_k, wp, rm);
term = term.mul(&x_half_sq, wp, rm).div(&denom, wp, rm);
sum = sum.add(&term, wp, rm);
if let (Some(t_exp), Some(s_exp)) = (term.exponent(), sum.exponent())
&& (s_exp as i64 - t_exp as i64) > wp as i64
{
break;
}
}
Ok(prefix.mul(&sum, wp, rm))
}
fn arb_bessel_i(
order: &BigFloat,
x: &BigFloat,
prec: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
let wp = prec + 32;
let r = bessel_i_series(order, x, wp, rm, cc)?;
Ok(round_to(r, prec, rm))
}
fn bessel_k_asymptotic(
order: &BigFloat,
x: &BigFloat,
wp: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> BigFloat {
let one = BigFloat::from_i32(1, wp);
let four = BigFloat::from_i32(4, wp);
let eight = BigFloat::from_i32(8, wp);
let four_nu_sq = four.mul(&order.mul(order, wp, rm), wp, rm);
let inv_x = one.div(x, wp, rm);
let mut a_k = one.clone();
let mut sum = one.clone();
let mut prev: Option<BigFloat> = None;
let max_terms = wp * 2 + 100;
for k in 1..max_terms {
let two_km1 = BigFloat::from_i32((2 * k as i32) - 1, wp);
let numer = four_nu_sq.sub(&two_km1.mul(&two_km1, wp, rm), wp, rm);
let denom = eight.mul(&BigFloat::from_i32(k as i32, wp), wp, rm);
a_k = a_k
.mul(&numer, wp, rm)
.div(&denom, wp, rm)
.mul(&inv_x, wp, rm);
if let Some(ref p) = prev
&& a_k.abs().cmp(p).is_some_and(|c| c > 0)
{
break;
}
sum = sum.add(&a_k, wp, rm);
if let (Some(t_exp), Some(s_exp)) = (a_k.exponent(), sum.exponent())
&& (s_exp as i64 - t_exp as i64) > wp as i64
{
break;
}
prev = Some(a_k.abs());
}
let pi = cc.pi(wp, rm).clone();
let two = BigFloat::from_i32(2, wp);
let amp = pi.div(&two.mul(x, wp, rm), wp, rm).sqrt(wp, rm);
let e_neg_x = x.neg().exp(wp, rm, cc);
amp.mul(&e_neg_x, wp, rm).mul(&sum, wp, rm)
}
fn arb_bessel_k(
order: &BigFloat,
x: &BigFloat,
prec: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
if x.is_zero() || x.is_negative() || x.is_nan() {
return Err(SymplexError::Unevaluable {
reason: "BesselK undefined at x ≤ 0".into(),
});
}
let base_wp = prec + 32;
let x_f64 = bigfloat_to_f64(x, rm, cc)?;
let order_f64 = bigfloat_to_f64(order, rm, cc)?;
let order_int = order_f64.round() as i64;
let is_int_order = (order_f64 - order_int as f64).abs() < 1e-12;
let order_abs = if is_int_order {
BigFloat::from_i128(order_int.unsigned_abs() as i128, base_wp)
} else {
order.abs()
};
if 2.0 * x_f64 > (base_wp as f64) * std::f64::consts::LN_2 {
let r = bessel_k_asymptotic(&order_abs, x, base_wp, rm, cc);
return Ok(round_to(r, prec, rm));
}
let wp = base_wp + 2 * cancellation_guard_bits(x_f64);
let mut xw = x.clone();
let _ = xw.set_precision(wp, rm);
if !is_int_order {
let neg_order = order.neg();
let i_neg = bessel_i_series(&neg_order, &xw, wp, rm, cc)?;
let i_pos = bessel_i_series(order, &xw, wp, rm, cc)?;
let pi = cc.pi(wp, rm).clone();
let two = BigFloat::from_i32(2, wp);
let sin_nu_pi = order.mul(&pi, wp, rm).sin(wp, rm, cc);
let r = pi
.div(&two, wp, rm)
.mul(&i_neg.sub(&i_pos, wp, rm), wp, rm)
.div(&sin_nu_pi, wp, rm);
return Ok(round_to(r, prec, rm));
}
let n = order_int.unsigned_abs() as usize;
let one = BigFloat::from_i32(1, wp);
let two = BigFloat::from_i32(2, wp);
let half = one.div(&two, wp, rm);
let x_half = xw.div(&two, wp, rm);
let x_half_sq = x_half.mul(&x_half, wp, rm);
let ln_x_half = x_half.ln(wp, rm, cc);
let gamma = arb_euler_gamma(wp, rm, cc)?;
let n_bf = BigFloat::from_i128(n as i128, wp);
let i_n = bessel_i_series(&n_bf, &xw, wp, rm, cc)?;
let mut finite = BigFloat::new(wp);
if n > 0 {
let mut f = Ratio::<BigInt>::from_integer(BigInt::from(1));
for i in 2..n {
f *= Ratio::from_integer(BigInt::from(i as u64));
}
let mut pow = one.clone(); let neg_x_half_sq = x_half_sq.neg();
for k in 0..n {
if k > 0 {
f /= Ratio::from_integer(BigInt::from(((n - k) * k) as u64));
pow = pow.mul(&neg_x_half_sq, wp, rm);
}
let coeff = ratio_to_bigfloat(&f, wp, rm);
finite = finite.add(&coeff.mul(&pow, wp, rm), wp, rm);
}
let x_half_pow_neg_n = one.div(&x_half.powi(n, wp, rm), wp, rm);
finite = half.mul(&x_half_pow_neg_n, wp, rm).mul(&finite, wp, rm);
}
let log_term = ln_x_half.mul(&i_n, wp, rm);
let log_term = if n.is_multiple_of(2) {
log_term.neg()
} else {
log_term
};
let mut h_k = BigFloat::new(wp); let mut h_nk = BigFloat::new(wp); for m in 1..=n {
h_nk = h_nk.add(
&one.div(&BigFloat::from_i128(m as i128, wp), wp, rm),
wp,
rm,
);
}
let mut n_fact = Ratio::<BigInt>::from_integer(BigInt::from(1));
for i in 2..=n {
n_fact *= Ratio::from_integer(BigInt::from(i as u64));
}
let mut inv_fact = one.div(&ratio_to_bigfloat(&n_fact, wp, rm), wp, rm);
let mut pow = one.clone();
let mut series = BigFloat::new(wp);
let two_gamma = gamma.mul(&two, wp, rm);
let max_terms = (x_f64 * 2.0) as usize + wp + 40;
for k in 0..=max_terms {
if k > 0 {
let k_bf = BigFloat::from_i128(k as i128, wp);
let nk_bf = BigFloat::from_i128((n + k) as i128, wp);
h_k = h_k.add(&one.div(&k_bf, wp, rm), wp, rm);
h_nk = h_nk.add(&one.div(&nk_bf, wp, rm), wp, rm);
inv_fact = inv_fact.div(&k_bf.mul(&nk_bf, wp, rm), wp, rm);
pow = pow.mul(&x_half_sq, wp, rm);
}
let psi_sum = h_k.add(&h_nk, wp, rm).sub(&two_gamma, wp, rm);
let term = psi_sum.mul(&pow, wp, rm).mul(&inv_fact, wp, rm);
series = series.add(&term, wp, rm);
if k > 2
&& let (Some(t_exp), Some(s_exp)) = (term.exponent(), series.exponent())
&& (s_exp as i64 - t_exp as i64) > wp as i64
{
break;
}
}
let series_term = half
.mul(&x_half.powi(n, wp, rm), wp, rm)
.mul(&series, wp, rm);
let series_term = if n.is_multiple_of(2) {
series_term
} else {
series_term.neg()
};
let r = finite.add(&log_term, wp, rm).add(&series_term, wp, rm);
Ok(round_to(r, prec, rm))
}
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
enum OrthoPoly {
Legendre,
ChebyshevT,
ChebyshevU,
Hermite,
Laguerre,
}
fn arb_orthopoly(kind: OrthoPoly, n: u64, x: &BigFloat, prec: usize, rm: RoundingMode) -> BigFloat {
let wp = prec + 32 + (n as f64).log2().ceil().max(0.0) as usize;
let one = BigFloat::from_i32(1, wp);
let two = BigFloat::from_i32(2, wp);
let mut xw = x.clone();
let _ = xw.set_precision(wp, rm);
if n == 0 {
return BigFloat::from_i32(1, prec);
}
let mut p_prev = one.clone();
let mut p_curr = match kind {
OrthoPoly::Legendre | OrthoPoly::ChebyshevT => xw.clone(),
OrthoPoly::ChebyshevU | OrthoPoly::Hermite => two.mul(&xw, wp, rm),
OrthoPoly::Laguerre => one.sub(&xw, wp, rm),
};
for k in 1..n {
let k_bf = BigFloat::from_i128(k as i128, wp);
let kp1 = BigFloat::from_i128(k as i128 + 1, wp);
let next = match kind {
OrthoPoly::Legendre => {
let two_k_plus_1 = BigFloat::from_i128(2 * k as i128 + 1, wp);
let a = two_k_plus_1.mul(&xw, wp, rm).mul(&p_curr, wp, rm);
let b = k_bf.mul(&p_prev, wp, rm);
a.sub(&b, wp, rm).div(&kp1, wp, rm)
}
OrthoPoly::ChebyshevT | OrthoPoly::ChebyshevU => two
.mul(&xw, wp, rm)
.mul(&p_curr, wp, rm)
.sub(&p_prev, wp, rm),
OrthoPoly::Hermite => {
let a = two.mul(&xw, wp, rm).mul(&p_curr, wp, rm);
let b = two.mul(&k_bf, wp, rm).mul(&p_prev, wp, rm);
a.sub(&b, wp, rm)
}
OrthoPoly::Laguerre => {
let two_k_plus_1 = BigFloat::from_i128(2 * k as i128 + 1, wp);
let a = two_k_plus_1.sub(&xw, wp, rm).mul(&p_curr, wp, rm);
let b = k_bf.mul(&p_prev, wp, rm);
a.sub(&b, wp, rm).div(&kp1, wp, rm)
}
};
p_prev = p_curr;
p_curr = next;
}
round_to(p_curr, prec, rm)
}
fn unevaluable(reason: impl Into<String>) -> SymplexError {
SymplexError::Unevaluable {
reason: reason.into(),
}
}
fn series_did_not_converge(operation: &'static str, max_terms: usize) -> SymplexError {
SymplexError::ComputationFailed {
operation,
reason: format!("series did not converge in {max_terms} terms"),
}
}
fn real_args(
name: &str,
args: &[ExprId],
cache: &FxHashMap<ExprId, Complex>,
) -> Result<Vec<BigFloat>, SymplexError> {
let mut out = Vec::with_capacity(args.len());
for &a in args {
let v = get_cached(cache, a)?;
if !v.1.is_zero() {
return Err(unevaluable(format!(
"{name} of complex argument not yet supported in evalf"
)));
}
out.push(v.0.clone());
}
Ok(out)
}
fn negligible(a: &BigFloat, b: &BigFloat, bits: usize) -> bool {
if a.is_zero() {
return true;
}
match (a.exponent(), b.exponent()) {
(Some(ae), Some(be)) => (be as i64 - ae as i64) > bits as i64,
_ => false,
}
}
fn bf_gt(a: &BigFloat, b: &BigFloat) -> bool {
a.cmp(b).is_some_and(|c| c > 0)
}
fn bf_lt(a: &BigFloat, b: &BigFloat) -> bool {
a.cmp(b).is_some_and(|c| c < 0)
}
fn bf_as_int(x: &BigFloat, rm: RoundingMode, cc: &mut Consts) -> Result<Option<i64>, SymplexError> {
let f = bigfloat_to_f64(x, rm, cc)?;
let r = f.round();
if (f - r).abs() < 1e-12 && r.abs() < 9.0e15 {
Ok(Some(r as i64))
} else {
Ok(None)
}
}
fn eval_lib_fn(
f: LibFn,
args: &[ExprId],
cache: &FxHashMap<ExprId, Complex>,
prec: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<Complex, SymplexError> {
match f {
LibFn::BesselJ => eval_bessel(arb_bessel_j, "BesselJ", args, cache, prec, rm, cc),
LibFn::BesselY => eval_bessel(arb_bessel_y, "BesselY", args, cache, prec, rm, cc),
LibFn::BesselI => eval_bessel(arb_bessel_i, "BesselI", args, cache, prec, rm, cc),
LibFn::BesselK => eval_bessel(arb_bessel_k, "BesselK", args, cache, prec, rm, cc),
LibFn::Legendre => eval_orthopoly(OrthoPoly::Legendre, f, args, cache, prec, rm, cc),
LibFn::ChebyshevT => eval_orthopoly(OrthoPoly::ChebyshevT, f, args, cache, prec, rm, cc),
LibFn::ChebyshevU => eval_orthopoly(OrthoPoly::ChebyshevU, f, args, cache, prec, rm, cc),
LibFn::Hermite => eval_orthopoly(OrthoPoly::Hermite, f, args, cache, prec, rm, cc),
LibFn::Laguerre => eval_orthopoly(OrthoPoly::Laguerre, f, args, cache, prec, rm, cc),
LibFn::Erfi
| LibFn::ErfInv
| LibFn::ErfcInv
| LibFn::ExpInt
| LibFn::Shi
| LibFn::Chi
| LibFn::FresnelS
| LibFn::FresnelC
| LibFn::LowerGamma
| LibFn::UpperGamma
| LibFn::PolyLog
| LibFn::DirichletEta
| LibFn::AiryAi
| LibFn::AiryBi
| LibFn::AiryAiPrime
| LibFn::AiryBiPrime
| LibFn::EllipticK
| LibFn::EllipticE
| LibFn::EllipticF
| LibFn::EllipticPi
| LibFn::Gegenbauer
| LibFn::Jacobi
| LibFn::AssocLegendre
| LibFn::AssocLaguerre
| LibFn::BetaInc
| LibFn::BetaIncRegularized => eval_special_09(f, args, cache, prec, rm, cc),
LibFn::Factorial2
| LibFn::Subfactorial
| LibFn::RisingFactorial
| LibFn::FallingFactorial
| LibFn::Fibonacci
| LibFn::Lucas
| LibFn::Bernoulli
| LibFn::Harmonic
| LibFn::Catalan
| LibFn::Bell
| LibFn::EulerNumber
| LibFn::Stirling1
| LibFn::Stirling2
| LibFn::PartitionCount
| LibFn::LambertW => Err(unevaluable(format!("cannot evaluate function '{f}'"))),
}
}
fn eval_bessel(
arb: fn(
&BigFloat,
&BigFloat,
usize,
RoundingMode,
&mut Consts,
) -> Result<BigFloat, SymplexError>,
what: &str,
args: &[ExprId],
cache: &FxHashMap<ExprId, Complex>,
prec: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<Complex, SymplexError> {
let order = get_cached(cache, args[0])?;
let arg = get_cached(cache, args[1])?;
if !order.1.is_zero() || !arg.1.is_zero() {
return Err(SymplexError::Unevaluable {
reason: "Bessel of complex argument not yet supported in evalf".into(),
});
}
tracing::debug!(prec, what, "evalf: Bessel function via series/asymptotic");
let result = arb(&order.0, &arg.0, prec, rm, cc)?;
Ok((result, BigFloat::new(prec)))
}
fn eval_orthopoly(
kind: OrthoPoly,
f: LibFn,
args: &[ExprId],
cache: &FxHashMap<ExprId, Complex>,
prec: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<Complex, SymplexError> {
let n_val = get_cached(cache, args[0])?;
let x_val = get_cached(cache, args[1])?;
if !n_val.1.is_zero() || !x_val.1.is_zero() {
return Err(SymplexError::Unevaluable {
reason: "orthogonal polynomial of complex arguments not supported".into(),
});
}
let n_f = bigfloat_to_f64(&n_val.0, rm, cc)?;
let n_round = n_f.round();
if (n_f - n_round).abs() > 1e-12 || !(0.0..=1.0e7).contains(&n_round) {
return Err(SymplexError::Unevaluable {
reason: format!("{f}: degree must be a non-negative integer"),
});
}
tracing::debug!(
prec,
n = n_round,
"evalf: orthogonal polynomial via recurrence"
);
let result = arb_orthopoly(kind, n_round as u64, &x_val.0, prec, rm);
Ok((result, BigFloat::new(prec)))
}
fn eval_special_09(
f: LibFn,
args: &[ExprId],
cache: &FxHashMap<ExprId, Complex>,
prec: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<Complex, SymplexError> {
let name = f.name();
let v = real_args(name, args, cache)?;
let real = |r: BigFloat| Ok((r, BigFloat::new(prec)));
debug!(prec, name, "evalf: special function (0.9)");
match f {
LibFn::Erfi => real(arb_erfi(&v[0], prec, rm, cc)?),
LibFn::ErfInv => real(arb_erfinv(&v[0], prec, rm, cc)?),
LibFn::ErfcInv => real(arb_erfcinv(&v[0], prec, rm, cc)?),
LibFn::ExpInt => real(arb_expint(&v[0], &v[1], prec, rm, cc)?),
LibFn::Shi => real(arb_shi_chi(&v[0], true, prec, rm, cc)?),
LibFn::Chi => {
if v[0].is_zero() {
return Err(unevaluable("Chi(0) is -∞"));
}
let chi = arb_shi_chi(&v[0].abs(), false, prec, rm, cc)?;
if v[0].is_negative() {
Ok((chi, cc.pi(prec, rm).clone()))
} else {
real(chi)
}
}
LibFn::FresnelS => real(arb_fresnel(&v[0], true, prec, rm, cc)?),
LibFn::FresnelC => real(arb_fresnel(&v[0], false, prec, rm, cc)?),
LibFn::LowerGamma => real(arb_lowergamma(&v[0], &v[1], prec, rm, cc)?),
LibFn::UpperGamma => real(arb_uppergamma(&v[0], &v[1], prec, rm, cc)?),
LibFn::PolyLog => real(arb_polylog(&v[0], &v[1], prec, rm, cc)?),
LibFn::DirichletEta => real(arb_dirichlet_eta(&v[0], prec, rm, cc)?),
LibFn::AiryAi => real(arb_airy(&v[0], AiryKind::Ai, prec, rm, cc)?),
LibFn::AiryBi => real(arb_airy(&v[0], AiryKind::Bi, prec, rm, cc)?),
LibFn::AiryAiPrime => real(arb_airy(&v[0], AiryKind::AiPrime, prec, rm, cc)?),
LibFn::AiryBiPrime => real(arb_airy(&v[0], AiryKind::BiPrime, prec, rm, cc)?),
LibFn::EllipticK => real(arb_elliptic_k(&v[0], prec, rm, cc)?),
LibFn::EllipticE => real(arb_elliptic_e(&v[0], prec, rm, cc)?),
LibFn::EllipticF => real(arb_elliptic_f(&v[0], &v[1], prec, rm, cc)?),
LibFn::EllipticPi => real(arb_elliptic_pi(&v[0], &v[1], prec, rm, cc)?),
LibFn::Gegenbauer => real(arb_gegenbauer(&v[0], &v[1], &v[2], prec, rm, cc)?),
LibFn::Jacobi => real(arb_jacobi(&v[0], &v[1], &v[2], &v[3], prec, rm, cc)?),
LibFn::AssocLegendre => real(arb_assoc_legendre(&v[0], &v[1], &v[2], prec, rm, cc)?),
LibFn::AssocLaguerre => real(arb_assoc_laguerre(&v[0], &v[1], &v[2], prec, rm, cc)?),
LibFn::BetaInc | LibFn::BetaIncRegularized => real(arb_betainc(
&v[0],
&v[1],
&v[2],
&v[3],
f == LibFn::BetaIncRegularized,
prec,
rm,
cc,
)?),
LibFn::Factorial2
| LibFn::Subfactorial
| LibFn::RisingFactorial
| LibFn::FallingFactorial
| LibFn::Fibonacci
| LibFn::Lucas
| LibFn::Bernoulli
| LibFn::Harmonic
| LibFn::Catalan
| LibFn::Bell
| LibFn::EulerNumber
| LibFn::Stirling1
| LibFn::Stirling2
| LibFn::PartitionCount
| LibFn::LambertW
| LibFn::BesselJ
| LibFn::BesselY
| LibFn::BesselI
| LibFn::BesselK
| LibFn::Legendre
| LibFn::ChebyshevT
| LibFn::ChebyshevU
| LibFn::Hermite
| LibFn::Laguerre => Err(unevaluable(format!("cannot evaluate function '{name}'"))),
}
}
fn arb_erfi(
x: &BigFloat,
prec: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
if x.is_zero() {
return Ok(BigFloat::new(prec));
}
if x.is_nan() || x.is_inf() {
return Err(unevaluable("erfi of special float value"));
}
let wp = prec + 32;
let x_f = bigfloat_to_f64(x, rm, cc)?;
let ax = x_f.abs();
let two = BigFloat::from_i32(2, wp);
let sqrt_pi = cc.pi(wp, rm).clone().sqrt(wp, rm);
let x_sq = x.mul(x, wp, rm);
if ax * ax > (wp as f64) * std::f64::consts::LN_2 + 4.0 {
let two_x_sq = x_sq.mul(&two, wp, rm);
let mut sum = BigFloat::from_i32(1, wp);
let mut term = BigFloat::from_i32(1, wp);
let mut prev: Option<BigFloat> = None;
for k in 1..(wp * 2 + 100) {
let f = BigFloat::from_i128(2 * k as i128 - 1, wp);
term = term.mul(&f, wp, rm).div(&two_x_sq, wp, rm);
if let Some(ref p) = prev
&& bf_gt(&term, p)
{
break;
}
sum = sum.add(&term, wp, rm);
if negligible(&term, &sum, wp) {
break;
}
prev = Some(term.clone());
}
let e = x_sq.exp(wp, rm, cc);
let denom = x.abs().mul(&sqrt_pi, wp, rm);
let r = e.mul(&sum, wp, rm).div(&denom, wp, rm);
let r = if x.is_negative() { r.neg() } else { r };
return Ok(round_to(r, prec, rm));
}
let mut prod = x.clone();
let _ = prod.set_precision(wp, rm);
let mut sum = prod.clone();
let max_terms = (3.0 * ax * ax) as usize + wp + 60;
let mut converged = false;
for n in 1..=max_terms {
let n_bf = BigFloat::from_i128(n as i128, wp);
prod = prod.mul(&x_sq, wp, rm).div(&n_bf, wp, rm);
let d = BigFloat::from_i128(2 * n as i128 + 1, wp);
let term = prod.div(&d, wp, rm);
sum = sum.add(&term, wp, rm);
if negligible(&term, &sum, wp) {
converged = true;
break;
}
}
if !converged {
return Err(series_did_not_converge("erfi", max_terms));
}
let r = sum.mul(&two, wp, rm).div(&sqrt_pi, wp, rm);
Ok(round_to(r, prec, rm))
}
fn arb_erfc(
x: &BigFloat,
prec: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
let wp = prec + 32;
let x_f = bigfloat_to_f64(x, rm, cc)?;
if x_f <= 0.5 {
let erf = arb_erf(x, prec, rm, cc)?;
let one = BigFloat::from_i32(1, wp);
return Ok(round_to(one.sub(&erf, wp, rm), prec, rm));
}
let x_sq_bits = (x_f * x_f * std::f64::consts::LOG2_E).ceil() as usize;
let mut xw = x.clone();
if x_sq_bits >= wp + 8 {
let _ = xw.set_precision(wp, rm);
return Ok(round_to(erfc_asymptotic(&xw, wp, rm, cc), prec, rm));
}
let wp2 = wp + x_sq_bits + 8;
let _ = xw.set_precision(wp2, rm);
let erf = arb_erf(&xw, wp2, rm, cc)?;
let one = BigFloat::from_i32(1, wp2);
Ok(round_to(one.sub(&erf, wp2, rm), prec, rm))
}
fn erf_inverse_halley(
y: &BigFloat,
x0: f64,
solve_erfc: bool,
prec: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
let wp = prec + 32;
let mut yw = y.clone();
let _ = yw.set_precision(wp, rm);
let half_sqrt_pi = cc
.pi(wp, rm)
.clone()
.sqrt(wp, rm)
.div(&BigFloat::from_i32(2, wp), wp, rm);
let one = BigFloat::from_i32(1, wp);
let mut x = BigFloat::from_f64(x0, wp);
for _ in 0..200 {
let f = if solve_erfc {
arb_erfc(&x, wp, rm, cc)?.sub(&yw, wp, rm)
} else {
arb_erf(&x, wp, rm, cc)?.sub(&yw, wp, rm)
};
let x_sq = x.mul(&x, wp, rm);
let g = f
.mul(&half_sqrt_pi, wp, rm)
.mul(&x_sq.exp(wp, rm, cc), wp, rm);
let g = if solve_erfc { g.neg() } else { g };
let denom = one.add(&x.mul(&g, wp, rm), wp, rm);
let delta = if denom.is_zero() {
g.clone()
} else {
g.div(&denom, wp, rm)
};
x = x.sub(&delta, wp, rm);
if negligible(&delta, &x, wp) {
break;
}
}
Ok(round_to(x, prec, rm))
}
fn arb_erfinv(
y: &BigFloat,
prec: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
if y.is_zero() {
return Ok(BigFloat::new(prec));
}
let wp = prec + 32;
let one = BigFloat::from_i32(1, wp);
let ay = y.abs();
let one_minus_ay = one.sub(&ay, wp, rm);
if one_minus_ay.is_zero() || one_minus_ay.is_negative() {
return Err(unevaluable("erfinv: argument must satisfy |y| < 1"));
}
let y_f = bigfloat_to_f64(y, rm, cc)?;
let r = if y_f.abs() > 0.5 {
arb_erfcinv_small(&one_minus_ay, prec, rm, cc)?
} else {
let a = 0.147_f64;
let l = (1.0 - y_f * y_f).ln();
let t = 2.0 / (std::f64::consts::PI * a) + l / 2.0;
let mut x0 = ((t * t - l / a).sqrt() - t).sqrt().copysign(y_f);
for _ in 0..3 {
let f = erf_f64(x0) - y_f;
x0 -= f * (std::f64::consts::PI.sqrt() / 2.0) * (x0 * x0).exp();
}
erf_inverse_halley(y, x0, false, prec, rm, cc)?
};
Ok(if y.is_negative() && !r.is_negative() && y_f.abs() > 0.5 {
r.neg()
} else {
r
})
}
fn arb_erfcinv_small(
u: &BigFloat,
prec: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
let wp = prec + 32;
let l = u.ln(wp, rm, cc).neg();
let l_f = bigfloat_to_f64(&l, rm, cc)?;
let mut x0 = l_f.max(0.5).sqrt();
for _ in 0..4 {
let v = l_f - (x0 * std::f64::consts::PI.sqrt()).ln();
x0 = v.max(0.25).sqrt();
}
let u_f = bigfloat_to_f64(u, rm, cc)?;
if u_f > 1e-300 {
for _ in 0..3 {
let f = (1.0 - erf_f64(x0)) - u_f;
x0 += f * (std::f64::consts::PI.sqrt() / 2.0) * (x0 * x0).exp();
}
}
erf_inverse_halley(u, x0, true, prec, rm, cc)
}
fn arb_erfcinv(
u: &BigFloat,
prec: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
let wp = prec + 32;
let one = BigFloat::from_i32(1, wp);
let two = BigFloat::from_i32(2, wp);
if u.is_zero() || u.is_negative() || !bf_lt(u, &two) {
return Err(unevaluable("erfcinv: argument must satisfy 0 < u < 2"));
}
if bf_gt(u, &one) {
let v = two.sub(u, wp, rm);
return Ok(arb_erfcinv(&v, prec, rm, cc)?.neg());
}
let u_f = bigfloat_to_f64(u, rm, cc)?;
if u_f >= 0.5 {
let y = one.sub(u, wp, rm);
return arb_erfinv(&y, prec, rm, cc);
}
arb_erfcinv_small(u, prec, rm, cc)
}
fn uppergamma_cf(
s: &BigFloat,
x: &BigFloat,
wp: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
let one = BigFloat::from_i32(1, wp);
let two = BigFloat::from_i32(2, wp);
let mut tiny = BigFloat::from_i32(1, wp);
tiny.set_exponent(-(4 * wp as i32));
let mut b = x.add(&one, wp, rm).sub(s, wp, rm);
let mut c = one.div(&tiny, wp, rm);
let mut d = if b.is_zero() {
one.div(&tiny, wp, rm)
} else {
one.div(&b, wp, rm)
};
let mut h = d.clone();
let max_iter = 40 * wp + 2000;
let mut converged = false;
for i in 1..=max_iter {
let i_bf = BigFloat::from_i128(i as i128, wp);
let an = i_bf.mul(&i_bf.sub(s, wp, rm), wp, rm).neg();
b = b.add(&two, wp, rm);
d = an.mul(&d, wp, rm).add(&b, wp, rm);
if d.is_zero() {
d = tiny.clone();
}
c = b.add(&an.div(&c, wp, rm), wp, rm);
if c.is_zero() {
c = tiny.clone();
}
d = one.div(&d, wp, rm);
let del = d.mul(&c, wp, rm);
h = h.mul(&del, wp, rm);
let dev = del.sub(&one, wp, rm);
if negligible(&dev, &one, wp) {
converged = true;
break;
}
}
if !converged {
return Err(SymplexError::ComputationFailed {
operation: "evalf",
reason: "continued fraction for the incomplete gamma function did not converge".into(),
});
}
let e_neg_x = x.neg().exp(wp, rm, cc);
let x_pow_s = bf_pow(x, s, wp, rm, cc);
Ok(e_neg_x.mul(&x_pow_s, wp, rm).mul(&h, wp, rm))
}
fn lowergamma_series(
s: &BigFloat,
x: &BigFloat,
wp: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
if x.is_zero() {
return Ok(BigFloat::new(wp));
}
let one = BigFloat::from_i32(1, wp);
let x_f = bigfloat_to_f64(x, rm, cc)?;
let mut term = one.div(s, wp, rm);
let mut sum = term.clone();
let max_terms = (2.0 * x_f) as usize + wp + 50;
let mut converged = false;
for k in 1..=max_terms {
let s_plus_k = s.add(&BigFloat::from_i128(k as i128, wp), wp, rm);
term = term.mul(x, wp, rm).div(&s_plus_k, wp, rm);
sum = sum.add(&term, wp, rm);
if negligible(&term, &sum, wp) {
converged = true;
break;
}
}
if !converged {
return Err(series_did_not_converge("lowergamma", max_terms));
}
let e_neg_x = x.neg().exp(wp, rm, cc);
let x_pow_s = bf_pow(x, s, wp, rm, cc);
Ok(x_pow_s.mul(&e_neg_x, wp, rm).mul(&sum, wp, rm))
}
fn uppergamma_wp(
s: &BigFloat,
x: &BigFloat,
wp: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
let s_f = bigfloat_to_f64(s, rm, cc)?;
let x_f = bigfloat_to_f64(x, rm, cc)?;
if x_f >= (s_f + 1.0).max(1.0) {
return uppergamma_cf(s, x, wp, rm, cc);
}
if s_f > 0.0 {
let wp2 = wp + 32 + cancellation_guard_bits(x_f);
let mut sw = s.clone();
let _ = sw.set_precision(wp2, rm);
let mut xw = x.clone();
let _ = xw.set_precision(wp2, rm);
let g = arb_gamma_real(&sw, wp2, rm, cc)?;
let lower = lowergamma_series(&sw, &xw, wp2, rm, cc)?;
return Ok(round_to(g.sub(&lower, wp2, rm), wp, rm));
}
let s_int = bf_as_int(s, rm, cc)?;
let steps = match s_int {
Some(k) => (-k) as usize,
None => (-s_f).ceil() as usize,
};
let wp2 = wp + 16 * steps + 32;
let mut sw = s.clone();
let _ = sw.set_precision(wp2, rm);
let mut xw = x.clone();
let _ = xw.set_precision(wp2, rm);
let e_neg_x = xw.neg().exp(wp2, rm, cc);
let (mut cur_s, mut val) = match s_int {
Some(_) => (BigFloat::new(wp2), arb_ei(&xw.neg(), wp2, rm, cc)?.neg()),
None => {
let s0 = sw.add(&BigFloat::from_i128(steps as i128, wp2), wp2, rm);
let g = arb_gamma_real(&s0, wp2, rm, cc)?;
let lower = lowergamma_series(&s0, &xw, wp2, rm, cc)?;
(s0, g.sub(&lower, wp2, rm))
}
};
let one = BigFloat::from_i32(1, wp2);
for _ in 0..steps {
let next_s = cur_s.sub(&one, wp2, rm);
let x_pow = bf_pow(&xw, &next_s, wp2, rm, cc);
let t = x_pow.mul(&e_neg_x, wp2, rm);
val = val.sub(&t, wp2, rm).div(&next_s, wp2, rm);
cur_s = next_s;
}
Ok(round_to(val, wp, rm))
}
fn arb_uppergamma(
s: &BigFloat,
x: &BigFloat,
prec: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
if x.is_nan() || s.is_nan() || x.is_inf_neg() || s.is_inf() {
return Err(unevaluable("uppergamma of special float value"));
}
if x.is_inf_pos() {
return Ok(BigFloat::new(prec));
}
if x.is_negative() {
return Err(unevaluable(
"uppergamma of negative argument not yet supported in evalf",
));
}
let wp = prec + 32;
if x.is_zero() {
if !s.is_positive() {
return Err(unevaluable("uppergamma(s, 0) diverges for s ≤ 0"));
}
return Ok(round_to(arb_gamma_real(s, wp, rm, cc)?, prec, rm));
}
let r = uppergamma_wp(s, x, wp, rm, cc)?;
Ok(round_to(r, prec, rm))
}
fn arb_lowergamma(
s: &BigFloat,
x: &BigFloat,
prec: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
if x.is_nan() || s.is_nan() || s.is_inf() || x.is_inf_neg() {
return Err(unevaluable("lowergamma of special float value"));
}
if !s.is_positive() {
return Err(unevaluable(
"lowergamma requires s > 0 for numerical evaluation",
));
}
if x.is_negative() {
return Err(unevaluable(
"lowergamma of negative argument not yet supported in evalf",
));
}
let wp = prec + 32;
if x.is_inf_pos() {
return Ok(round_to(arb_gamma_real(s, wp, rm, cc)?, prec, rm));
}
let s_f = bigfloat_to_f64(s, rm, cc)?;
let x_f = bigfloat_to_f64(x, rm, cc)?;
if x_f < s_f + 1.0 {
let r = lowergamma_series(s, x, wp, rm, cc)?;
return Ok(round_to(r, prec, rm));
}
let g = arb_gamma_real(s, wp, rm, cc)?;
let upper = uppergamma_cf(s, x, wp, rm, cc)?;
Ok(round_to(g.sub(&upper, wp, rm), prec, rm))
}
fn arb_expint(
n: &BigFloat,
x: &BigFloat,
prec: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
if x.is_nan() || n.is_nan() || n.is_inf() {
return Err(unevaluable("expint of special float value"));
}
if x.is_inf_pos() {
return Ok(BigFloat::new(prec));
}
let wp = prec + 32;
let one = BigFloat::from_i32(1, wp);
if x.is_zero() {
let n_minus_1 = n.sub(&one, wp, rm);
if !n_minus_1.is_positive() {
return Err(unevaluable("expint(n, 0) diverges for n ≤ 1"));
}
return Ok(round_to(one.div(&n_minus_1, wp, rm), prec, rm));
}
if x.is_negative() {
return Err(unevaluable(
"expint of negative argument is complex; not yet supported in evalf",
));
}
let s = one.sub(n, wp, rm);
let g = uppergamma_wp(&s, x, wp, rm, cc)?;
let x_pow = bf_pow(x, &n.sub(&one, wp, rm), wp, rm, cc);
Ok(round_to(x_pow.mul(&g, wp, rm), prec, rm))
}
fn arb_shi_chi(
x: &BigFloat,
want_shi: bool,
prec: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
if x.is_nan() || x.is_inf() {
return Err(unevaluable("Shi/Chi of special float value"));
}
if x.is_zero() {
return if want_shi {
Ok(BigFloat::new(prec))
} else {
Err(unevaluable("Chi(0) is -∞"))
};
}
let wp = prec + 32;
let ax = x.abs();
let x_f = bigfloat_to_f64(&ax, rm, cc)?;
let r = if x_f > 40.0 {
let ei_pos = arb_ei(&ax, wp, rm, cc)?;
let ei_neg = arb_ei(&ax.neg(), wp, rm, cc)?;
let two = BigFloat::from_i32(2, wp);
if want_shi {
ei_pos.sub(&ei_neg, wp, rm).div(&two, wp, rm)
} else {
ei_pos.add(&ei_neg, wp, rm).div(&two, wp, rm)
}
} else {
let mut xw = ax.clone();
let _ = xw.set_precision(wp, rm);
let x2 = xw.mul(&xw, wp, rm);
let max_terms = (x_f * 1.5) as usize + wp + 40;
let name = if want_shi { "Shi" } else { "Chi" };
let mut converged = false;
let series = if want_shi {
let mut t = xw.clone();
let mut sum = BigFloat::new(wp);
for k in 0..max_terms {
if k > 0 {
let d = BigFloat::from_i128((2 * k as i128) * (2 * k as i128 + 1), wp);
t = t.mul(&x2, wp, rm).div(&d, wp, rm);
}
let term = t.div(&BigFloat::from_i128(2 * k as i128 + 1, wp), wp, rm);
sum = sum.add(&term, wp, rm);
if negligible(&term, &sum, wp) {
converged = true;
break;
}
}
sum
} else {
let gamma = arb_euler_gamma(wp, rm, cc)?;
let ln_x = xw.ln(wp, rm, cc);
let mut u = BigFloat::from_i32(1, wp); let mut sum = BigFloat::new(wp);
for k in 1..max_terms {
let d = BigFloat::from_i128((2 * k as i128 - 1) * (2 * k as i128), wp);
u = u.mul(&x2, wp, rm).div(&d, wp, rm);
let term = u.div(&BigFloat::from_i128(2 * k as i128, wp), wp, rm);
sum = sum.add(&term, wp, rm);
if negligible(&term, &sum, wp) {
converged = true;
break;
}
}
gamma.add(&ln_x, wp, rm).add(&sum, wp, rm)
};
if !converged {
return Err(series_did_not_converge(name, max_terms));
}
series
};
let r = if want_shi && x.is_negative() {
r.neg()
} else {
r
};
Ok(round_to(r, prec, rm))
}
fn arb_fresnel(
x: &BigFloat,
want_s: bool,
prec: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
if x.is_nan() || x.is_inf() {
return Err(unevaluable("Fresnel integral of special float value"));
}
if x.is_zero() {
return Ok(BigFloat::new(prec));
}
let base_wp = prec + 32;
let ax = x.abs();
let x_f = bigfloat_to_f64(&ax, rm, cc)?;
let arg_f = std::f64::consts::FRAC_PI_2 * x_f * x_f;
let r = if arg_f > (base_wp as f64) * std::f64::consts::LN_2 {
let wp = base_wp;
let mut xw = ax.clone();
let _ = xw.set_precision(wp, rm);
let one = BigFloat::from_i32(1, wp);
let pi = cc.pi(wp, rm).clone();
let pi_x2 = pi.mul(&xw, wp, rm).mul(&xw, wp, rm);
let inv = one.div(&pi_x2, wp, rm); let inv2 = inv.mul(&inv, wp, rm);
let mut tf = one.clone();
let mut tg = inv.clone();
let mut f_sum = BigFloat::new(wp);
let mut g_sum = BigFloat::new(wp);
let mut prev: Option<BigFloat> = None;
for m in 0..(wp * 2 + 100) {
if m > 0 {
let a = BigFloat::from_i128((4 * m as i128 - 3) * (4 * m as i128 - 1), wp);
let b = BigFloat::from_i128((4 * m as i128 - 1) * (4 * m as i128 + 1), wp);
tf = tf.mul(&a, wp, rm).mul(&inv2, wp, rm);
tg = tg.mul(&b, wp, rm).mul(&inv2, wp, rm);
}
if let Some(ref p) = prev
&& bf_gt(&tf, p)
{
break;
}
if m % 2 == 0 {
f_sum = f_sum.add(&tf, wp, rm);
g_sum = g_sum.add(&tg, wp, rm);
} else {
f_sum = f_sum.sub(&tf, wp, rm);
g_sum = g_sum.sub(&tg, wp, rm);
}
if negligible(&tf, &one, wp) {
break;
}
prev = Some(tf.clone());
}
let pi_x = pi.mul(&xw, wp, rm);
let f = f_sum.div(&pi_x, wp, rm);
let g = g_sum.div(&pi_x, wp, rm);
let two = BigFloat::from_i32(2, wp);
let theta = pi_x2.div(&two, wp, rm);
let s = theta.sin(wp, rm, cc);
let c = theta.cos(wp, rm, cc);
let half = one.div(&two, wp, rm);
if want_s {
half.sub(&f.mul(&c, wp, rm), wp, rm)
.sub(&g.mul(&s, wp, rm), wp, rm)
} else {
half.add(&f.mul(&s, wp, rm), wp, rm)
.sub(&g.mul(&c, wp, rm), wp, rm)
}
} else {
let wp = base_wp + (arg_f * std::f64::consts::LOG2_E).ceil() as usize + 16;
let mut xw = ax.clone();
let _ = xw.set_precision(wp, rm);
let pi = cc.pi(wp, rm).clone();
let two = BigFloat::from_i32(2, wp);
let half_pi = pi.div(&two, wp, rm);
let x2 = xw.mul(&xw, wp, rm);
let u = half_pi.mul(&x2, wp, rm);
let u2 = u.mul(&u, wp, rm);
let mut t = if want_s {
u.mul(&xw, wp, rm)
} else {
xw.clone()
};
let mut sum = BigFloat::new(wp);
let max_terms = (arg_f * 2.0) as usize + wp + 40;
let mut converged = false;
for k in 0..max_terms {
if k > 0 {
let d = if want_s {
BigFloat::from_i128((2 * k as i128) * (2 * k as i128 + 1), wp)
} else {
BigFloat::from_i128((2 * k as i128 - 1) * (2 * k as i128), wp)
};
t = t.mul(&u2, wp, rm).div(&d, wp, rm);
}
let div = BigFloat::from_i128(4 * k as i128 + if want_s { 3 } else { 1 }, wp);
let term = t.div(&div, wp, rm);
if k % 2 == 0 {
sum = sum.add(&term, wp, rm);
} else {
sum = sum.sub(&term, wp, rm);
}
if negligible(&term, &sum, wp) {
converged = true;
break;
}
}
if !converged {
let name = if want_s { "fresnels" } else { "fresnelc" };
return Err(series_did_not_converge(name, max_terms));
}
sum
};
let r = if x.is_negative() { r.neg() } else { r };
Ok(round_to(r, prec, rm))
}
fn beta_wp(
a: &BigFloat,
b: &BigFloat,
wp: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
let ga = arb_gamma_real(a, wp, rm, cc)?;
let gb = arb_gamma_real(b, wp, rm, cc)?;
let a_plus_b = a.add(b, wp, rm);
let gab = arb_gamma_real(&a_plus_b, wp, rm, cc)?;
if gab.is_zero() {
return Err(unevaluable("betainc: Gamma(a + b) is zero"));
}
Ok(ga.mul(&gb, wp, rm).div(&gab, wp, rm))
}
fn betainc_cf(
a: &BigFloat,
b: &BigFloat,
x: &BigFloat,
wp: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
let one = BigFloat::from_i32(1, wp);
let mut tiny = BigFloat::from_i32(1, wp);
tiny.set_exponent(-(4 * wp as i32));
let qab = a.add(b, wp, rm);
let qap = a.add(&one, wp, rm);
let qam = a.sub(&one, wp, rm);
let mut c = one.clone();
let mut d = one.sub(&qab.mul(x, wp, rm).div(&qap, wp, rm), wp, rm);
if d.is_zero() {
d = tiny.clone();
}
d = one.div(&d, wp, rm);
let mut h = d.clone();
let ab_f = bigfloat_to_f64(&qab, rm, cc)?;
let max_iter = 40 * wp + 2000 + (10.0 * ab_f.abs().sqrt()).ceil() as usize;
let mut converged = false;
let step = |aa: &BigFloat, c: &mut BigFloat, d: &mut BigFloat| -> BigFloat {
*d = one.add(&aa.mul(d, wp, rm), wp, rm);
if d.is_zero() {
*d = tiny.clone();
}
*c = one.add(&aa.div(c, wp, rm), wp, rm);
if c.is_zero() {
*c = tiny.clone();
}
*d = one.div(d, wp, rm);
d.mul(c, wp, rm)
};
for m in 1..=max_iter {
let m_bf = BigFloat::from_i128(m as i128, wp);
let m2_bf = BigFloat::from_i128(2 * m as i128, wp);
let num = m_bf.mul(&b.sub(&m_bf, wp, rm), wp, rm).mul(x, wp, rm);
let den = qam.add(&m2_bf, wp, rm).mul(&a.add(&m2_bf, wp, rm), wp, rm);
let aa = num.div(&den, wp, rm);
let del = step(&aa, &mut c, &mut d);
h = h.mul(&del, wp, rm);
let num = a
.add(&m_bf, wp, rm)
.mul(&qab.add(&m_bf, wp, rm), wp, rm)
.mul(x, wp, rm)
.neg();
let den = a.add(&m2_bf, wp, rm).mul(&qap.add(&m2_bf, wp, rm), wp, rm);
let aa = num.div(&den, wp, rm);
let del = step(&aa, &mut c, &mut d);
h = h.mul(&del, wp, rm);
let dev = del.sub(&one, wp, rm);
if negligible(&dev, &one, wp) {
converged = true;
break;
}
}
if !converged {
return Err(SymplexError::ComputationFailed {
operation: "evalf",
reason: format!(
"continued fraction for the incomplete beta function did not converge in \
{max_iter} steps"
),
});
}
Ok(h)
}
fn betainc_reg_direct(
a: &BigFloat,
b: &BigFloat,
x: &BigFloat,
one_minus_x: &BigFloat,
wp: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
let h = betainc_cf(a, b, x, wp, rm, cc)?;
let xa = bf_pow(x, a, wp, rm, cc);
let yb = bf_pow(one_minus_x, b, wp, rm, cc);
let beta = beta_wp(a, b, wp, rm, cc)?;
let denom = a.mul(&beta, wp, rm);
Ok(xa.mul(&yb, wp, rm).div(&denom, wp, rm).mul(&h, wp, rm))
}
fn betainc_switch_point(a: &BigFloat, b: &BigFloat, wp: usize, rm: RoundingMode) -> BigFloat {
let one = BigFloat::from_i32(1, wp);
let two = BigFloat::from_i32(2, wp);
let a_plus_1 = a.add(&one, wp, rm);
let a_plus_b_plus_2 = a.add(b, wp, rm).add(&two, wp, rm);
a_plus_1.div(&a_plus_b_plus_2, wp, rm)
}
fn betainc_reg_wp(
a: &BigFloat,
b: &BigFloat,
x: &BigFloat,
wp: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
if x.is_zero() {
return Ok(BigFloat::new(wp));
}
let one = BigFloat::from_i32(1, wp);
let one_minus_x = one.sub(x, wp, rm);
if one_minus_x.is_zero() {
return Ok(one);
}
if bf_gt(x, &betainc_switch_point(a, b, wp, rm)) {
let v = betainc_reg_direct(b, a, &one_minus_x, x, wp, rm, cc)?;
return Ok(one.sub(&v, wp, rm));
}
betainc_reg_direct(a, b, x, &one_minus_x, wp, rm, cc)
}
fn betainc_reg_diff_wp(
a: &BigFloat,
b: &BigFloat,
lo: &BigFloat,
hi: &BigFloat,
wp: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<(BigFloat, usize), SymplexError> {
let one = BigFloat::from_i32(1, wp);
let (p, q) = if bf_gt(lo, &betainc_switch_point(a, b, wp, rm)) {
let y_lo = one.sub(hi, wp, rm); let y_hi = one.sub(lo, wp, rm); let p = betainc_reg_direct(b, a, &y_hi, lo, wp, rm, cc)?;
let q = if y_lo.is_zero() {
BigFloat::new(wp)
} else {
betainc_reg_direct(b, a, &y_lo, hi, wp, rm, cc)?
};
(p, q)
} else {
let p = betainc_reg_wp(a, b, hi, wp, rm, cc)?;
let q = betainc_reg_wp(a, b, lo, wp, rm, cc)?;
(p, q)
};
let diff = p.sub(&q, wp, rm);
let lost = match (p.exponent(), q.exponent(), diff.exponent()) {
_ if q.is_zero() => 0,
_ if diff.is_zero() => wp,
(Some(ep), Some(eq), Some(ed)) => (ep.max(eq) as i64 - ed as i64).max(0) as usize,
_ => 0,
};
Ok((diff, lost))
}
#[allow(clippy::too_many_arguments)]
fn arb_betainc(
a: &BigFloat,
b: &BigFloat,
x1: &BigFloat,
x2: &BigFloat,
regularized: bool,
prec: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
let name = if regularized {
"betainc_regularized"
} else {
"betainc"
};
if [a, b, x1, x2].iter().any(|v| v.is_nan() || v.is_inf()) {
return Err(unevaluable(format!("{name} of special float value")));
}
if !a.is_positive() || !b.is_positive() {
return Err(unevaluable(format!(
"{name} requires a, b > 0 for numerical evaluation"
)));
}
let one = BigFloat::from_i32(1, prec);
if [x1, x2].iter().any(|x| x.is_negative() || bf_gt(x, &one)) {
return Err(unevaluable(format!(
"{name}: the limits x1, x2 must lie in [0, 1] for numerical evaluation"
)));
}
let (lo, hi, negate) = match x1.cmp(x2) {
Some(c) if c < 0 => (x1, x2, false),
Some(c) if c > 0 => (x2, x1, true),
_ => return Ok(BigFloat::new(prec)),
};
let mut wp = prec + 32;
let mut diff;
let mut attempts = 0;
loop {
let (d, lost) = betainc_reg_diff_wp(a, b, lo, hi, wp, rm, cc)?;
diff = d;
attempts += 1;
if lost + prec + 16 <= wp || attempts >= 3 {
break;
}
wp = prec + 32 + lost;
}
if !regularized {
let beta = beta_wp(a, b, wp, rm, cc)?;
diff = diff.mul(&beta, wp, rm);
}
if negate {
diff = diff.neg();
}
Ok(round_to(diff, prec, rm))
}
fn zeta_nonpositive_int(m: i64, wp: usize, rm: RoundingMode) -> BigFloat {
if m == 0 {
return BigFloat::from_f64(-0.5, wp);
}
let k = (-m) as usize;
let b = crate::base::bernoulli::bernoulli(k + 1);
if b.is_zero() {
return BigFloat::new(wp);
}
let v = -b / Ratio::from_integer(BigInt::from(k as i64 + 1));
ratio_to_bigfloat(&v, wp, rm)
}
fn int_pow_s(
k: usize,
s: &BigFloat,
s_int: Option<i64>,
wp: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> BigFloat {
let kb = BigFloat::from_i128(k as i128, wp);
match s_int {
Some(n) if n >= 0 => kb.powi(n as usize, wp, rm),
Some(n) => BigFloat::from_i32(1, wp).div(&kb.powi((-n) as usize, wp, rm), wp, rm),
None => bf_pow(&kb, s, wp, rm, cc),
}
}
fn polylog_series_terms(s_f: f64, z_f: f64, wp: usize) -> usize {
if !(z_f > 0.0 && z_f < 1.0) {
return wp * 2 + 60;
}
let a = -z_f.ln(); let target = (wp as f64) * std::f64::consts::LN_2 + a;
let mut k = target / a;
let m = (-s_f).max(0.0);
if m > 0.0 {
k = k.max(m / a);
for _ in 0..200 {
if k * a - m * k.ln() >= target {
break;
}
k *= 2.0;
}
}
(k.ceil() as usize).saturating_add(20)
}
fn polylog_series(
s: &BigFloat,
s_int: Option<i64>,
z: &BigFloat,
wp: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
let (sum, lost) = polylog_series_raw(s, s_int, z, wp, rm, cc)?;
if lost <= 4 {
return Ok(sum);
}
let wp2 = wp + lost + 8;
let mut sw = s.clone();
let _ = sw.set_precision(wp2, rm);
let mut zw = z.clone();
let _ = zw.set_precision(wp2, rm);
let (sum, _) = polylog_series_raw(&sw, s_int, &zw, wp2, rm, cc)?;
Ok(round_to(sum, wp, rm))
}
fn polylog_series_raw(
s: &BigFloat,
s_int: Option<i64>,
z: &BigFloat,
wp: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<(BigFloat, usize), SymplexError> {
let z_f = bigfloat_to_f64(z, rm, cc).unwrap_or(0.5).abs();
let s_f = bigfloat_to_f64(s, rm, cc).unwrap_or(0.0);
let max_terms = polylog_series_terms(s_f, z_f, wp);
let mut zk = z.clone();
let mut sum = BigFloat::new(wp);
let mut max_exp: Option<i64> = None;
let mut converged = false;
for k in 1..=max_terms {
if k > 1 {
zk = zk.mul(z, wp, rm);
}
let term = zk.div(&int_pow_s(k, s, s_int, wp, rm, cc), wp, rm);
if let Some(e) = term.exponent() {
max_exp = Some(max_exp.map_or(e as i64, |m| m.max(e as i64)));
}
sum = sum.add(&term, wp, rm);
if negligible(&term, &sum, wp) {
converged = true;
break;
}
}
if !converged {
return Err(series_did_not_converge("polylog", max_terms));
}
let lost = match (max_exp, sum.exponent()) {
(Some(m), Some(se)) if m > se as i64 => (m - se as i64) as usize,
_ => 0,
};
Ok((sum, lost))
}
fn polylog_nonpositive(
n: usize,
z: &BigFloat,
wp: usize,
rm: RoundingMode,
) -> Result<BigFloat, SymplexError> {
let one = BigFloat::from_i32(1, wp);
let one_minus_z = one.sub(z, wp, rm);
if one_minus_z.is_zero() {
return Err(unevaluable("polylog(s, 1) diverges for s ≤ 1"));
}
let w = z.div(&one_minus_z, wp, rm);
let mut k_fact = BigInt::from(1);
let mut w_pow = w.clone();
let mut sum = BigFloat::new(wp);
for k in 0..=n {
if k > 0 {
k_fact *= BigInt::from(k as u64);
w_pow = w_pow.mul(&w, wp, rm);
}
let s2 = crate::domains::combinatorics::stirling2(n as u64 + 1, k as u64 + 1)
.ok_or_else(|| unevaluable("polylog: Stirling number overflow"))?;
let c = ratio_to_bigfloat(&Ratio::from_integer(&k_fact * s2), wp, rm);
sum = sum.add(&c.mul(&w_pow, wp, rm), wp, rm);
}
Ok(sum)
}
fn polylog_near_one(
s: &BigFloat,
s_int: Option<i64>,
z: &BigFloat,
wp: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
let one = BigFloat::from_i32(1, wp);
let mu = z.ln(wp, rm, cc); let neg_mu = mu.neg();
let max_terms = wp * 2 + 100;
let not_converged = || series_did_not_converge("polylog", max_terms);
let mut sum = BigFloat::new(wp);
let mut mu_pow_over_fact = one.clone(); let mut small_run = 0u8;
let mut converged = false;
match s_int {
Some(n) if n >= 1 => {
let n = n as usize;
for k in 0..max_terms {
if k > 0 {
mu_pow_over_fact = mu_pow_over_fact.mul(&mu, wp, rm).div(
&BigFloat::from_i128(k as i128, wp),
wp,
rm,
);
}
if k + 1 == n {
let mut h = BigFloat::new(wp);
for j in 1..n {
h = h.add(
&one.div(&BigFloat::from_i128(j as i128, wp), wp, rm),
wp,
rm,
);
}
let ln_neg_mu = neg_mu.ln(wp, rm, cc);
let t = mu_pow_over_fact.mul(&h.sub(&ln_neg_mu, wp, rm), wp, rm);
sum = sum.add(&t, wp, rm);
continue;
}
let arg = n as i64 - k as i64;
let zeta = if arg <= 0 {
zeta_nonpositive_int(arg, wp, rm)
} else {
arb_zeta(&BigFloat::from_i128(arg as i128, wp), wp, rm, cc)?
};
let term = zeta.mul(&mu_pow_over_fact, wp, rm);
sum = sum.add(&term, wp, rm);
if k >= n {
if negligible(&term, &sum, wp) {
small_run += 1;
if small_run >= 2 {
converged = true;
break;
}
} else {
small_run = 0;
}
}
}
if !converged {
return Err(not_converged());
}
Ok(sum)
}
_ => {
let one_minus_s = one.sub(s, wp, rm);
let g = arb_gamma_real(&one_minus_s, wp, rm, cc)?;
let s_minus_1 = s.sub(&one, wp, rm);
let lead = g.mul(&bf_pow(&neg_mu, &s_minus_1, wp, rm, cc), wp, rm);
for k in 0..max_terms {
if k > 0 {
mu_pow_over_fact = mu_pow_over_fact.mul(&mu, wp, rm).div(
&BigFloat::from_i128(k as i128, wp),
wp,
rm,
);
}
let arg = s.sub(&BigFloat::from_i128(k as i128, wp), wp, rm);
let zeta = arb_zeta(&arg, wp, rm, cc)?;
let term = zeta.mul(&mu_pow_over_fact, wp, rm);
sum = sum.add(&term, wp, rm);
if k >= 2 {
if negligible(&term, &sum, wp) {
small_run += 1;
if small_run >= 2 {
converged = true;
break;
}
} else {
small_run = 0;
}
}
}
if !converged {
return Err(not_converged());
}
Ok(lead.add(&sum, wp, rm))
}
}
}
fn polylog_series_limit(s_int: Option<i64>, wp: usize) -> BigFloat {
if s_int.is_some() {
BigFloat::from_f64(0.5, wp)
} else {
BigFloat::from_f64(0.95, wp)
}
}
fn polylog_unit_interval(
s: &BigFloat,
s_int: Option<i64>,
z: &BigFloat,
wp: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
let limit = polylog_series_limit(s_int, wp);
let az = z.abs();
if !bf_gt(&az, &limit) {
return polylog_series(s, s_int, z, wp, rm, cc);
}
if z.is_negative() {
let one = BigFloat::from_i32(1, wp);
let x2 = az.mul(&az, wp, rm);
let a = polylog_unit_interval(s, s_int, &x2, wp, rm, cc)?;
let b = polylog_near_one(s, s_int, &az, wp, rm, cc)?;
let two = BigFloat::from_i32(2, wp);
let one_minus_s = one.sub(s, wp, rm);
let scale = bf_pow(&two, &one_minus_s, wp, rm, cc);
return Ok(scale.mul(&a, wp, rm).sub(&b, wp, rm));
}
polylog_near_one(s, s_int, z, wp, rm, cc)
}
fn arb_polylog(
s: &BigFloat,
z: &BigFloat,
prec: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
if s.is_nan() || z.is_nan() || s.is_inf() || z.is_inf() {
return Err(unevaluable("polylog of special float value"));
}
if z.is_zero() {
return Ok(BigFloat::new(prec));
}
let wp = prec + 32;
let one = BigFloat::from_i32(1, wp);
let s_int = bf_as_int(s, rm, cc)?;
if let Some(n) = s_int
&& n <= 0
{
let r = polylog_nonpositive((-n) as usize, z, wp, rm)?;
return Ok(round_to(r, prec, rm));
}
let az = z.abs();
if bf_gt(&az, &one) {
return Err(unevaluable(
"polylog outside the unit disc not yet supported in evalf",
));
}
if az.sub(&one, wp, rm).is_zero() {
if z.is_negative() {
let eta = arb_eta(s, wp, rm, cc)?;
return Ok(round_to(eta.neg(), prec, rm));
}
if !bf_gt(s, &one) {
return Err(unevaluable("polylog(s, 1) diverges for s ≤ 1"));
}
let zeta = arb_zeta(s, wp, rm, cc)?;
return Ok(round_to(zeta, prec, rm));
}
let r = polylog_unit_interval(s, s_int, z, wp, rm, cc)?;
Ok(round_to(r, prec, rm))
}
fn arb_eta(
s: &BigFloat,
wp: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
if s.is_positive() {
return arb_eta_borwein(s, wp, rm, cc);
}
let one = BigFloat::from_i32(1, wp);
let two = BigFloat::from_i32(2, wp);
let one_minus_s = one.sub(s, wp, rm);
let coeff = one.sub(&bf_pow(&two, &one_minus_s, wp, rm, cc), wp, rm);
let zeta = arb_zeta(s, wp, rm, cc)?;
Ok(coeff.mul(&zeta, wp, rm))
}
fn arb_dirichlet_eta(
s: &BigFloat,
prec: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
if s.is_nan() || s.is_inf_neg() {
return Err(unevaluable("dirichlet_eta of special float value"));
}
if s.is_inf_pos() {
return Ok(BigFloat::from_i32(1, prec));
}
let wp = prec + 32;
let mut sw = s.clone();
let _ = sw.set_precision(wp, rm);
Ok(round_to(arb_eta(&sw, wp, rm, cc)?, prec, rm))
}
#[derive(Clone, Copy, PartialEq, Eq, Debug)]
enum AiryKind {
Ai,
Bi,
AiPrime,
BiPrime,
}
fn arb_airy(
x: &BigFloat,
kind: AiryKind,
prec: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
if x.is_nan() || x.is_inf() {
return Err(unevaluable("Airy function of special float value"));
}
let base_wp = prec + 32;
let x_f = bigfloat_to_f64(x, rm, cc)?;
let ax_f = x_f.abs();
let zeta_f = 2.0 / 3.0 * ax_f.powf(1.5);
if 2.0 * zeta_f > (base_wp as f64) * std::f64::consts::LN_2 {
return Ok(round_to(
airy_asymptotic(x, kind, base_wp, rm, cc),
prec,
rm,
));
}
let wp = base_wp + (2.0 * zeta_f * std::f64::consts::LOG2_E).ceil() as usize + 16;
let mut xw = x.clone();
let _ = xw.set_precision(wp, rm);
let one = BigFloat::from_i32(1, wp);
let three = BigFloat::from_i32(3, wp);
let third = one.div(&three, wp, rm);
let two_thirds = third.mul(&BigFloat::from_i32(2, wp), wp, rm);
let g13 = arb_gamma_real(&third, wp, rm, cc)?;
let g23 = arb_gamma_real(&two_thirds, wp, rm, cc)?;
let c1 = one.div(
&bf_pow(&three, &two_thirds, wp, rm, cc).mul(&g23, wp, rm),
wp,
rm,
);
let c2 = one.div(
&bf_pow(&three, &third, wp, rm, cc).mul(&g13, wp, rm),
wp,
rm,
);
let x3 = xw.mul(&xw, wp, rm).mul(&xw, wp, rm);
let derivative = matches!(kind, AiryKind::AiPrime | AiryKind::BiPrime);
let (mut f_sum, mut g_sum) = if derivative {
(BigFloat::new(wp), one.clone()) } else {
(one.clone(), xw.clone()) };
let mut fk = one.clone(); let mut gk = xw.clone(); let max_terms = if xw.is_zero() {
0
} else {
(ax_f * 2.0) as usize + wp + 40
};
for k in 1..=max_terms {
let k3 = 3 * k as i128;
fk = fk
.mul(&x3, wp, rm)
.div(&BigFloat::from_i128((k3 - 1) * k3, wp), wp, rm);
gk = gk
.mul(&x3, wp, rm)
.div(&BigFloat::from_i128(k3 * (k3 + 1), wp), wp, rm);
let (tf, tg) = if derivative {
(
fk.mul(&BigFloat::from_i128(k3, wp), wp, rm)
.div(&xw, wp, rm),
gk.mul(&BigFloat::from_i128(k3 + 1, wp), wp, rm)
.div(&xw, wp, rm),
)
} else {
(fk.clone(), gk.clone())
};
f_sum = f_sum.add(&tf, wp, rm);
g_sum = g_sum.add(&tg, wp, rm);
if negligible(&tf, &one, wp + 4) && negligible(&tg, &one, wp + 4) {
break;
}
}
let a = c1.mul(&f_sum, wp, rm);
let b = c2.mul(&g_sum, wp, rm);
let r = match kind {
AiryKind::Ai | AiryKind::AiPrime => a.sub(&b, wp, rm),
AiryKind::Bi | AiryKind::BiPrime => three.sqrt(wp, rm).mul(&a.add(&b, wp, rm), wp, rm),
};
Ok(round_to(r, prec, rm))
}
fn airy_asymptotic(
x: &BigFloat,
kind: AiryKind,
wp: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> BigFloat {
let one = BigFloat::from_i32(1, wp);
let two = BigFloat::from_i32(2, wp);
let ax = x.abs();
let sqrt_ax = ax.sqrt(wp, rm);
let x_quarter = sqrt_ax.sqrt(wp, rm); let zeta = two
.mul(&ax, wp, rm)
.mul(&sqrt_ax, wp, rm)
.div(&BigFloat::from_i32(3, wp), wp, rm);
let inv_zeta = one.div(&zeta, wp, rm);
let sqrt_pi = cc.pi(wp, rm).clone().sqrt(wp, rm);
let derivative = matches!(kind, AiryKind::AiPrime | AiryKind::BiPrime);
let mut coeffs: Vec<BigFloat> = Vec::new();
let mut u = one.clone();
let mut zeta_pow = one.clone();
let mut prev: Option<BigFloat> = None;
for k in 0..(wp * 2 + 100) {
if k > 0 {
let k_i = k as i128;
let num = BigFloat::from_i128((6 * k_i - 5) * (6 * k_i - 3) * (6 * k_i - 1), wp);
let den = BigFloat::from_i128(216 * k_i * (2 * k_i - 1), wp);
u = u.mul(&num, wp, rm).div(&den, wp, rm);
zeta_pow = zeta_pow.mul(&inv_zeta, wp, rm);
}
let mut c = u.mul(&zeta_pow, wp, rm);
if derivative {
let k_i = k as i128;
let ratio = BigFloat::from_i128(6 * k_i + 1, wp).div(
&BigFloat::from_i128(6 * k_i - 1, wp),
wp,
rm,
);
c = c.mul(&ratio, wp, rm).neg();
}
if let Some(ref p) = prev
&& bf_gt(&c.abs(), p)
{
break;
}
let done = negligible(&c, &one, wp);
prev = Some(c.abs());
coeffs.push(c);
if done {
break;
}
}
if !x.is_negative() {
let mut sum = BigFloat::new(wp);
let alternate = matches!(kind, AiryKind::Ai | AiryKind::AiPrime);
for (k, c) in coeffs.iter().enumerate() {
if alternate && k % 2 == 1 {
sum = sum.sub(c, wp, rm);
} else {
sum = sum.add(c, wp, rm);
}
}
let amp = match kind {
AiryKind::Ai => zeta.neg().exp(wp, rm, cc).div(
&two.mul(&sqrt_pi, wp, rm).mul(&x_quarter, wp, rm),
wp,
rm,
),
AiryKind::Bi => zeta
.exp(wp, rm, cc)
.div(&sqrt_pi.mul(&x_quarter, wp, rm), wp, rm),
AiryKind::AiPrime => x_quarter
.mul(&zeta.neg().exp(wp, rm, cc), wp, rm)
.div(&two.mul(&sqrt_pi, wp, rm), wp, rm)
.neg(),
AiryKind::BiPrime => x_quarter
.mul(&zeta.exp(wp, rm, cc), wp, rm)
.div(&sqrt_pi, wp, rm),
};
return amp.mul(&sum, wp, rm);
}
let mut p = BigFloat::new(wp);
let mut q = BigFloat::new(wp);
for (k, c) in coeffs.iter().enumerate() {
let sign_neg = (k / 2) % 2 == 1;
let target = if k % 2 == 0 { &mut p } else { &mut q };
*target = if sign_neg {
target.sub(c, wp, rm)
} else {
target.add(c, wp, rm)
};
}
let pi = cc.pi(wp, rm).clone();
let theta = zeta.add(&pi.div(&BigFloat::from_i32(4, wp), wp, rm), wp, rm);
let s = theta.sin(wp, rm, cc);
let c = theta.cos(wp, rm, cc);
let sp_minus_cq = s.mul(&p, wp, rm).sub(&c.mul(&q, wp, rm), wp, rm);
let cp_plus_sq = c.mul(&p, wp, rm).add(&s.mul(&q, wp, rm), wp, rm);
match kind {
AiryKind::Ai => sp_minus_cq.div(&sqrt_pi.mul(&x_quarter, wp, rm), wp, rm),
AiryKind::Bi => cp_plus_sq.div(&sqrt_pi.mul(&x_quarter, wp, rm), wp, rm),
AiryKind::AiPrime => x_quarter
.mul(&cp_plus_sq, wp, rm)
.div(&sqrt_pi, wp, rm)
.neg(),
AiryKind::BiPrime => x_quarter.mul(&sp_minus_cq, wp, rm).div(&sqrt_pi, wp, rm),
}
}
fn all_close(vals: &[&BigFloat], bits: usize, wp: usize, rm: RoundingMode) -> bool {
let first = vals[0];
vals[1..]
.iter()
.all(|v| negligible(&v.sub(first, wp, rm), first, bits))
}
fn carlson_bits(wp: usize) -> usize {
wp / 2 + 8
}
fn carlson_rc(
x: &BigFloat,
y: &BigFloat,
wp: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> BigFloat {
let one = BigFloat::from_i32(1, wp);
let diff = y.sub(x, wp, rm);
let inv_sqrt_y = one.div(&y.sqrt(wp, rm), wp, rm);
if diff.is_zero() {
return inv_sqrt_y;
}
let w = diff.div(y, wp, rm);
let w_f = bigfloat_to_f64(&w, rm, cc).unwrap_or(1.0);
if w_f.abs() < 0.1 {
let mut sum = one.clone();
let mut c = one.clone();
let mut wk = one.clone();
for k in 1..(wp * 2 + 100) {
let k_i = k as i128;
let num = BigFloat::from_i128((2 * k_i - 1) * (2 * k_i - 1), wp);
let den = BigFloat::from_i128(2 * k_i * (2 * k_i + 1), wp);
c = c.mul(&num, wp, rm).div(&den, wp, rm);
wk = wk.mul(&w, wp, rm);
let term = c.mul(&wk, wp, rm);
sum = sum.add(&term, wp, rm);
if negligible(&term, &sum, wp) {
break;
}
}
return inv_sqrt_y.mul(&sum, wp, rm);
}
let ratio = x.div(y, wp, rm).sqrt(wp, rm); if diff.is_positive() {
ratio.acos(wp, rm, cc).div(&diff.sqrt(wp, rm), wp, rm)
} else {
ratio
.acosh(wp, rm, cc)
.div(&diff.neg().sqrt(wp, rm), wp, rm)
}
}
fn carlson_rf(
x: &BigFloat,
y: &BigFloat,
z: &BigFloat,
wp: usize,
rm: RoundingMode,
) -> Result<BigFloat, SymplexError> {
let four = BigFloat::from_i32(4, wp);
let (mut x, mut y, mut z) = (x.clone(), y.clone(), z.clone());
for _ in 0..(wp * 2 + 200) {
if all_close(&[&x, &y, &z], carlson_bits(wp), wp, rm) {
let one = BigFloat::from_i32(1, wp);
let mu = x
.add(&y, wp, rm)
.add(&z, wp, rm)
.div(&BigFloat::from_i32(3, wp), wp, rm);
return Ok(one.div(&mu.sqrt(wp, rm), wp, rm));
}
let (sx, sy, sz) = (x.sqrt(wp, rm), y.sqrt(wp, rm), z.sqrt(wp, rm));
let lambda =
sx.mul(&sy, wp, rm)
.add(&sy.mul(&sz, wp, rm), wp, rm)
.add(&sz.mul(&sx, wp, rm), wp, rm);
x = x.add(&lambda, wp, rm).div(&four, wp, rm);
y = y.add(&lambda, wp, rm).div(&four, wp, rm);
z = z.add(&lambda, wp, rm).div(&four, wp, rm);
}
Err(SymplexError::ComputationFailed {
operation: "evalf",
reason: "Carlson R_F did not converge".into(),
})
}
fn carlson_rd(
x: &BigFloat,
y: &BigFloat,
z: &BigFloat,
wp: usize,
rm: RoundingMode,
) -> Result<BigFloat, SymplexError> {
let one = BigFloat::from_i32(1, wp);
let three = BigFloat::from_i32(3, wp);
let four = BigFloat::from_i32(4, wp);
let (mut x, mut y, mut z) = (x.clone(), y.clone(), z.clone());
let mut acc = BigFloat::new(wp);
let mut scale = one.clone(); for _ in 0..(wp * 2 + 200) {
if all_close(&[&x, &y, &z], carlson_bits(wp), wp, rm) {
let mu = x.add(&y, wp, rm).add(&three.mul(&z, wp, rm), wp, rm).div(
&BigFloat::from_i32(5, wp),
wp,
rm,
);
let tail = one.div(&mu.mul(&mu.sqrt(wp, rm), wp, rm), wp, rm);
return Ok(acc.add(&scale.mul(&tail, wp, rm), wp, rm));
}
let (sx, sy, sz) = (x.sqrt(wp, rm), y.sqrt(wp, rm), z.sqrt(wp, rm));
let lambda =
sx.mul(&sy, wp, rm)
.add(&sy.mul(&sz, wp, rm), wp, rm)
.add(&sz.mul(&sx, wp, rm), wp, rm);
let term = three.div(&sz.mul(&z.add(&lambda, wp, rm), wp, rm), wp, rm);
acc = acc.add(&scale.mul(&term, wp, rm), wp, rm);
scale = scale.div(&four, wp, rm);
x = x.add(&lambda, wp, rm).div(&four, wp, rm);
y = y.add(&lambda, wp, rm).div(&four, wp, rm);
z = z.add(&lambda, wp, rm).div(&four, wp, rm);
}
Err(SymplexError::ComputationFailed {
operation: "evalf",
reason: "Carlson R_D did not converge".into(),
})
}
fn carlson_rj(
x: &BigFloat,
y: &BigFloat,
z: &BigFloat,
p: &BigFloat,
wp: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
if !p.is_positive() {
return Err(unevaluable(
"elliptic_pi: the characteristic n must satisfy n < 1 for a real value",
));
}
let one = BigFloat::from_i32(1, wp);
let three = BigFloat::from_i32(3, wp);
let four = BigFloat::from_i32(4, wp);
let (mut x, mut y, mut z, mut p) = (x.clone(), y.clone(), z.clone(), p.clone());
let mut acc = BigFloat::new(wp);
let mut scale = one.clone();
for _ in 0..(wp * 2 + 200) {
if all_close(&[&x, &y, &z, &p], carlson_bits(wp), wp, rm) {
let two_p = p.mul(&BigFloat::from_i32(2, wp), wp, rm);
let mu = x.add(&y, wp, rm).add(&z, wp, rm).add(&two_p, wp, rm).div(
&BigFloat::from_i32(5, wp),
wp,
rm,
);
let tail = one.div(&mu.mul(&mu.sqrt(wp, rm), wp, rm), wp, rm);
return Ok(acc.add(&scale.mul(&tail, wp, rm), wp, rm));
}
let (sx, sy, sz, sp) = (
x.sqrt(wp, rm),
y.sqrt(wp, rm),
z.sqrt(wp, rm),
p.sqrt(wp, rm),
);
let lambda =
sx.mul(&sy, wp, rm)
.add(&sy.mul(&sz, wp, rm), wp, rm)
.add(&sz.mul(&sx, wp, rm), wp, rm);
let alpha = p.mul(&sx.add(&sy, wp, rm).add(&sz, wp, rm), wp, rm).add(
&sx.mul(&sy, wp, rm).mul(&sz, wp, rm),
wp,
rm,
);
let beta = sp.mul(&p.add(&lambda, wp, rm), wp, rm);
let rc = carlson_rc(
&alpha.mul(&alpha, wp, rm),
&beta.mul(&beta, wp, rm),
wp,
rm,
cc,
);
acc = acc.add(&scale.mul(&three.mul(&rc, wp, rm), wp, rm), wp, rm);
scale = scale.div(&four, wp, rm);
x = x.add(&lambda, wp, rm).div(&four, wp, rm);
y = y.add(&lambda, wp, rm).div(&four, wp, rm);
z = z.add(&lambda, wp, rm).div(&four, wp, rm);
p = p.add(&lambda, wp, rm).div(&four, wp, rm);
}
Err(SymplexError::ComputationFailed {
operation: "evalf",
reason: "Carlson R_J did not converge".into(),
})
}
fn one_minus_m(
m: &BigFloat,
wp: usize,
rm: RoundingMode,
what: &str,
) -> Result<BigFloat, SymplexError> {
let one = BigFloat::from_i32(1, wp);
let v = one.sub(m, wp, rm);
if v.is_negative() {
return Err(unevaluable(format!(
"{what}: parameter m must satisfy m ≤ 1 for a real value"
)));
}
Ok(v)
}
fn arb_elliptic_k(
m: &BigFloat,
prec: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
if m.is_nan() || m.is_inf() {
return Err(unevaluable("elliptic_k of special float value"));
}
let wp = prec + 32;
let y = one_minus_m(m, wp, rm, "elliptic_k")?;
if y.is_zero() {
return Err(unevaluable("elliptic_k(1) is infinite"));
}
let _ = cc;
let r = carlson_rf(&BigFloat::new(wp), &y, &BigFloat::from_i32(1, wp), wp, rm)?;
Ok(round_to(r, prec, rm))
}
fn arb_elliptic_e(
m: &BigFloat,
prec: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
if m.is_nan() || m.is_inf() {
return Err(unevaluable("elliptic_e of special float value"));
}
let wp = prec + 32;
let one = BigFloat::from_i32(1, wp);
let y = one_minus_m(m, wp, rm, "elliptic_e")?;
if y.is_zero() {
return Ok(BigFloat::from_i32(1, prec));
}
let _ = cc;
let zero = BigFloat::new(wp);
let rf = carlson_rf(&zero, &y, &one, wp, rm)?;
let rd = carlson_rd(&zero, &y, &one, wp, rm)?;
let three = BigFloat::from_i32(3, wp);
let r = rf.sub(&m.mul(&rd, wp, rm).div(&three, wp, rm), wp, rm);
Ok(round_to(r, prec, rm))
}
fn arb_elliptic_f(
phi: &BigFloat,
m: &BigFloat,
prec: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
if phi.is_nan() || phi.is_inf() || m.is_nan() || m.is_inf() {
return Err(unevaluable("elliptic_f of special float value"));
}
if phi.is_zero() {
return Ok(BigFloat::new(prec));
}
let wp = prec + 32;
let one = BigFloat::from_i32(1, wp);
let pi = cc.pi(wp, rm).clone();
let k_bf = phi
.div(&pi, wp, rm)
.add(&BigFloat::from_f64(0.5, wp), wp, rm)
.floor();
let phi0 = phi.sub(&k_bf.mul(&pi, wp, rm), wp, rm);
let s = phi0.sin(wp, rm, cc);
let c = phi0.cos(wp, rm, cc);
let s2 = s.mul(&s, wp, rm);
let c2 = c.mul(&c, wp, rm);
let y = one.sub(&m.mul(&s2, wp, rm), wp, rm);
if !y.is_positive() {
return Err(unevaluable(
"elliptic_f: 1 − m sin²φ must be positive for a real value",
));
}
let mut r = s.mul(&carlson_rf(&c2, &y, &one, wp, rm)?, wp, rm);
if !k_bf.is_zero() {
let k = arb_elliptic_k(m, wp, rm, cc)?;
let two = BigFloat::from_i32(2, wp);
r = r.add(&two.mul(&k_bf, wp, rm).mul(&k, wp, rm), wp, rm);
}
Ok(round_to(r, prec, rm))
}
fn arb_elliptic_pi(
n: &BigFloat,
m: &BigFloat,
prec: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
if n.is_nan() || n.is_inf() || m.is_nan() || m.is_inf() {
return Err(unevaluable("elliptic_pi of special float value"));
}
let wp = prec + 32;
let one = BigFloat::from_i32(1, wp);
let zero = BigFloat::new(wp);
let y = one_minus_m(m, wp, rm, "elliptic_pi")?;
if y.is_zero() {
return Err(unevaluable("elliptic_pi(n, 1) is infinite"));
}
let p = one.sub(n, wp, rm);
let rf = carlson_rf(&zero, &y, &one, wp, rm)?;
if n.is_zero() {
return Ok(round_to(rf, prec, rm));
}
let rj = carlson_rj(&zero, &y, &one, &p, wp, rm, cc)?;
let three = BigFloat::from_i32(3, wp);
let r = rf.add(&n.mul(&rj, wp, rm).div(&three, wp, rm), wp, rm);
Ok(round_to(r, prec, rm))
}
fn poly_degree(
n: &BigFloat,
what: &str,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<u64, SymplexError> {
match bf_as_int(n, rm, cc)? {
Some(d) if (0..=10_000_000).contains(&d) => Ok(d as u64),
_ => Err(unevaluable(format!(
"{what}: degree must be a non-negative integer"
))),
}
}
fn arb_gegenbauer(
n: &BigFloat,
a: &BigFloat,
x: &BigFloat,
prec: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
let deg = poly_degree(n, "gegenbauer", rm, cc)?;
let wp = prec + 32 + (deg as f64).log2().ceil().max(0.0) as usize;
let one = BigFloat::from_i32(1, wp);
let two = BigFloat::from_i32(2, wp);
if deg == 0 {
return Ok(BigFloat::from_i32(1, prec));
}
let mut prev = one.clone();
let mut curr = two.mul(a, wp, rm).mul(x, wp, rm);
for k in 1..deg {
let k_bf = BigFloat::from_i128(k as i128, wp);
let t1 = two
.mul(&k_bf.add(a, wp, rm), wp, rm)
.mul(x, wp, rm)
.mul(&curr, wp, rm);
let coeff = k_bf.add(&two.mul(a, wp, rm), wp, rm).sub(&one, wp, rm);
let t2 = coeff.mul(&prev, wp, rm);
let next = t1
.sub(&t2, wp, rm)
.div(&BigFloat::from_i128(k as i128 + 1, wp), wp, rm);
prev = curr;
curr = next;
}
Ok(round_to(curr, prec, rm))
}
fn arb_jacobi(
n: &BigFloat,
a: &BigFloat,
b: &BigFloat,
x: &BigFloat,
prec: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
let deg = poly_degree(n, "jacobi", rm, cc)? as usize;
if deg > 100_000 {
return Err(unevaluable("jacobi: degree too large for evalf"));
}
let wp = prec + 32 + 2 * deg;
let one = BigFloat::from_i32(1, wp);
let two = BigFloat::from_i32(2, wp);
if deg == 0 {
return Ok(BigFloat::from_i32(1, prec));
}
let lo = x.sub(&one, wp, rm).div(&two, wp, rm);
let hi = x.add(&one, wp, rm).div(&two, wp, rm);
let shifted = |p: &BigFloat, lo_j: usize, hi_j: usize| -> BigFloat {
let mut acc = one.clone();
for j in lo_j..=hi_j {
acc = acc.mul(&p.add(&BigFloat::from_i128(j as i128, wp), wp, rm), wp, rm);
}
acc
};
let mut sum = BigFloat::new(wp);
for s in 0..=deg {
let pa = shifted(a, s + 1, deg);
let pb = shifted(b, deg - s + 1, deg);
let denom = ratio_to_bigfloat(&(factorial_bigint(deg - s) * factorial_bigint(s)), wp, rm);
let term = pa
.mul(&pb, wp, rm)
.div(&denom, wp, rm)
.mul(&lo.powi(s, wp, rm), wp, rm)
.mul(&hi.powi(deg - s, wp, rm), wp, rm);
sum = sum.add(&term, wp, rm);
}
Ok(round_to(sum, prec, rm))
}
fn factorial_bigint(n: usize) -> Ratio<BigInt> {
let mut f = BigInt::from(1);
for i in 2..=n {
f *= BigInt::from(i as u64);
}
Ratio::from_integer(f)
}
fn arb_assoc_legendre(
n: &BigFloat,
m: &BigFloat,
x: &BigFloat,
prec: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
let deg = poly_degree(n, "assoc_legendre", rm, cc)? as i64;
let order = bf_as_int(m, rm, cc)?
.ok_or_else(|| unevaluable("assoc_legendre: order must be an integer"))?;
if order.abs() > deg {
return Ok(BigFloat::new(prec));
}
let wp = prec + 32 + 2 * (deg as f64).log2().ceil().max(0.0) as usize;
let one = BigFloat::from_i32(1, wp);
let mu = order.unsigned_abs() as usize;
let x2 = x.mul(x, wp, rm);
let one_minus_x2 = one.sub(&x2, wp, rm);
let root_pow = if mu.is_multiple_of(2) {
one_minus_x2.powi(mu / 2, wp, rm)
} else {
if one_minus_x2.is_negative() {
return Err(unevaluable(
"assoc_legendre with odd order is complex for |x| > 1",
));
}
one_minus_x2
.sqrt(wp, rm)
.mul(&one_minus_x2.powi(mu / 2, wp, rm), wp, rm)
};
let mut dfact = BigInt::from(1);
for j in 1..=mu {
dfact *= BigInt::from((2 * j - 1) as u64);
}
if mu % 2 == 1 {
dfact = -dfact;
}
let mut prev = ratio_to_bigfloat(&Ratio::from_integer(dfact), wp, rm).mul(&root_pow, wp, rm);
let deg_u = deg as usize;
let mut curr = if deg_u > mu {
BigFloat::from_i128(2 * mu as i128 + 1, wp)
.mul(x, wp, rm)
.mul(&prev, wp, rm)
} else {
prev.clone()
};
for k in (mu + 1)..deg_u {
let t1 = BigFloat::from_i128(2 * k as i128 + 1, wp)
.mul(x, wp, rm)
.mul(&curr, wp, rm);
let t2 = BigFloat::from_i128((k + mu) as i128, wp).mul(&prev, wp, rm);
let next = t1
.sub(&t2, wp, rm)
.div(&BigFloat::from_i128((k - mu + 1) as i128, wp), wp, rm);
prev = curr;
curr = next;
}
if order < 0 {
let mut ratio = factorial_bigint(deg_u - mu) / factorial_bigint(deg_u + mu);
if mu % 2 == 1 {
ratio = -ratio;
}
curr = curr.mul(&ratio_to_bigfloat(&ratio, wp, rm), wp, rm);
}
Ok(round_to(curr, prec, rm))
}
fn arb_assoc_laguerre(
n: &BigFloat,
a: &BigFloat,
x: &BigFloat,
prec: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<BigFloat, SymplexError> {
let deg = poly_degree(n, "assoc_laguerre", rm, cc)?;
let wp = prec + 32 + (deg as f64).log2().ceil().max(0.0) as usize;
let one = BigFloat::from_i32(1, wp);
if deg == 0 {
return Ok(BigFloat::from_i32(1, prec));
}
let mut prev = one.clone();
let mut curr = one.add(a, wp, rm).sub(x, wp, rm);
for k in 1..deg {
let k_bf = BigFloat::from_i128(k as i128, wp);
let coeff = BigFloat::from_i128(2 * k as i128 + 1, wp)
.add(a, wp, rm)
.sub(x, wp, rm);
let t1 = coeff.mul(&curr, wp, rm);
let t2 = k_bf.add(a, wp, rm).mul(&prev, wp, rm);
let next = t1
.sub(&t2, wp, rm)
.div(&BigFloat::from_i128(k as i128 + 1, wp), wp, rm);
prev = curr;
curr = next;
}
Ok(round_to(curr, prec, rm))
}
fn get_cached(cache: &FxHashMap<ExprId, Complex>, id: ExprId) -> Result<&Complex, SymplexError> {
cache.get(&id).ok_or_else(|| SymplexError::Unevaluable {
reason: format!("sub-expression {id:?} not in cache (likely contains free symbols)"),
})
}
fn ratio_to_bigfloat(r: &Ratio<BigInt>, prec: usize, rm: RoundingMode) -> BigFloat {
crate::base::numeric::ratio_to_bigfloat(r, prec, rm)
}
fn try_as_small_int(arena: &Arena, id: ExprId) -> Option<i32> {
let r = arena.as_num(id)?;
if !r.is_integer() {
return None;
}
let n = r.to_integer();
let val: i64 = n.try_into().ok()?;
if val >= i32::MIN as i64 && val <= i32::MAX as i64 {
Some(val as i32)
} else {
None
}
}
fn is_negligible_part(part: &BigFloat, other: &BigFloat, digits: u32) -> bool {
if part.is_zero() {
return true;
}
if other.is_zero() {
return false;
}
match (part.exponent(), other.exponent()) {
(Some(p_exp), Some(o_exp)) => {
let bit_threshold = (digits as i64) * 34 / 10 + 4;
(o_exp as i64 - p_exp as i64) > bit_threshold
}
_ => false,
}
}
fn format_complex(
z: &Complex,
digits: u32,
_prec: usize,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<String, SymplexError> {
let re = &z.0;
let im = &z.1;
let re_negligible = is_negligible_part(re, im, digits);
let im_negligible = is_negligible_part(im, re, digits);
if re_negligible && im_negligible {
return Ok("0".to_string());
}
if im_negligible {
return format_decimal(re, digits, rm, cc);
}
if re_negligible {
let im_str = format_decimal(im, digits, rm, cc)?;
return if im_str == "1" {
Ok("i".to_string())
} else if im_str == "-1" {
Ok("-i".to_string())
} else {
Ok(format!("{im_str}*i"))
};
}
let re_str = format_decimal(re, digits, rm, cc)?;
if im.is_negative() {
let im_abs_str = format_decimal(&im.abs(), digits, rm, cc)?;
if im_abs_str == "1" {
Ok(format!("{re_str} - i"))
} else {
Ok(format!("{re_str} - {im_abs_str}*i"))
}
} else {
let im_str = format_decimal(im, digits, rm, cc)?;
if im_str == "1" {
Ok(format!("{re_str} + i"))
} else {
Ok(format!("{re_str} + {im_str}*i"))
}
}
}
fn round_decimal_digits(mantissa: &[u8], digits: usize) -> (Vec<u8>, i64) {
if digits == 0 || mantissa.len() <= digits {
return (mantissa.to_vec(), 0);
}
let mut kept: Vec<u8> = mantissa[..digits].to_vec();
let next = mantissa[digits];
let rest_nonzero = mantissa[digits + 1..].iter().any(|&d| d != 0);
let round_up = match next.cmp(&5) {
std::cmp::Ordering::Greater => true,
std::cmp::Ordering::Less => false,
std::cmp::Ordering::Equal => rest_nonzero || kept[digits - 1] % 2 == 1,
};
if !round_up {
return (kept, 0);
}
for d in kept.iter_mut().rev() {
if *d == 9 {
*d = 0;
} else {
*d += 1;
return (kept, 0);
}
}
kept[0] = 1;
(kept, 1)
}
fn format_decimal(
val: &BigFloat,
digits: u32,
rm: RoundingMode,
cc: &mut Consts,
) -> Result<String, SymplexError> {
if val.is_zero() {
return Ok("0".to_string());
}
if val.is_inf_pos() {
return Ok("oo".to_string());
}
if val.is_inf_neg() {
return Ok("-oo".to_string());
}
let (sign, mantissa, exponent) = val
.convert_to_radix(Radix::Dec, rm, cc)
.map_err(|e| SymplexError::NotImplemented(format!("decimal conversion failed: {e:?}")))?;
let sign_str = if sign == Sign::Neg { "-" } else { "" };
let (rounded, exp_carry) = round_decimal_digits(&mantissa, digits as usize);
let n = rounded.len();
if n == 0 {
return Ok("0".to_string());
}
let digit_chars: Vec<char> = rounded.iter().map(|&d| (b'0' + d) as char).collect();
let adj_exp = exponent as i64 - 1 + exp_carry;
if adj_exp >= 0 && (adj_exp as usize) < n {
let dot_pos = adj_exp as usize + 1;
let mut s = String::with_capacity(n + 2);
s.push_str(sign_str);
for (i, &ch) in digit_chars.iter().enumerate() {
if i == dot_pos && dot_pos < n {
s.push('.');
}
s.push(ch);
}
if s.contains('.') {
let trimmed = s.trim_end_matches('0').trim_end_matches('.');
Ok(trimmed.to_string())
} else {
Ok(s)
}
} else if adj_exp >= 0 && (adj_exp as usize) < n + 6 {
let mut s = String::with_capacity(adj_exp as usize + 2);
s.push_str(sign_str);
for &ch in &digit_chars {
s.push(ch);
}
let needed_zeros = (adj_exp as usize + 1).saturating_sub(digit_chars.len());
for _ in 0..needed_zeros {
s.push('0');
}
Ok(s)
} else if (-4..0).contains(&adj_exp) {
let mut s = String::with_capacity(n + (-adj_exp) as usize + 3);
s.push_str(sign_str);
s.push_str("0.");
let leading_zeros = (-adj_exp - 1) as usize;
for _ in 0..leading_zeros {
s.push('0');
}
for &ch in &digit_chars {
s.push(ch);
}
let trimmed = s.trim_end_matches('0');
Ok(trimmed.to_string())
} else {
let mut s = String::with_capacity(n + 10);
s.push_str(sign_str);
s.push(digit_chars[0]);
if n > 1 {
s.push('.');
for &ch in &digit_chars[1..] {
s.push(ch);
}
while s.ends_with('0') {
s.pop();
}
if s.ends_with('.') {
s.pop();
}
}
s.push_str(&format!("e{adj_exp}"));
Ok(s)
}
}
#[cfg(test)]
mod tests {
use super::*;
use crate::base::arena::Arena;
fn assert_evalf_starts_with(arena: &Arena, expr: ExprId, digits: u32, prefix: &str) {
let result = evalf(arena, expr, digits).unwrap();
assert!(
result.starts_with(prefix),
"expected to start with '{prefix}', got: '{result}'"
);
}
fn assert_evalf_eq(arena: &Arena, expr: ExprId, digits: u32, expected: &str) {
let result = evalf(arena, expr, digits).unwrap();
assert_eq!(result, expected, "evalf mismatch");
}
#[test]
fn integer_zero() {
let a = Arena::new();
assert_evalf_eq(&a, a.zero, 10, "0");
}
#[test]
fn integer_one() {
let a = Arena::new();
assert_evalf_eq(&a, a.one, 10, "1");
}
#[test]
fn integer_negative_one() {
let a = Arena::new();
assert_evalf_eq(&a, a.neg_one, 10, "-1");
}
#[test]
fn integer_five() {
let mut a = Arena::new();
let five = a.int(5);
assert_evalf_eq(&a, five, 10, "5");
}
#[test]
fn integer_large() {
let mut a = Arena::new();
let n = a.int(123456);
assert_evalf_eq(&a, n, 10, "123456");
}
#[test]
fn rational_half() {
let mut a = Arena::new();
let half = a.rational(1, 2);
assert_evalf_eq(&a, half, 10, "0.5");
}
#[test]
fn rational_one_third() {
let mut a = Arena::new();
let third = a.rational(1, 3);
assert_evalf_starts_with(&a, third, 10, "0.333333333");
}
#[test]
fn rational_negative() {
let mut a = Arena::new();
let r = a.rational(-3, 4);
assert_evalf_starts_with(&a, r, 10, "-0.75");
}
#[test]
fn constant_pi_15_digits() {
let a = Arena::new();
assert_evalf_starts_with(&a, a.pi, 15, "3.14159265358979");
}
#[test]
fn constant_pi_50_digits() {
let a = Arena::new();
let result = evalf(&a, a.pi, 50).unwrap();
assert!(
result.starts_with("3.1415926535897932384626433832795"),
"pi to 50 digits: {result}"
);
}
#[test]
fn constant_e_15_digits() {
let a = Arena::new();
assert_evalf_starts_with(&a, a.e_const, 15, "2.71828182845905");
}
#[test]
fn add_integers() {
let mut a = Arena::new();
let two = a.int(2);
let three = a.int(3);
let expr = a.add(&[two, three]);
assert_evalf_eq(&a, expr, 10, "5");
}
#[test]
fn mul_integers() {
let mut a = Arena::new();
let two = a.int(2);
let three = a.int(3);
let expr = a.mul(&[two, three]);
assert_evalf_eq(&a, expr, 10, "6");
}
#[test]
fn pow_integer() {
let mut a = Arena::new();
let two = a.int(2);
let ten = a.int(10);
let expr = a.pow(two, ten);
assert_evalf_eq(&a, expr, 10, "1024");
}
#[test]
fn pow_negative_exponent() {
let mut a = Arena::new();
let two = a.int(2);
let neg_one = a.int(-1);
let expr = a.pow(two, neg_one);
assert_evalf_eq(&a, expr, 10, "0.5");
}
#[test]
fn pi_plus_one() {
let mut a = Arena::new();
let pi = a.pi;
let one = a.one;
let expr = a.add(&[pi, one]);
assert_evalf_starts_with(&a, expr, 15, "4.14159265358979");
}
#[test]
fn sin_zero() {
let mut a = Arena::new();
let zero = a.zero;
let expr = a.sin(zero);
assert_evalf_eq(&a, expr, 10, "0");
}
#[test]
fn cos_zero() {
let mut a = Arena::new();
let zero = a.zero;
let expr = a.cos(zero);
assert_evalf_eq(&a, expr, 10, "1");
}
#[test]
fn exp_zero() {
let mut a = Arena::new();
let zero = a.zero;
let expr = a.exp(zero);
assert_evalf_eq(&a, expr, 10, "1");
}
#[test]
fn exp_one() {
let mut a = Arena::new();
let one = a.one;
let expr = a.exp(one);
assert_evalf_starts_with(&a, expr, 15, "2.71828182845905");
}
#[test]
fn ln_one() {
let mut a = Arena::new();
let one = a.one;
let expr = a.ln(one);
assert_evalf_eq(&a, expr, 10, "0");
}
#[test]
fn ln_e() {
let mut a = Arena::new();
let e = a.e_const;
let expr = a.ln(e);
let result = evalf(&a, expr, 15).unwrap();
let val: f64 = result.parse().unwrap_or(999.0);
assert!(
(val - 1.0).abs() < 1e-10,
"ln(e) should be ~1, got: {result}"
);
}
#[test]
fn sqrt_four() {
let mut a = Arena::new();
let four = a.int(4);
let expr = a.sqrt(four);
assert_evalf_eq(&a, expr, 10, "2");
}
#[test]
fn sqrt_two_30_digits() {
let mut a = Arena::new();
let two = a.int(2);
let expr = a.sqrt(two);
assert_evalf_starts_with(&a, expr, 30, "1.4142135623730950488");
}
#[test]
fn sin_pi_is_near_zero() {
let mut a = Arena::new();
let pi = a.pi;
let expr = a.sin(pi);
let result = evalf(&a, expr, 20).unwrap();
let val: f64 = result.parse().unwrap_or(999.0);
assert!(val.abs() < 1e-15, "sin(pi) should be ~0, got: {result}");
}
#[test]
fn neg_pi() {
let mut a = Arena::new();
let pi = a.pi;
let expr = a.neg(pi);
assert_evalf_starts_with(&a, expr, 10, "-3.14159265");
}
#[test]
fn abs_neg_pi() {
let mut a = Arena::new();
let pi = a.pi;
let neg_pi = a.neg(pi);
let expr = a.abs(neg_pi);
assert_evalf_starts_with(&a, expr, 10, "3.14159265");
}
#[test]
fn abs_negative_integer() {
let mut a = Arena::new();
let neg_five = a.int(-5);
let expr = a.abs(neg_five);
assert_evalf_eq(&a, expr, 10, "5");
}
#[test]
fn free_symbol_errors() {
let mut a = Arena::new();
let x = a.symbol("x");
let result = evalf(&a, x, 10);
assert!(result.is_err(), "evalf of free symbol should error");
}
#[test]
fn infinity_errors() {
let a = Arena::new();
let result = evalf(&a, a.infinity, 10);
assert!(result.is_err(), "evalf of infinity should error");
}
#[test]
fn nan_errors() {
let a = Arena::new();
let result = evalf(&a, a.nan, 10);
assert!(result.is_err(), "evalf of NaN should error");
}
#[test]
fn imaginary_unit_works() {
let a = Arena::new();
let result = evalf(&a, a.i_unit, 10);
assert!(result.is_ok(), "evalf of i should work now");
let s = result.unwrap();
assert!(
s.contains("i") || s.contains("I"),
"should display as imaginary: {s}"
);
}
#[test]
fn evalf_one_plus_i() {
let mut a = Arena::new();
let one = a.one;
let i = a.i_unit;
let expr = a.add(&[one, i]);
let result = evalf(&a, expr, 10);
assert!(result.is_ok(), "evalf(1+i) should work: {:?}", result.err());
let s = result.unwrap();
assert!(
s.contains("i") || s.contains("I"),
"1+i should contain 'i': {s}"
);
}
#[test]
fn evalf_exp_i_pi() {
let mut a = Arena::new();
let i = a.i_unit;
let pi = a.pi;
let i_times_pi = a.mul(&[i, pi]);
let expr = a.exp(i_times_pi);
let result = evalf(&a, expr, 15);
assert!(
result.is_ok(),
"evalf(exp(i*pi)) should work: {:?}",
result.err()
);
let s = result.unwrap();
let val: f64 = s.parse().unwrap_or(999.0);
assert!(
(val - (-1.0)).abs() < 1e-10,
"exp(i*pi) should be ~-1, got: {s}"
);
}
#[test]
fn evalf_i_squared() {
let mut a = Arena::new();
let i = a.i_unit;
let two = a.int(2);
let expr = a.pow(i, two);
let result = evalf(&a, expr, 10);
assert!(result.is_ok(), "evalf(i^2) should work: {:?}", result.err());
let s = result.unwrap();
let val: f64 = s.parse().unwrap_or(999.0);
assert!((val - (-1.0)).abs() < 1e-10, "i^2 should be -1, got: {s}");
}
#[test]
fn evalf_sqrt_neg_one() {
let mut a = Arena::new();
let neg_one = a.neg_one;
let expr = a.sqrt(neg_one);
let result = evalf(&a, expr, 10);
assert!(
result.is_ok(),
"evalf(sqrt(-1)) should work: {:?}",
result.err()
);
let s = result.unwrap();
assert!(
s.contains("i") || s.contains("I"),
"sqrt(-1) should be imaginary: {s}"
);
}
#[test]
fn evalf_abs_3_plus_4i() {
let mut a = Arena::new();
let three = a.int(3);
let four = a.int(4);
let i = a.i_unit;
let four_i = a.mul(&[four, i]);
let sum = a.add(&[three, four_i]);
let expr = a.abs(sum);
let result = evalf(&a, expr, 10);
assert!(
result.is_ok(),
"evalf(|3+4i|) should work: {:?}",
result.err()
);
let s = result.unwrap();
let val: f64 = s.parse().unwrap_or(999.0);
assert!((val - 5.0).abs() < 1e-10, "|3+4i| should be 5, got: {s}");
}
#[test]
fn polynomial_at_value() {
let mut a = Arena::new();
let x = a.symbol("x");
let two = a.int(2);
let x_sq = a.pow(x, two);
let two_x = a.mul(&[two, x]);
let expr = a.add(&[x_sq, two_x, a.one]);
let three = a.int(3);
let evaluated = crate::transforms::subs::subs(&mut a, expr, x, three);
assert_evalf_eq(&a, evaluated, 10, "16");
}
#[test]
fn sin_cos_identity() {
let mut a = Arena::new();
let one = a.one;
let two = a.int(2);
let sin_1 = a.sin(one);
let cos_1 = a.cos(one);
let sin_sq = a.pow(sin_1, two);
let cos_sq = a.pow(cos_1, two);
let expr = a.add(&[sin_sq, cos_sq]);
let result = evalf(&a, expr, 15).unwrap();
let val: f64 = result.parse().unwrap_or(0.0);
assert!(
(val - 1.0).abs() < 1e-12,
"sin^2(1) + cos^2(1) should be ~1, got: {result}"
);
}
#[test]
fn same_precision_same_result() {
let a = Arena::new();
let r1 = evalf(&a, a.pi, 20).unwrap();
let r2 = evalf(&a, a.pi, 20).unwrap();
assert_eq!(r1, r2, "same precision should give same result");
}
#[test]
fn bessel_j0_at_zero() {
let mut a = Arena::new();
let zero = a.zero;
let j0_0 = a.besselj(zero, zero);
let result = evalf(&a, j0_0, 15).unwrap();
let val: f64 = result.parse().unwrap_or(f64::NAN);
assert!(
(val - 1.0).abs() < 1e-10,
"J_0(0) should be 1.0, got {result}"
);
}
#[test]
fn bessel_j0_at_one() {
let mut a = Arena::new();
let zero = a.zero;
let one = a.one;
let j0_1 = a.besselj(zero, one);
let result = evalf(&a, j0_1, 15).unwrap();
let val: f64 = result.parse().unwrap_or(f64::NAN);
assert!(
(val - 0.7651976865579666).abs() < 1e-8,
"J_0(1) should be ~0.7652, got {result}"
);
}
#[test]
fn bessel_j1_at_zero() {
let mut a = Arena::new();
let one = a.one;
let zero = a.zero;
let j1_0 = a.besselj(one, zero);
let result = evalf(&a, j1_0, 15).unwrap();
let val: f64 = result.parse().unwrap_or(f64::NAN);
assert!(val.abs() < 1e-10, "J_1(0) should be 0, got {result}");
}
#[test]
fn bessel_j1_at_one() {
let mut a = Arena::new();
let one = a.one;
let j1_1 = a.besselj(one, one);
let result = evalf(&a, j1_1, 15).unwrap();
let val: f64 = result.parse().unwrap_or(f64::NAN);
assert!(
(val - 0.44005058574493355).abs() < 1e-8,
"J_1(1) should be ~0.4401, got {result}"
);
}
#[test]
fn bessel_y0_at_one() {
let mut a = Arena::new();
let zero = a.zero;
let one = a.one;
let y0_1 = a.bessely(zero, one);
let result = evalf(&a, y0_1, 15).unwrap();
let val: f64 = result.parse().unwrap_or(f64::NAN);
assert!(
(val - 0.08825696421567696).abs() < 1e-6,
"Y_0(1) should be ~0.0883, got {result}"
);
}
#[test]
fn euler_gamma_and_catalan_high_precision() {
let a = Arena::new();
assert_evalf_starts_with(
&a,
a.euler_gamma,
60,
"0.57721566490153286060651209008240243104215933593992359880576",
);
assert_evalf_starts_with(
&a,
a.catalan,
60,
"0.91596559417721901505460351493238411077414937428167213426649",
);
assert_evalf_starts_with(
&a,
a.golden_ratio,
40,
"1.61803398874989484820458683436563811772",
);
}
#[test]
fn bigint_to_bigfloat_is_exact_for_huge_integers() {
let n = (BigInt::from(1) << 200) + BigInt::from(1);
let bf = crate::base::numeric::bigint_to_bigfloat(&n, 256);
let two200 = BigFloat::from_i32(2, 256).powi(200, 256, RoundingMode::ToEven);
let diff = bf.sub(&two200, 256, RoundingMode::ToEven);
assert_eq!(diff, BigFloat::from_i32(1, 256));
}
#[test]
fn bf_pow_handles_exactly_representable_results() {
let mut cc = Consts::new().unwrap();
let rm = RoundingMode::ToEven;
let four = BigFloat::from_i32(4, 160);
let half = BigFloat::from_f64(0.5, 160);
let r = bf_pow(&four, &half, 160, rm, &mut cc);
let diff = r.sub(&BigFloat::from_i32(2, 160), 160, rm);
assert!(
diff.is_zero() || diff.exponent().is_none_or(|e| e < -140),
"4^0.5 = {r}"
);
let three = BigFloat::from_i32(3, 160);
let r = bf_pow(&three, &BigFloat::from_i32(5, 160), 160, rm, &mut cc);
assert_eq!(r, BigFloat::from_i32(243, 160));
let r = bf_pow(&three, &BigFloat::from_i32(-1, 160), 160, rm, &mut cc);
let back = r.mul(&three, 160, rm);
let d = back.sub(&BigFloat::from_i32(1, 160), 160, rm);
assert!(d.is_zero() || d.exponent().is_none_or(|e| e < -140));
}
#[test]
fn si_ci_ei_li_at_one_and_two() {
let mut a = Arena::new();
let one = a.one;
let two = a.int(2);
let si = a.si(one);
assert_evalf_starts_with(&a, si, 20, "0.94608307036718301494");
let ci = a.ci(one);
assert_evalf_starts_with(&a, ci, 20, "0.33740392290096813466");
let ei = a.ei(one);
assert_evalf_starts_with(&a, ei, 20, "1.8951178163559367555");
let li = a.li(two);
assert_evalf_starts_with(&a, li, 20, "1.0451637801174927848");
}
#[test]
fn si_asymptotic_and_series_regimes_agree() {
let mut a = Arena::new();
let x = a.int(200);
let si = a.si(x);
let lo = evalf(&a, si, 12).unwrap(); let hi = evalf(&a, si, 120).unwrap(); assert!(hi.starts_with(&lo[..12]), "{lo} vs {hi}");
}
#[test]
fn zeta_borwein_and_functional_equation() {
let mut a = Arena::new();
let three = a.int(3);
let z3 = a.zeta(three);
assert_evalf_starts_with(&a, z3, 30, "1.20205690315959428539973816151");
let half = a.rational(1, 2);
let zh = a.zeta(half);
assert_evalf_starts_with(&a, zh, 20, "-1.4603545088095868129");
let neg = a.rational(-5, 2);
let zn = a.zeta(neg);
assert_evalf_starts_with(&a, zn, 15, "0.00851692877785");
let neg_half = a.rational(-1, 2);
let znh = a.zeta(neg_half);
assert_evalf_starts_with(&a, znh, 15, "-0.20788622497735");
}
#[test]
fn polygamma_matches_trigamma_closed_forms() {
let mut a = Arena::new();
let one = a.one;
let x = a.rational(7, 3);
let pg = a.intern(ExprNode::Polygamma(one, x));
assert_evalf_starts_with(&a, pg, 20, "0.53309712542709408179");
let half = a.rational(1, 2);
let pg_half = a.intern(ExprNode::Polygamma(one, half));
assert_evalf_starts_with(&a, pg_half, 20, "4.9348022005446793094");
}
#[test]
fn bessel_i_k_reference() {
let mut a = Arena::new();
let zero = a.zero;
let one = a.one;
let i0 = a.besseli(zero, one);
assert_evalf_starts_with(&a, i0, 20, "1.2660658777520083356");
let k0 = a.besselk(zero, one);
assert_evalf_starts_with(&a, k0, 20, "0.42102443824070833334");
let k1 = a.besselk(one, one);
assert_evalf_starts_with(&a, k1, 20, "0.60190723019723457474");
}
#[test]
fn orthogonal_polynomials_large_degree() {
let mut a = Arena::new();
let n = a.int(100);
let third = a.rational(1, 3);
let t = a.chebyshev_t(n, third);
let s = evalf(&a, t, 15).unwrap();
let v: f64 = s.parse().unwrap();
let expected = (100.0 * (1.0f64 / 3.0).acos()).cos();
assert!((v - expected).abs() < 1e-12, "{v} vs {expected}");
let p = a.legendre(n, a.one);
assert_evalf_starts_with(&a, p, 10, "1");
}
#[test]
fn complex_nodes_evaluate_numerically() {
let mut a = Arena::new();
let z = a.symbol("z");
let re_z = a.intern(ExprNode::Re(z));
let arg_z = a.intern(ExprNode::Arg(z));
let conj_z = a.intern(ExprNode::Conjugate(z));
let three = a.int(3);
let four = a.int(4);
let four_i = a.mul(&[four, a.i_unit]);
let w = a.add(&[three, four_i]);
let _ = re_z;
let re_w = a.intern(ExprNode::Re(w));
let arg_w = a.intern(ExprNode::Arg(w));
let conj_w = a.intern(ExprNode::Conjugate(w));
let _ = (arg_z, conj_z);
assert_eq!(evalf(&a, re_w, 10).unwrap(), "3");
assert_evalf_starts_with(&a, arg_w, 15, "0.92729521800161");
let s = evalf(&a, conj_w, 10).unwrap();
assert!(s.contains('3') && s.contains('4') && s.contains('-'), "{s}");
}
fn apply(a: &mut Arena, f: LibFn, args: &[ExprId]) -> ExprId {
a.lib_apply(f, args)
}
#[test]
fn erfc_node_large_argument_is_relative_precision() {
let mut a = Arena::new();
let seven = a.int(7);
let e7 = a.erfc(seven);
assert_evalf_starts_with(&a, e7, 16, "4.183825607779414e-23");
let ten = a.int(10);
let e10 = a.erfc(ten);
assert_evalf_starts_with(
&a,
e10,
50,
"2.0884875837625447570007862949577886115608181193212e-45",
);
let thirty = a.int(30);
let e30 = a.erfc(thirty);
assert_evalf_starts_with(&a, e30, 16, "2.564656203756112e-393");
let m3 = a.int(-3);
let em3 = a.erfc(m3);
assert_evalf_starts_with(&a, em3, 16, "1.999977909503001");
}
#[test]
fn erf_taylor_branch_carries_cancellation_guard_bits() {
let mut a = Arena::new();
let x = a.rational(74, 10);
let e = a.erf(x);
assert_evalf_starts_with(&a, e, 30, "0.999999999999999999999999875161");
let ten = a.int(10);
let e10 = a.erf(ten);
assert_evalf_starts_with(
&a,
e10,
50,
"0.99999999999999999999999999999999999999999999791151",
);
let three = a.int(3);
let e3 = a.erf(three);
assert_evalf_starts_with(&a, e3, 40, "0.999977909503001414558627223870417679620");
}
#[test]
fn polylog_at_minus_one_is_minus_eta_for_all_s() {
let mut a = Arena::new();
let neg_one = a.neg_one;
let half = a.rational(1, 2);
let li = apply(&mut a, LibFn::PolyLog, &[half, neg_one]);
assert_evalf_starts_with(&a, li, 30, "-0.60489864342163037024726591423");
let one = a.one;
let li1 = apply(&mut a, LibFn::PolyLog, &[one, neg_one]);
assert_evalf_starts_with(&a, li1, 30, "-0.69314718055994530941723212145");
let neg_half = a.rational(-1, 2);
let lim = apply(&mut a, LibFn::PolyLog, &[neg_half, neg_one]);
assert_evalf_starts_with(&a, lim, 30, "-0.38010481260968401677754215655");
let three_halves = a.rational(3, 2);
let li32 = apply(&mut a, LibFn::PolyLog, &[three_halves, neg_one]);
assert_evalf_starts_with(&a, li32, 30, "-0.76514702462540794536726875860");
let one_id = a.one;
let div = apply(&mut a, LibFn::PolyLog, &[half, one_id]);
assert!(evalf(&a, div, 16).is_err());
}
#[test]
fn dirichlet_eta_near_one_uses_alternating_sum() {
let mut a = Arena::new();
let ten = a.int(10);
let e20 = a.int(-20);
let tiny = a.pow(ten, e20);
let s = a.add(&[a.one, tiny]);
let eta = apply(&mut a, LibFn::DirichletEta, &[s]);
assert_evalf_starts_with(&a, eta, 30, "0.69314718055994530941883081049");
let s2 = a.sub(a.one, tiny);
let eta2 = apply(&mut a, LibFn::DirichletEta, &[s2]);
assert_evalf_starts_with(&a, eta2, 30, "0.69314718055994530941563343242");
let one = a.one;
let eta1 = apply(&mut a, LibFn::DirichletEta, &[one]);
assert_evalf_starts_with(&a, eta1, 30, "0.69314718055994530941723212145");
let neg_half = a.rational(-1, 2);
let etam = apply(&mut a, LibFn::DirichletEta, &[neg_half]);
assert_evalf_starts_with(&a, etam, 30, "0.38010481260968401677754215655");
}
#[test]
fn airy_derivatives_at_exact_zero() {
let mut a = Arena::new();
let zero = a.zero;
let aip = apply(&mut a, LibFn::AiryAiPrime, &[zero]);
assert_evalf_starts_with(&a, aip, 30, "-0.25881940379280679840518356018");
let bip = apply(&mut a, LibFn::AiryBiPrime, &[zero]);
assert_evalf_starts_with(&a, bip, 30, "0.44828835735382635791482371039");
}
#[test]
fn polylog_series_bound_accounts_for_negative_order() {
let n = polylog_series_terms(2.0, 0.5, 100);
assert!((120..=122).contains(&n), "{n}");
let n = polylog_series_terms(-40.5, 0.5, 132);
let a = -(0.5f64).ln();
let k = (n - 20) as f64;
assert!(
k * a - 40.5 * k.ln() >= 132.0 * std::f64::consts::LN_2,
"{n}"
);
assert!(n > 200, "{n}");
}
}