use std::sync::LazyLock;
use crate::{
atom::{Atom, AtomCore, AtomOrView, AtomView, EvaluationInfo, FunctionBuilder, Symbol},
coefficient::{Coefficient, CoefficientView},
domains::{
backend::float::Constant,
float::{Complex, ErrorPropagatingFloat, Float, FloatLike, Real, RealLike, SingleFloat},
integer::Integer,
rational::Rational,
},
function, get_symbol,
state::{State, StateInitializer},
symbol,
utils::Settable,
};
static SPECIALS: LazyLock<SpecialSymbols> = LazyLock::new(|| SpecialSymbols {
euler_gamma: get_symbol!("euler_gamma").expect("Euler gamma not defined"),
gamma: get_symbol!("gamma").expect("gamma not defined"),
erf: get_symbol!("erf").expect("erf not defined"),
polygamma: get_symbol!("polygamma").expect("polygamma not defined"),
polylog: get_symbol!("polylog").expect("polylog not defined"),
zeta: get_symbol!("zeta").expect("zeta not defined"),
});
static GEOMETRICS: LazyLock<GeometricSymbols> = LazyLock::new(|| GeometricSymbols {
tan: get_symbol!("tan").expect("tan not defined"),
cot: get_symbol!("cot").expect("cot not defined"),
sec: get_symbol!("sec").expect("sec not defined"),
csc: get_symbol!("csc").expect("csc not defined"),
asin: get_symbol!("asin").expect("asin not defined"),
acos: get_symbol!("acos").expect("acos not defined"),
atan: get_symbol!("atan").expect("atan not defined"),
acot: get_symbol!("acot").expect("acot not defined"),
asec: get_symbol!("asec").expect("asec not defined"),
acsc: get_symbol!("acsc").expect("acsc not defined"),
sinh: get_symbol!("sinh").expect("sinh not defined"),
cosh: get_symbol!("cosh").expect("cosh not defined"),
tanh: get_symbol!("tanh").expect("tanh not defined"),
coth: get_symbol!("coth").expect("coth not defined"),
sech: get_symbol!("sech").expect("sech not defined"),
csch: get_symbol!("csch").expect("csch not defined"),
asinh: get_symbol!("asinh").expect("asinh not defined"),
acosh: get_symbol!("acosh").expect("acosh not defined"),
atanh: get_symbol!("atanh").expect("atanh not defined"),
acoth: get_symbol!("acoth").expect("acoth not defined"),
asech: get_symbol!("asech").expect("asech not defined"),
acsch: get_symbol!("acsch").expect("acsch not defined"),
});
static BESSELS: LazyLock<BesselSymbols> = LazyLock::new(|| BesselSymbols {
bessel_j: get_symbol!("bessel_j").expect("bessel_j not defined"),
bessel_y: get_symbol!("bessel_y").expect("bessel_y not defined"),
bessel_i: get_symbol!("bessel_i").expect("bessel_i not defined"),
bessel_k: get_symbol!("bessel_k").expect("bessel_k not defined"),
});
struct SpecialSymbols {
euler_gamma: Symbol,
gamma: Symbol,
erf: Symbol,
polygamma: Symbol,
polylog: Symbol,
zeta: Symbol,
}
struct GeometricSymbols {
tan: Symbol,
cot: Symbol,
sec: Symbol,
csc: Symbol,
asin: Symbol,
acos: Symbol,
atan: Symbol,
acot: Symbol,
asec: Symbol,
acsc: Symbol,
sinh: Symbol,
cosh: Symbol,
tanh: Symbol,
coth: Symbol,
sech: Symbol,
csch: Symbol,
asinh: Symbol,
acosh: Symbol,
atanh: Symbol,
acoth: Symbol,
asech: Symbol,
acsch: Symbol,
}
struct BesselSymbols {
bessel_j: Symbol,
bessel_y: Symbol,
bessel_i: Symbol,
bessel_k: Symbol,
}
impl SpecialSymbols {
fn new() -> Self {
let euler_gamma = symbol!(
"γ",
aliases = ["euler_gamma"],
der = |_x, _i, out| {
out.to_num(Coefficient::zero());
},
eval = EvaluationInfo::constant(|_tags, prec| {
Ok(Complex::new(
Float::with_val(prec, Constant::Euler),
Float::new(prec),
))
})
);
let mut symbols = crate::symbol_group!(
"gamma";;
|symbols, b| {
let gamma_symbol = symbols[0];
let polygamma_symbol = symbols[1];
b.with_normalization_function(|x, out| {
if let Some([arg]) = function_arguments::<1>(x) {
if maybe_eval_unary_float_in_norm(arg, out, gamma_numeric_eval) {
return;
}
if let Ok(rat) = Rational::try_from(arg)
&& let Some(exact) = gamma_exact_rational(&rat)
{
**out = exact;
}
}
})
.with_derivative_function(move |x, i, out| {
if i == 0
&& let Some([arg]) = function_arguments::<1>(x)
{
**out = function!(gamma_symbol, arg) * function!(polygamma_symbol, 0, arg);
}
})
.with_series_function(
move |args| {
let [arg] = args else { return None; };
let point = arg.coefficient(0.into());
let Ok(pole) = Integer::try_from(&point) else {
return None;
};
if pole > 0 {
return None;
}
let arg = arg.to_atom();
let shift = (-pole).to_i64().unwrap();
let mut regularized = function!(gamma_symbol, arg.to_owned() + shift + 1);
for k in 0..shift {
regularized /= arg.to_owned() + k;
}
Some((Atom::num(1) / (arg + shift), regularized))
},
)
.with_evaluation_info(
EvaluationInfo::new()
.register(|args: &[Complex<Float>]| {
unary_eval_complex_float(args, gamma_numeric_eval)
})
.register(|args: &[Float]| {
let prec = args.first().map(|x| x.prec()).unwrap_or(53);
let [arg] = args else {
return Float::with_val(53, f64::NAN);
};
gamma_numeric_eval(&Complex::new(arg.clone(), Float::new(prec)), prec).re
})
.register(|args: &[f64]| unary_eval_f64(args, gamma_numeric_eval))
.register(|args: &[Complex<f64>]| {
unary_eval_complex_f64(args, gamma_numeric_eval)
}),
)
},
"polygamma";;
|symbols, b| {
let polygamma_symbol = symbols[1];
b.with_normalization_function(move |x, out| {
if let Some([n, z]) = function_arguments::<2>(x)
&& let Ok(order) = u32::try_from(n) {
if maybe_eval_polygamma_float_in_norm(order, z, out) {
return;
}
if let Ok(rat) = Rational::try_from(z)
&& let Some(exact) = polygamma_exact_rational(order, &rat)
{
**out = exact;
return;
}
if let Ok(rat) = Rational::try_from(z)
&& rat.denominator() == 1
&& rat.numerator() <= 0
{
out.to_num(Coefficient::complex_infinity());
}
}
})
.with_derivative_function(move |x, i, out| {
if i == 1
&& let Some([n, z]) = function_arguments::<2>(x)
&& let Ok(order) = u32::try_from(n)
{
**out = function!(polygamma_symbol, Atom::num(order + 1), z);
}
})
.with_series_function(move |args| {
let [n, arg] = args else { return None; };
let order = u32::try_from(n.coefficient(0.into())).ok()?;
let point = arg.coefficient(0.into());
let Ok(pole) = Integer::try_from(&point) else {
return None;
};
if pole > 0 {
return None;
}
let arg = arg.to_atom();
let shift = (-pole).to_i64().unwrap();
let pole = arg.to_owned() + shift;
let pole_power = pole.pow(Atom::num(order + 1));
let sign = if order % 2 == 0 {
Atom::num(Integer::factorial(order))
} else {
-Atom::num(Integer::factorial(order))
};
let mut recurrence_sum = Atom::new();
for k in 0..=shift {
recurrence_sum +=
Atom::num(1) / (arg.to_owned() + k).pow(Atom::num(order + 1));
}
let regularized = pole_power.clone()
* (function!(polygamma_symbol, Atom::num(order), arg + shift + 1)
- sign * recurrence_sum);
Some((Atom::num(1) / pole_power, regularized))
})
.with_evaluation_info(
EvaluationInfo::new().with_tags(1)
.register_tagged(|tags| {
let order = nonnegative_integer_tag("polygamma", tags).ok();
Box::new(move |args: &[Complex<Float>]| {
let Some(order) = order else {
return Complex::new(
Float::with_val(53, f64::NAN),
Float::with_val(53, f64::NAN),
);
};
unary_eval_complex_float(args, |arg, prec| {
polygamma_numeric_eval(order, arg, prec)
})
})
})
.register_tagged(|tags| {
let order = nonnegative_integer_tag("polygamma", tags).ok();
Box::new(move |args: &[Float]| {
let Some(order) = order else {
return Float::with_val(53, f64::NAN);
};
let prec = args.first().map(|x| x.prec()).unwrap_or(53);
let [arg] = args else {
return Float::with_val(53, f64::NAN);
};
polygamma_numeric_eval(
order,
&Complex::new(arg.clone(), Float::new(prec)),
prec,
)
.re
})
})
.register_tagged(|tags| {
let order = nonnegative_integer_tag("polygamma", tags).ok();
Box::new(move |args: &[f64]| polygamma_eval_f64(order, args))
})
.register_tagged(|tags| {
let order = nonnegative_integer_tag("polygamma", tags).ok();
Box::new(move |args: &[Complex<f64>]| {
polygamma_eval_complex_f64(order, args)
})
}),
)
},
"polylog";;
|symbols, b| {
let polylog_symbol = symbols[2];
b.with_normalization_function(|x, out| {
if let Some([s, z]) = function_arguments::<2>(x) {
if let Some(exact) = polylog_exact(s, z) {
**out = exact;
return;
}
if maybe_eval_binary_float_in_norm(s, z, out, polylog_numeric_eval) {
}
}
})
.with_derivative_function(move |x, i, out| {
if i == 1
&& let Some([s, z]) = function_arguments::<2>(x)
{
**out = function!(polylog_symbol, s.to_owned() - 1, z) / z;
}
})
.with_evaluation_info(
EvaluationInfo::new().with_tags(1)
.register_tagged(|tags| {
let tags = tags.iter().map(|x| x.to_owned()).collect::<Vec<_>>();
Box::new(move |args: &[Complex<Float>]| {
let prec = args
.first()
.map(|x| x.re.prec().max(x.im.prec()))
.unwrap_or(53);
let tag_views = tags.iter().map(|x| x.as_view()).collect::<Vec<_>>();
let Ok(order) = complex_float_tag("polylog", &tag_views, prec) else {
return Complex::new(
Float::with_val(53, f64::NAN),
Float::with_val(53, f64::NAN),
);
};
let [arg] = args else {
return Complex::new(
Float::with_val(53, f64::NAN),
Float::with_val(53, f64::NAN),
);
};
polylog_numeric_eval(&order, arg, prec).unwrap_or_else(|| {
Complex::new(
Float::with_val(53, f64::NAN),
Float::with_val(53, f64::NAN),
)
})
})
})
.register_tagged(|tags| {
let tags = tags.iter().map(|x| x.to_owned()).collect::<Vec<_>>();
Box::new(move |args: &[Float]| {
let prec = args.first().map(|x| x.prec()).unwrap_or(53);
let tag_views = tags.iter().map(|x| x.as_view()).collect::<Vec<_>>();
let Ok(order) = complex_float_tag("polylog", &tag_views, prec) else {
return Float::with_val(53, f64::NAN);
};
let [arg] = args else {
return Float::with_val(53, f64::NAN);
};
polylog_numeric_eval(
&order,
&Complex::new(arg.clone(), Float::new(prec)),
prec,
)
.map(|z| z.re)
.unwrap_or_else(|| Float::with_val(53, f64::NAN))
})
})
.register_tagged(|tags| {
let order = complex_float_tag("polylog", tags, 53).ok();
Box::new(move |args: &[f64]| polylog_eval_f64(order.clone(), args))
})
.register_tagged(|tags| {
let order = complex_float_tag("polylog", tags, 53).ok();
Box::new(move |args: &[Complex<f64>]| {
polylog_eval_complex_f64(order.clone(), args)
})
}),
)
},
"dzeta";;
|symbols, b| {
let dzeta = symbols[3];
let zeta = symbols[4];
b.with_normalization_function(move |x, out| {
if let Some([depth, arg]) = function_arguments::<2>(x) {
if depth == 0 {
**out = function!(zeta, arg);
return;
}
if depth == 1 && arg == 0 {
**out = -(Atom::var(Symbol::PI) * 2).log() / 2;
}
}
})
.with_derivative_function(move |x, _, out| {
if let Some([depth, arg]) = function_arguments::<2>(x)
{
**out = function!(dzeta, depth.to_owned() + 1, arg);
}
})
}
,
"zeta";;
|symbols, b| {
let gamma = symbols[0];
let dzeta = symbols[3];
let zeta_symbol = symbols[4];
b.with_normalization_function(|x, out| {
if let Some([arg]) = function_arguments::<1>(x) {
if maybe_eval_unary_float_in_norm(arg, out, zeta_numeric_eval) {
return;
}
if let Ok(rat) = Rational::try_from(arg)
&& let Some(exact) = zeta_exact_rational(&rat)
{
**out = exact;
}
}
})
.with_derivative_function(move |x, _, out| {
if let Some([arg]) = function_arguments::<1>(x)
{
**out = function!(dzeta, 1, arg);
}
})
.with_series_function(move |args| {
let [arg] = args else { return None; };
let point = arg.coefficient(0.into());
if point != Coefficient::one() {
return None;
}
let arg = arg.to_atom();
let pole = arg.to_owned() - 1;
let regularized = (Symbol::PI.to_atom() * 2).pow(&pole) * -2 * (pole.clone() * Symbol::PI / 2).cos() * function!(gamma, -&pole + 1) * function!(zeta_symbol, -&pole);
Some((Atom::num(1) / pole, regularized))
})
.with_evaluation_info(
EvaluationInfo::new()
.register(|args: &[Complex<Float>]| {
unary_eval_complex_float(args, zeta_numeric_eval)
})
.register(|args: &[Float]| {
let prec = args.first().map(|x| x.prec()).unwrap_or(53);
let [arg] = args else {
return Float::with_val(53, f64::NAN);
};
zeta_numeric_eval(&Complex::new(arg.clone(), Float::new(prec)), prec).re
})
.register(|args: &[f64]| zeta_eval_f64(args))
.register(|args: &[Complex<f64>]| zeta_eval_complex_f64(args)),
)
}
,
"erf";;
|symbols, b| {
let erf_symbol = symbols[5];
b.with_normalization_function(move |x, out| {
if let Some([arg]) = function_arguments::<1>(x) {
if arg.is_zero() {
out.to_num(Coefficient::zero());
return;
}
if maybe_eval_unary_float_in_norm(arg, out, erf_numeric_eval) {
return;
}
if is_negative_atom(arg) {
**out = -function!(erf_symbol, -arg.to_owned());
}
}
})
.with_derivative_function(move |x, i, out| {
if i == 0
&& let Some([arg]) = function_arguments::<1>(x)
{
**out = Atom::num(2)
* (-arg.to_owned().pow(Atom::num(2))).exp()
/ Atom::from(State::PI).sqrt();
}
})
.with_evaluation_info(
EvaluationInfo::new()
.register(|args: &[Complex<Float>]| {
unary_eval_complex_float(args, erf_numeric_eval)
})
.register(|args: &[Float]| {
let prec = args.first().map(|x| x.prec()).unwrap_or(53);
let [arg] = args else {
return Float::with_val(53, f64::NAN);
};
erf_numeric_eval(&Complex::new(arg.clone(), Float::new(prec)), prec).re
})
.register(|args: &[f64]| unary_eval_f64(args, erf_numeric_eval))
.register(|args: &[Complex<f64>]| {
unary_eval_complex_f64(args, erf_numeric_eval)
}),
)
}
);
let erf = symbols.pop().unwrap();
let zeta = symbols.pop().unwrap();
let _dzeta = symbols.pop().unwrap();
let polylog = symbols.pop().unwrap();
let polygamma = symbols.pop().unwrap();
let gamma = symbols.pop().unwrap();
Self {
euler_gamma,
gamma,
erf,
polygamma,
polylog,
zeta,
}
}
}
impl GeometricSymbols {
fn new() -> Self {
let mut symbols = crate::symbol_group!(
"tan";;
|symbols, b| {
let tan = symbols[0];
let sec = symbols[2];
b.with_normalization_function(|x, out| {
if let Some([arg]) = function_arguments::<1>(x) {
if arg.is_zero() {
out.to_num(Coefficient::zero());
return;
}
let _ = maybe_eval_unary_float_in_norm(arg, out, |z, _| z.tan());
}
})
.with_derivative_function(move |x, i, out| {
if i == 0 && let Some([arg]) = function_arguments::<1>(x) {
**out = function!(sec, arg).pow(Atom::num(2));
}
})
.with_series_function(move |args| {
let [arg] = args else { return None; };
let point = arg.coefficient(0.into());
if !is_tan_pole(point.as_view()) {
return None;
}
let delta = arg.to_atom() - &point;
Some((Atom::num(1) / &delta, -&delta / function!(tan, delta)))
})
.with_evaluation_info(
EvaluationInfo::new().register(|args: &[Complex<Float>]| {
unary_eval_complex_float(args, |z, _| z.tan())
})
.register(|args: &[Float]| unary_eval_real(args, |x| x.tan()))
.register(|args: &[f64]| unary_eval_real_f64(args, f64::tan))
.register(|args: &[Complex<f64>]| unary_eval_complex_real_f64(args, |z| z.tan()))
.register(|args: &[ErrorPropagatingFloat<f64>]| {
unary_eval_real(args, |x| x.tan())
})
.register(|args: &[ErrorPropagatingFloat<Float>]| {
unary_eval_real(args, |x| x.tan())
}),
)
},
"cot";;
|symbols, b| {
let tan = symbols[0];
let csc = symbols[3];
b.with_normalization_function(|x, out| {
if let Some([arg]) = function_arguments::<1>(x) {
if arg.is_zero() {
out.to_num(Coefficient::complex_infinity());
return;
}
let _ = maybe_eval_unary_float_in_norm(arg, out, |z, _| z.tan().inv());
}
})
.with_derivative_function(move |x, i, out| {
if i == 0 && let Some([arg]) = function_arguments::<1>(x) {
**out = -function!(csc, arg).pow(Atom::num(2));
}
})
.with_series_function(move |args| {
let [arg] = args else { return None; };
let point = arg.coefficient(0.into());
if !is_cot_csc_pole(point.as_view()) {
return None;
}
let delta = arg.to_atom() - &point;
Some((Atom::num(1) / &delta, &delta / function!(tan, delta)))
})
.with_evaluation_info(
EvaluationInfo::new().register(|args: &[Complex<Float>]| {
unary_eval_complex_float(args, |z, _| z.tan().inv())
})
.register(|args: &[Float]| unary_eval_real(args, |x| x.tan().inv()))
.register(|args: &[f64]| unary_eval_real_f64(args, |x| x.tan().recip()))
.register(|args: &[Complex<f64>]| unary_eval_complex_real_f64(args, |z| z.tan().inv()))
.register(|args: &[ErrorPropagatingFloat<f64>]| {
unary_eval_real(args, |x| x.tan().inv())
})
.register(|args: &[ErrorPropagatingFloat<Float>]| {
unary_eval_real(args, |x| x.tan().inv())
}),
)
},
"sec";;
|symbols, b| {
let tan = symbols[0];
let sec = symbols[2];
b.with_normalization_function(|x, out| {
if let Some([arg]) = function_arguments::<1>(x) {
if arg.is_zero() {
out.to_num(Coefficient::one());
return;
}
let _ = maybe_eval_unary_float_in_norm(arg, out, |z, _| z.cos().inv());
}
})
.with_derivative_function(move |x, i, out| {
if i == 0 && let Some([arg]) = function_arguments::<1>(x) {
**out = function!(sec, arg) * function!(tan, arg);
}
})
.with_series_function(move |args| {
let [arg] = args else { return None; };
let point = arg.coefficient(0.into());
let residue = sec_residue(point.as_view())?;
let delta = arg.to_atom() - &point;
Some((
Atom::num(1) / &delta,
residue * &delta / function!(State::SIN, delta),
))
})
.with_evaluation_info(
EvaluationInfo::new().register(|args: &[Complex<Float>]| {
unary_eval_complex_float(args, |z, _| z.cos().inv())
})
.register(|args: &[Float]| unary_eval_real(args, |x| x.cos().inv()))
.register(|args: &[f64]| unary_eval_real_f64(args, |x| x.cos().recip()))
.register(|args: &[Complex<f64>]| unary_eval_complex_real_f64(args, |z| z.cos().inv()))
.register(|args: &[ErrorPropagatingFloat<f64>]| {
unary_eval_real(args, |x| x.cos().inv())
})
.register(|args: &[ErrorPropagatingFloat<Float>]| {
unary_eval_real(args, |x| x.cos().inv())
}),
)
},
"csc";;
|symbols, b| {
let cot = symbols[1];
let csc = symbols[3];
b.with_normalization_function(|x, out| {
if let Some([arg]) = function_arguments::<1>(x) {
if arg.is_zero() {
out.to_num(Coefficient::complex_infinity());
return;
}
let _ = maybe_eval_unary_float_in_norm(arg, out, |z, _| z.sin().inv());
}
})
.with_derivative_function(move |x, i, out| {
if i == 0 && let Some([arg]) = function_arguments::<1>(x) {
**out = -function!(csc, arg) * function!(cot, arg);
}
})
.with_series_function(move |args| {
let [arg] = args else { return None; };
let point = arg.coefficient(0.into());
let residue = csc_residue(point.as_view())?;
let delta = arg.to_atom() - &point;
Some((
Atom::num(1) / &delta,
residue * &delta / function!(State::SIN, delta),
))
})
.with_evaluation_info(
EvaluationInfo::new().register(|args: &[Complex<Float>]| {
unary_eval_complex_float(args, |z, _| z.sin().inv())
})
.register(|args: &[Float]| unary_eval_real(args, |x| x.sin().inv()))
.register(|args: &[f64]| unary_eval_real_f64(args, |x| x.sin().recip()))
.register(|args: &[Complex<f64>]| unary_eval_complex_real_f64(args, |z| z.sin().inv()))
.register(|args: &[ErrorPropagatingFloat<f64>]| {
unary_eval_real(args, |x| x.sin().inv())
})
.register(|args: &[ErrorPropagatingFloat<Float>]| {
unary_eval_real(args, |x| x.sin().inv())
}),
)
},
"sinh";;
|symbols, b| {
let cosh = symbols[5];
b.with_normalization_function(|x, out| {
if let Some([arg]) = function_arguments::<1>(x) {
if arg.is_zero() {
out.to_num(Coefficient::zero());
return;
}
let _ = maybe_eval_unary_float_in_norm(arg, out, |z, _| z.sinh());
}
})
.with_derivative_function(move |x, i, out| {
if i == 0 && let Some([arg]) = function_arguments::<1>(x) {
**out = function!(cosh, arg);
}
})
.with_evaluation_info(
EvaluationInfo::new().register(|args: &[Complex<Float>]| {
unary_eval_complex_float(args, |z, _| z.sinh())
})
.register(|args: &[Float]| unary_eval_real(args, |x| x.sinh()))
.register(|args: &[f64]| unary_eval_real_f64(args, f64::sinh))
.register(|args: &[Complex<f64>]| unary_eval_complex_real_f64(args, |z| z.sinh()))
.register(|args: &[ErrorPropagatingFloat<f64>]| {
unary_eval_real(args, |x| x.sinh())
})
.register(|args: &[ErrorPropagatingFloat<Float>]| {
unary_eval_real(args, |x| x.sinh())
}),
)
},
"cosh";;
|symbols, b| {
let sinh = symbols[4];
b.with_normalization_function(|x, out| {
if let Some([arg]) = function_arguments::<1>(x) {
if arg.is_zero() {
out.to_num(Coefficient::one());
return;
}
let _ = maybe_eval_unary_float_in_norm(arg, out, |z, _| z.cosh());
}
})
.with_derivative_function(move |x, i, out| {
if i == 0 && let Some([arg]) = function_arguments::<1>(x) {
**out = function!(sinh, arg);
}
})
.with_evaluation_info(
EvaluationInfo::new().register(|args: &[Complex<Float>]| {
unary_eval_complex_float(args, |z, _| z.cosh())
})
.register(|args: &[Float]| unary_eval_real(args, |x| x.cosh()))
.register(|args: &[f64]| unary_eval_real_f64(args, f64::cosh))
.register(|args: &[Complex<f64>]| unary_eval_complex_real_f64(args, |z| z.cosh()))
.register(|args: &[ErrorPropagatingFloat<f64>]| {
unary_eval_real(args, |x| x.cosh())
})
.register(|args: &[ErrorPropagatingFloat<Float>]| {
unary_eval_real(args, |x| x.cosh())
}),
)
},
"tanh";;
|symbols, b| {
let tanh = symbols[6];
let sech = symbols[8];
b.with_normalization_function(|x, out| {
if let Some([arg]) = function_arguments::<1>(x) {
if arg.is_zero() {
out.to_num(Coefficient::zero());
return;
}
let _ = maybe_eval_unary_float_in_norm(arg, out, |z, _| z.tanh());
}
})
.with_derivative_function(move |x, i, out| {
if i == 0 && let Some([arg]) = function_arguments::<1>(x) {
**out = function!(sech, arg).pow(Atom::num(2));
}
})
.with_series_function(move |args| {
let [arg] = args else { return None; };
let point = arg.coefficient(0.into());
if !is_tanh_sech_pole(point.as_view()) {
return None;
}
let delta = arg.to_atom() - &point;
Some((Atom::num(1) / &delta, &delta / function!(tanh, delta)))
})
.with_evaluation_info(
EvaluationInfo::new().register(|args: &[Complex<Float>]| {
unary_eval_complex_float(args, |z, _| z.tanh())
})
.register(|args: &[Float]| unary_eval_real(args, |x| x.tanh()))
.register(|args: &[f64]| unary_eval_real_f64(args, f64::tanh))
.register(|args: &[Complex<f64>]| unary_eval_complex_real_f64(args, |z| z.tanh()))
.register(|args: &[ErrorPropagatingFloat<f64>]| {
unary_eval_real(args, |x| x.tanh())
})
.register(|args: &[ErrorPropagatingFloat<Float>]| {
unary_eval_real(args, |x| x.tanh())
}),
)
},
"coth";;
|symbols, b| {
let tanh = symbols[6];
let csch = symbols[9];
b.with_normalization_function(|x, out| {
if let Some([arg]) = function_arguments::<1>(x) {
if arg.is_zero() {
out.to_num(Coefficient::complex_infinity());
return;
}
let _ = maybe_eval_unary_float_in_norm(arg, out, |z, _| z.tanh().inv());
}
})
.with_derivative_function(move |x, i, out| {
if i == 0 && let Some([arg]) = function_arguments::<1>(x) {
**out = -function!(csch, arg).pow(Atom::num(2));
}
})
.with_series_function(move |args| {
let [arg] = args else { return None; };
let point = arg.coefficient(0.into());
if !is_coth_csch_pole(point.as_view()) {
return None;
}
let delta = arg.to_atom() - &point;
Some((Atom::num(1) / &delta, &delta / function!(tanh, delta)))
})
.with_evaluation_info(
EvaluationInfo::new().register(|args: &[Complex<Float>]| {
unary_eval_complex_float(args, |z, _| z.tanh().inv())
})
.register(|args: &[Float]| unary_eval_real(args, |x| x.tanh().inv()))
.register(|args: &[f64]| unary_eval_real_f64(args, |x| x.tanh().recip()))
.register(|args: &[Complex<f64>]| unary_eval_complex_real_f64(args, |z| z.tanh().inv()))
.register(|args: &[ErrorPropagatingFloat<f64>]| {
unary_eval_real(args, |x| x.tanh().inv())
})
.register(|args: &[ErrorPropagatingFloat<Float>]| {
unary_eval_real(args, |x| x.tanh().inv())
}),
)
},
"sech";;
|symbols, b| {
let sinh = symbols[4];
let sech = symbols[8];
let tanh = symbols[6];
b.with_normalization_function(|x, out| {
if let Some([arg]) = function_arguments::<1>(x) {
if arg.is_zero() {
out.to_num(Coefficient::one());
return;
}
let _ = maybe_eval_unary_float_in_norm(arg, out, |z, _| z.cosh().inv());
}
})
.with_derivative_function(move |x, i, out| {
if i == 0 && let Some([arg]) = function_arguments::<1>(x) {
**out = -function!(sech, arg) * function!(tanh, arg);
}
})
.with_series_function(move |args| {
let [arg] = args else { return None; };
let point = arg.coefficient(0.into());
let residue = sech_residue(point.as_view())?;
let delta = arg.to_atom() - &point;
Some((Atom::num(1) / &delta, residue * &delta / function!(sinh, delta)))
})
.with_evaluation_info(
EvaluationInfo::new().register(|args: &[Complex<Float>]| {
unary_eval_complex_float(args, |z, _| z.cosh().inv())
})
.register(|args: &[Float]| unary_eval_real(args, |x| x.cosh().inv()))
.register(|args: &[f64]| unary_eval_real_f64(args, |x| x.cosh().recip()))
.register(|args: &[Complex<f64>]| unary_eval_complex_real_f64(args, |z| z.cosh().inv()))
.register(|args: &[ErrorPropagatingFloat<f64>]| {
unary_eval_real(args, |x| x.cosh().inv())
})
.register(|args: &[ErrorPropagatingFloat<Float>]| {
unary_eval_real(args, |x| x.cosh().inv())
}),
)
},
"csch";;
|symbols, b| {
let sinh = symbols[4];
let csch = symbols[9];
let coth = symbols[7];
b.with_normalization_function(|x, out| {
if let Some([arg]) = function_arguments::<1>(x) {
if arg.is_zero() {
out.to_num(Coefficient::complex_infinity());
return;
}
let _ = maybe_eval_unary_float_in_norm(arg, out, |z, _| z.sinh().inv());
}
})
.with_derivative_function(move |x, i, out| {
if i == 0 && let Some([arg]) = function_arguments::<1>(x) {
**out = -function!(csch, arg) * function!(coth, arg);
}
})
.with_series_function(move |args| {
let [arg] = args else { return None; };
let point = arg.coefficient(0.into());
let residue = csch_residue(point.as_view())?;
let delta = arg.to_atom() - &point;
Some((Atom::num(1) / &delta, residue * &delta / function!(sinh, delta)))
})
.with_evaluation_info(
EvaluationInfo::new().register(|args: &[Complex<Float>]| {
unary_eval_complex_float(args, |z, _| z.sinh().inv())
})
.register(|args: &[Float]| unary_eval_real(args, |x| x.sinh().inv()))
.register(|args: &[f64]| unary_eval_real_f64(args, |x| x.sinh().recip()))
.register(|args: &[Complex<f64>]| unary_eval_complex_real_f64(args, |z| z.sinh().inv()))
.register(|args: &[ErrorPropagatingFloat<f64>]| {
unary_eval_real(args, |x| x.sinh().inv())
})
.register(|args: &[ErrorPropagatingFloat<Float>]| {
unary_eval_real(args, |x| x.sinh().inv())
}),
)
}
);
let csch = symbols.pop().unwrap();
let sech = symbols.pop().unwrap();
let coth = symbols.pop().unwrap();
let tanh = symbols.pop().unwrap();
let cosh = symbols.pop().unwrap();
let sinh = symbols.pop().unwrap();
let csc = symbols.pop().unwrap();
let sec = symbols.pop().unwrap();
let cot = symbols.pop().unwrap();
let tan = symbols.pop().unwrap();
let asin = symbol!(
"asin",
norm = |x, out| {
if let Some([arg]) = function_arguments::<1>(x) {
if arg.is_zero() {
out.to_num(Coefficient::zero());
return;
}
let _ = maybe_eval_unary_float_in_norm(arg, out, |z, _| z.asin());
}
},
der = move |x, i, out| {
if i == 0
&& let Some([arg]) = function_arguments::<1>(x)
{
let arg = arg.to_owned();
**out =
Atom::num(1) / function!(State::SQRT, Atom::num(1) - arg.pow(Atom::num(2)));
}
},
eval = EvaluationInfo::new()
.register(|args: &[Complex<Float>]| {
unary_eval_complex_float(args, |z, _| z.asin())
})
.register(|args: &[Float]| unary_eval_real(args, |x| x.asin()))
.register(|args: &[f64]| unary_eval_real_f64(args, f64::asin))
.register(|args: &[Complex<f64>]| {
unary_eval_complex_real_f64(args, |z| z.asin())
})
.register(|args: &[ErrorPropagatingFloat<f64>]| {
unary_eval_real(args, |x| x.asin())
})
.register(|args: &[ErrorPropagatingFloat<Float>]| {
unary_eval_real(args, |x| x.asin())
})
);
let acos = symbol!(
"acos",
norm = |x, out| {
if let Some([arg]) = function_arguments::<1>(x) {
if arg.is_one() {
out.to_num(Coefficient::zero());
return;
}
let _ = maybe_eval_unary_float_in_norm(arg, out, |z, _| z.acos());
}
},
der = move |x, i, out| {
if i == 0
&& let Some([arg]) = function_arguments::<1>(x)
{
let arg = arg.to_owned();
**out = -Atom::num(1)
/ function!(State::SQRT, Atom::num(1) - arg.pow(Atom::num(2)));
}
},
eval = EvaluationInfo::new()
.register(|args: &[Complex<Float>]| {
unary_eval_complex_float(args, |z, _| z.acos())
})
.register(|args: &[Float]| unary_eval_real(args, |x| x.acos()))
.register(|args: &[f64]| unary_eval_real_f64(args, f64::acos))
.register(|args: &[Complex<f64>]| {
unary_eval_complex_real_f64(args, |z| z.acos())
})
.register(|args: &[ErrorPropagatingFloat<f64>]| {
unary_eval_real(args, |x| x.acos())
})
.register(|args: &[ErrorPropagatingFloat<Float>]| {
unary_eval_real(args, |x| x.acos())
})
);
let atan = symbol!(
"atan",
norm = |x, out| {
if let Some([arg]) = function_arguments::<1>(x) {
if arg.is_zero() {
out.to_num(Coefficient::zero());
return;
}
let _ = maybe_eval_unary_float_in_norm(arg, out, atan_numeric_eval);
}
if let Some([x, y]) = function_arguments::<2>(x)
&& maybe_eval_binary_float_in_norm(x, y, out, |x, y, prec| {
Some(atan2_numeric_eval(x, y, prec))
})
{}
},
der = move |x, i, out| {
if i == 0
&& let Some([arg]) = function_arguments::<1>(x)
{
let arg = arg.to_owned();
**out = Atom::num(1) / (Atom::num(1) + arg.pow(Atom::num(2)));
}
if let Some([x, y]) = function_arguments::<2>(x) {
let x = x.to_owned();
let y = y.to_owned();
let denom = x.clone().pow(Atom::num(2)) + y.clone().pow(Atom::num(2));
if i == 0 {
**out = -y / denom;
} else if i == 1 {
**out = x / denom;
}
}
},
eval = EvaluationInfo::new()
.register(|args: &[Complex<Float>]| match args {
[z] => atan_numeric_eval(z, z.re.prec().max(z.im.prec())),
[x, y] => atan2_numeric_eval(
x,
y,
y.re.prec()
.max(y.im.prec())
.max(x.re.prec().max(x.im.prec())),
),
_ =>
Complex::new(Float::with_val(53, f64::NAN), Float::with_val(53, f64::NAN),),
})
.register(|args: &[Float]| atan_eval_real(args))
.register(|args: &[f64]| match args {
[z] => z.atan(),
[x, y] => (*y).atan2(*x),
_ => f64::NAN,
})
.register(|args: &[Complex<f64>]| match args {
[z] => atan_eval_complex_f64(*z),
[x, y] => {
if y.im == 0.0 && x.im == 0.0 {
Complex::new(y.re.atan2(x.re), 0.0)
} else {
y.atan2(x)
}
}
_ => Complex::new(f64::NAN, f64::NAN),
})
.register(|args: &[ErrorPropagatingFloat<f64>]| atan_eval_real(args))
.register(|args: &[ErrorPropagatingFloat<Float>]| atan_eval_real(args))
);
let acot = symbol!(
"acot",
norm = |x, out| {
if let Some([arg]) = function_arguments::<1>(x) {
let _ = maybe_eval_unary_float_in_norm(arg, out, |z, prec| {
atan_numeric_eval(&z.inv(), prec)
});
}
},
der = move |x, i, out| {
if i == 0
&& let Some([arg]) = function_arguments::<1>(x)
{
let arg = arg.to_owned();
**out = -Atom::num(1) / (Atom::num(1) + arg.pow(Atom::num(2)));
}
},
eval = EvaluationInfo::new()
.register(|args: &[Complex<Float>]| {
unary_eval_complex_float(args, |z, prec| atan_numeric_eval(&z.inv(), prec))
})
.register(|args: &[Float]| {
unary_eval_real(args, |x| {
let inv = x.inv();
inv.atan2(&inv.one())
})
})
.register(|args: &[f64]| unary_eval_real_f64(args, |x| x.recip().atan()))
.register(|args: &[Complex<f64>]| {
unary_eval_complex_real_f64(args, |z| atan_eval_complex_f64(z.inv()))
})
.register(|args: &[ErrorPropagatingFloat<f64>]| {
unary_eval_real(args, |x| {
let inv = x.inv();
inv.atan2(&inv.one())
})
})
.register(|args: &[ErrorPropagatingFloat<Float>]| {
unary_eval_real(args, |x| {
let inv = x.inv();
inv.atan2(&inv.one())
})
})
);
let asec = symbol!(
"asec",
norm = |x, out| {
if let Some([arg]) = function_arguments::<1>(x) {
if arg.is_one() {
out.to_num(Coefficient::zero());
return;
}
let _ = maybe_eval_unary_float_in_norm(arg, out, |z, _| z.inv().acos());
}
},
der = move |x, i, out| {
if i == 0
&& let Some([arg]) = function_arguments::<1>(x)
{
let inv = Atom::num(1) / arg;
**out = Atom::num(1)
/ (arg.pow(Atom::num(2))
* function!(State::SQRT, Atom::num(1) - inv.pow(Atom::num(2))));
}
},
eval = EvaluationInfo::new()
.register(|args: &[Complex<Float>]| {
unary_eval_complex_float(args, |z, _| z.inv().acos())
})
.register(|args: &[Float]| unary_eval_real(args, |x| x.inv().acos()))
.register(|args: &[f64]| unary_eval_real_f64(args, |x| x.recip().acos()))
.register(|args: &[Complex<f64>]| {
unary_eval_complex_real_f64(args, |z| z.inv().acos())
})
.register(|args: &[ErrorPropagatingFloat<f64>]| {
unary_eval_real(args, |x| x.inv().acos())
})
.register(|args: &[ErrorPropagatingFloat<Float>]| {
unary_eval_real(args, |x| x.inv().acos())
})
);
let acsc = symbol!(
"acsc",
norm = |x, out| {
if let Some([arg]) = function_arguments::<1>(x) {
let _ = maybe_eval_unary_float_in_norm(arg, out, |z, _| z.inv().asin());
}
},
der = move |x, i, out| {
if i == 0
&& let Some([arg]) = function_arguments::<1>(x)
{
let inv = Atom::num(1) / arg;
**out = -Atom::num(1)
/ (arg.pow(Atom::num(2))
* function!(State::SQRT, Atom::num(1) - inv.pow(Atom::num(2))));
}
},
eval = EvaluationInfo::new()
.register(|args: &[Complex<Float>]| {
unary_eval_complex_float(args, |z, _| z.inv().asin())
})
.register(|args: &[Float]| unary_eval_real(args, |x| x.inv().asin()))
.register(|args: &[f64]| unary_eval_real_f64(args, |x| x.recip().asin()))
.register(|args: &[Complex<f64>]| {
unary_eval_complex_real_f64(args, |z| z.inv().asin())
})
.register(|args: &[ErrorPropagatingFloat<f64>]| {
unary_eval_real(args, |x| x.inv().asin())
})
.register(|args: &[ErrorPropagatingFloat<Float>]| {
unary_eval_real(args, |x| x.inv().asin())
})
);
let asinh = symbol!(
"asinh",
norm = |x, out| {
if let Some([arg]) = function_arguments::<1>(x) {
if arg.is_zero() {
out.to_num(Coefficient::zero());
return;
}
let _ = maybe_eval_unary_float_in_norm(arg, out, |z, _| z.asinh());
}
},
der = move |x, i, out| {
if i == 0
&& let Some([arg]) = function_arguments::<1>(x)
{
let arg = arg.to_owned();
**out =
Atom::num(1) / function!(State::SQRT, Atom::num(1) + arg.pow(Atom::num(2)));
}
},
eval = EvaluationInfo::new()
.register(|args: &[Complex<Float>]| {
unary_eval_complex_float(args, |z, _| z.asinh())
})
.register(|args: &[Float]| unary_eval_real(args, |x| x.asinh()))
.register(|args: &[f64]| unary_eval_real_f64(args, f64::asinh))
.register(|args: &[Complex<f64>]| {
unary_eval_complex_real_f64(args, |z| z.asinh())
})
.register(|args: &[ErrorPropagatingFloat<f64>]| {
unary_eval_real(args, |x| x.asinh())
})
.register(|args: &[ErrorPropagatingFloat<Float>]| {
unary_eval_real(args, |x| x.asinh())
})
);
let acosh = symbol!(
"acosh",
norm = |x, out| {
if let Some([arg]) = function_arguments::<1>(x) {
if arg.is_one() {
out.to_num(Coefficient::zero());
return;
}
let _ = maybe_eval_unary_float_in_norm(arg, out, |z, _| z.acosh());
}
},
der = move |x, i, out| {
if i == 0
&& let Some([arg]) = function_arguments::<1>(x)
{
let arg = arg.to_owned();
**out = Atom::num(1)
/ (function!(State::SQRT, &arg - Atom::num(1))
* function!(State::SQRT, arg + Atom::num(1)));
}
},
eval = EvaluationInfo::new()
.register(|args: &[Complex<Float>]| {
unary_eval_complex_float(args, |z, _| z.acosh())
})
.register(|args: &[Float]| unary_eval_real(args, |x| x.acosh()))
.register(|args: &[f64]| unary_eval_real_f64(args, f64::acosh))
.register(|args: &[Complex<f64>]| {
unary_eval_complex_real_f64(args, |z| z.acosh())
})
.register(|args: &[ErrorPropagatingFloat<f64>]| {
unary_eval_real(args, |x| x.acosh())
})
.register(|args: &[ErrorPropagatingFloat<Float>]| {
unary_eval_real(args, |x| x.acosh())
})
);
let atanh = symbol!(
"atanh",
norm = |x, out| {
if let Some([arg]) = function_arguments::<1>(x) {
if arg.is_zero() {
out.to_num(Coefficient::zero());
return;
}
let _ = maybe_eval_unary_float_in_norm(arg, out, |z, _| z.atanh());
}
},
der = move |x, i, out| {
if i == 0
&& let Some([arg]) = function_arguments::<1>(x)
{
let arg = arg.to_owned();
**out = Atom::num(1) / (Atom::num(1) - arg.pow(Atom::num(2)));
}
},
eval = EvaluationInfo::new()
.register(|args: &[Complex<Float>]| {
unary_eval_complex_float(args, |z, _| z.atanh())
})
.register(|args: &[Float]| unary_eval_real(args, |x| x.atanh()))
.register(|args: &[f64]| unary_eval_real_f64(args, f64::atanh))
.register(|args: &[Complex<f64>]| {
unary_eval_complex_real_f64(args, |z| z.atanh())
})
.register(|args: &[ErrorPropagatingFloat<f64>]| {
unary_eval_real(args, |x| x.atanh())
})
.register(|args: &[ErrorPropagatingFloat<Float>]| {
unary_eval_real(args, |x| x.atanh())
})
);
let acoth = symbol!(
"acoth",
norm = |x, out| {
if let Some([arg]) = function_arguments::<1>(x) {
let _ = maybe_eval_unary_float_in_norm(arg, out, |z, _| z.inv().atanh());
}
},
der = move |x, i, out| {
if i == 0
&& let Some([arg]) = function_arguments::<1>(x)
{
let arg = arg.to_owned();
**out = Atom::num(1) / (Atom::num(1) - arg.pow(Atom::num(2)));
}
},
eval = EvaluationInfo::new()
.register(|args: &[Complex<Float>]| {
unary_eval_complex_float(args, |z, _| z.inv().atanh())
})
.register(|args: &[Float]| unary_eval_real(args, |x| x.inv().atanh()))
.register(|args: &[f64]| unary_eval_real_f64(args, |x| x.recip().atanh()))
.register(|args: &[Complex<f64>]| {
unary_eval_complex_real_f64(args, |z| z.inv().atanh())
})
.register(|args: &[ErrorPropagatingFloat<f64>]| {
unary_eval_real(args, |x| x.inv().atanh())
})
.register(|args: &[ErrorPropagatingFloat<Float>]| {
unary_eval_real(args, |x| x.inv().atanh())
})
);
let asech = symbol!(
"asech",
norm = |x, out| {
if let Some([arg]) = function_arguments::<1>(x) {
if arg.is_one() {
out.to_num(Coefficient::zero());
return;
}
let _ = maybe_eval_unary_float_in_norm(arg, out, |z, _| z.inv().acosh());
}
},
der = move |x, i, out| {
if i == 0
&& let Some([arg]) = function_arguments::<1>(x)
{
let arg = arg.to_owned();
let inv = Atom::num(1) / &arg;
**out = -Atom::num(1)
/ (arg.pow(Atom::num(2))
* function!(State::SQRT, &inv - Atom::num(1))
* function!(State::SQRT, inv + Atom::num(1)));
}
},
eval = EvaluationInfo::new()
.register(|args: &[Complex<Float>]| {
unary_eval_complex_float(args, |z, _| z.inv().acosh())
})
.register(|args: &[Float]| unary_eval_real(args, |x| x.inv().acosh()))
.register(|args: &[f64]| unary_eval_real_f64(args, |x| x.recip().acosh()))
.register(|args: &[Complex<f64>]| {
unary_eval_complex_real_f64(args, |z| z.inv().acosh())
})
.register(|args: &[ErrorPropagatingFloat<f64>]| {
unary_eval_real(args, |x| x.inv().acosh())
})
.register(|args: &[ErrorPropagatingFloat<Float>]| {
unary_eval_real(args, |x| x.inv().acosh())
})
);
let acsch = symbol!(
"acsch",
norm = |x, out| {
if let Some([arg]) = function_arguments::<1>(x) {
let _ = maybe_eval_unary_float_in_norm(arg, out, |z, _| z.inv().asinh());
}
},
der = move |x, i, out| {
if i == 0
&& let Some([arg]) = function_arguments::<1>(x)
{
let arg = arg.to_owned();
let inv = Atom::num(1) / &arg;
**out = -Atom::num(1)
/ (arg.pow(Atom::num(2))
* function!(State::SQRT, Atom::num(1) + inv.pow(Atom::num(2))));
}
},
eval = EvaluationInfo::new()
.register(|args: &[Complex<Float>]| {
unary_eval_complex_float(args, |z, _| z.inv().asinh())
})
.register(|args: &[Float]| unary_eval_real(args, |x| x.inv().asinh()))
.register(|args: &[f64]| unary_eval_real_f64(args, |x| x.recip().asinh()))
.register(|args: &[Complex<f64>]| {
unary_eval_complex_real_f64(args, |z| z.inv().asinh())
})
.register(|args: &[ErrorPropagatingFloat<f64>]| {
unary_eval_real(args, |x| x.inv().asinh())
})
.register(|args: &[ErrorPropagatingFloat<Float>]| {
unary_eval_real(args, |x| x.inv().asinh())
})
);
Self {
tan,
cot,
sec,
csc,
asin,
acos,
atan,
acot,
asec,
acsc,
sinh,
cosh,
tanh,
coth,
sech,
csch,
asinh,
acosh,
atanh,
acoth,
asech,
acsch,
}
}
}
impl BesselSymbols {
fn new() -> Self {
let bessel_j = symbol!(
"bessel_j",
norm = |x, out| {
if let Some([nu, z]) = function_arguments::<2>(x) {
if z.is_zero()
&& let Some(n) = atom_to_integer(nu)
{
if n == 0 {
out.to_num(Coefficient::one());
} else {
out.to_num(Coefficient::zero());
}
return;
}
let _ = maybe_eval_binary_float_in_norm(nu, z, out, bessel_j_numeric_eval);
}
},
der = move |x, i, out| {
if i == 1
&& let Some([nu, z]) = function_arguments::<2>(x)
{
let symbol = x.as_fun_view().unwrap().get_symbol();
let nu = nu.to_owned();
**out = (function!(symbol, nu.clone() - Atom::num(1), z)
- function!(symbol, nu + Atom::num(1), z))
/ Atom::num(2);
}
},
eval = EvaluationInfo::new()
.with_tags(1)
.register_tagged(|tags| {
let tags = tags.iter().map(|x| x.to_owned()).collect::<Vec<_>>();
Box::new(move |args: &[Complex<Float>]| {
tagged_unary_eval_complex_float(
"bessel_j",
&tags,
args,
bessel_j_numeric_eval,
)
})
})
.register_tagged(|tags| {
let tags = tags.iter().map(|x| x.to_owned()).collect::<Vec<_>>();
Box::new(move |args: &[Float]| {
let prec = args.first().map(|x| x.prec()).unwrap_or(53);
let tag_views = tags.iter().map(|x| x.as_view()).collect::<Vec<_>>();
let Ok(order) = complex_float_tag("bessel_j", &tag_views, prec) else {
return Float::with_val(53, f64::NAN);
};
let [arg] = args else {
return Float::with_val(53, f64::NAN);
};
bessel_j_numeric_eval(
&order,
&Complex::new(arg.clone(), Float::new(prec)),
prec,
)
.map(|z| z.re)
.unwrap_or_else(|| Float::with_val(53, f64::NAN))
})
})
.register_tagged(|tags| {
let order = complex_float_tag("bessel_j", tags, 53).ok();
Box::new(move |args: &[f64]| {
tagged_unary_eval_real_f64(order.clone(), args, bessel_j_numeric_eval)
})
})
.register_tagged(|tags| {
let order = complex_float_tag("bessel_j", tags, 53).ok();
Box::new(move |args: &[Complex<f64>]| {
tagged_unary_eval_complex_f64(order.clone(), args, bessel_j_numeric_eval)
})
})
);
let bessel_y = symbol!(
"bessel_y",
norm = |x, out| {
if let Some([nu, z]) = function_arguments::<2>(x) {
if z.is_zero() {
out.to_num(Coefficient::complex_infinity());
return;
}
let _ = maybe_eval_binary_float_in_norm(nu, z, out, bessel_y_numeric_eval);
}
},
der = move |x, i, out| {
if i == 1
&& let Some([nu, z]) = function_arguments::<2>(x)
{
let symbol = x.as_fun_view().unwrap().get_symbol();
let nu = nu.to_owned();
**out = (function!(symbol, nu.clone() - Atom::num(1), z)
- function!(symbol, nu + Atom::num(1), z))
/ Atom::num(2);
}
},
eval = EvaluationInfo::new()
.with_tags(1)
.register_tagged(|tags| {
let tags = tags.iter().map(|x| x.to_owned()).collect::<Vec<_>>();
Box::new(move |args: &[Complex<Float>]| {
tagged_unary_eval_complex_float(
"bessel_y",
&tags,
args,
bessel_y_numeric_eval,
)
})
})
.register_tagged(|tags| {
let tags = tags.iter().map(|x| x.to_owned()).collect::<Vec<_>>();
Box::new(move |args: &[Float]| {
let prec = args.first().map(|x| x.prec()).unwrap_or(53);
let tag_views = tags.iter().map(|x| x.as_view()).collect::<Vec<_>>();
let Ok(order) = complex_float_tag("bessel_y", &tag_views, prec) else {
return Float::with_val(53, f64::NAN);
};
let [arg] = args else {
return Float::with_val(53, f64::NAN);
};
bessel_y_numeric_eval(
&order,
&Complex::new(arg.clone(), Float::new(prec)),
prec,
)
.map(|z| z.re)
.unwrap_or_else(|| Float::with_val(53, f64::NAN))
})
})
.register_tagged(|tags| {
let order = complex_float_tag("bessel_y", tags, 53).ok();
Box::new(move |args: &[f64]| {
tagged_unary_eval_real_f64(order.clone(), args, bessel_y_numeric_eval)
})
})
.register_tagged(|tags| {
let order = complex_float_tag("bessel_y", tags, 53).ok();
Box::new(move |args: &[Complex<f64>]| {
tagged_unary_eval_complex_f64(order.clone(), args, bessel_y_numeric_eval)
})
})
);
let bessel_i = symbol!(
"bessel_i",
norm = |x, out| {
if let Some([nu, z]) = function_arguments::<2>(x) {
if z.is_zero()
&& let Some(n) = atom_to_integer(nu)
&& n >= 0
{
if n == 0 {
out.to_num(Coefficient::one());
} else {
out.to_num(Coefficient::zero());
}
return;
}
let _ = maybe_eval_binary_float_in_norm(nu, z, out, bessel_i_numeric_eval);
}
},
der = move |x, i, out| {
if i == 1
&& let Some([nu, z]) = function_arguments::<2>(x)
{
let symbol = x.as_fun_view().unwrap().get_symbol();
let nu = nu.to_owned();
**out = (function!(symbol, nu.clone() - Atom::num(1), z)
+ function!(symbol, nu + Atom::num(1), z))
/ Atom::num(2);
}
},
eval = EvaluationInfo::new()
.with_tags(1)
.register_tagged(|tags| {
let tags = tags.iter().map(|x| x.to_owned()).collect::<Vec<_>>();
Box::new(move |args: &[Complex<Float>]| {
tagged_unary_eval_complex_float(
"bessel_i",
&tags,
args,
bessel_i_numeric_eval,
)
})
})
.register_tagged(|tags| {
let tags = tags.iter().map(|x| x.to_owned()).collect::<Vec<_>>();
Box::new(move |args: &[Float]| {
let prec = args.first().map(|x| x.prec()).unwrap_or(53);
let tag_views = tags.iter().map(|x| x.as_view()).collect::<Vec<_>>();
let Ok(order) = complex_float_tag("bessel_i", &tag_views, prec) else {
return Float::with_val(53, f64::NAN);
};
let [arg] = args else {
return Float::with_val(53, f64::NAN);
};
bessel_i_numeric_eval(
&order,
&Complex::new(arg.clone(), Float::new(prec)),
prec,
)
.map(|z| z.re)
.unwrap_or_else(|| Float::with_val(53, f64::NAN))
})
})
.register_tagged(|tags| {
let order = complex_float_tag("bessel_i", tags, 53).ok();
Box::new(move |args: &[f64]| {
tagged_unary_eval_real_f64(order.clone(), args, bessel_i_numeric_eval)
})
})
.register_tagged(|tags| {
let order = complex_float_tag("bessel_i", tags, 53).ok();
Box::new(move |args: &[Complex<f64>]| {
tagged_unary_eval_complex_f64(order.clone(), args, bessel_i_numeric_eval)
})
})
);
let bessel_k = symbol!(
"bessel_k",
norm = |x, out| {
if let Some([nu, z]) = function_arguments::<2>(x) {
if z.is_zero() {
out.to_num(Coefficient::complex_infinity());
return;
}
let _ = maybe_eval_binary_float_in_norm(nu, z, out, bessel_k_numeric_eval);
}
},
der = move |x, i, out| {
if i == 1
&& let Some([nu, z]) = function_arguments::<2>(x)
{
let symbol = x.as_fun_view().unwrap().get_symbol();
let nu = nu.to_owned();
**out = -(function!(symbol, nu.clone() - Atom::num(1), z)
+ function!(symbol, nu + Atom::num(1), z))
/ Atom::num(2);
}
},
eval = EvaluationInfo::new()
.with_tags(1)
.register_tagged(|tags| {
let tags = tags.iter().map(|x| x.to_owned()).collect::<Vec<_>>();
Box::new(move |args: &[Complex<Float>]| {
tagged_unary_eval_complex_float(
"bessel_k",
&tags,
args,
bessel_k_numeric_eval,
)
})
})
.register_tagged(|tags| {
let tags = tags.iter().map(|x| x.to_owned()).collect::<Vec<_>>();
Box::new(move |args: &[Float]| {
let prec = args.first().map(|x| x.prec()).unwrap_or(53);
let tag_views = tags.iter().map(|x| x.as_view()).collect::<Vec<_>>();
let Ok(order) = complex_float_tag("bessel_k", &tag_views, prec) else {
return Float::with_val(53, f64::NAN);
};
let [arg] = args else {
return Float::with_val(53, f64::NAN);
};
bessel_k_numeric_eval(
&order,
&Complex::new(arg.clone(), Float::new(prec)),
prec,
)
.map(|z| z.re)
.unwrap_or_else(|| Float::with_val(53, f64::NAN))
})
})
.register_tagged(|tags| {
let order = complex_float_tag("bessel_k", tags, 53).ok();
Box::new(move |args: &[f64]| {
tagged_unary_eval_real_f64(order.clone(), args, bessel_k_numeric_eval)
})
})
.register_tagged(|tags| {
let order = complex_float_tag("bessel_k", tags, 53).ok();
Box::new(move |args: &[Complex<f64>]| {
tagged_unary_eval_complex_f64(order.clone(), args, bessel_k_numeric_eval)
})
})
);
Self {
bessel_j,
bessel_y,
bessel_i,
bessel_k,
}
}
}
#[cfg(not(doctest))]
crate::_inventory::submit! {
StateInitializer::new("symbolica::special_functions", || {
let _ = GeometricSymbols::new();
let _ = SpecialSymbols::new();
let _ = BesselSymbols::new();
}, &["symbolica"])
}
pub fn gamma() -> Symbol {
SPECIALS.gamma
}
pub fn erf() -> Symbol {
SPECIALS.erf
}
pub fn euler_gamma() -> Symbol {
SPECIALS.euler_gamma
}
pub fn polygamma() -> Symbol {
SPECIALS.polygamma
}
pub fn polylog() -> Symbol {
SPECIALS.polylog
}
pub fn zeta() -> Symbol {
SPECIALS.zeta
}
pub fn tan() -> Symbol {
GEOMETRICS.tan
}
pub fn cot() -> Symbol {
GEOMETRICS.cot
}
pub fn sec() -> Symbol {
GEOMETRICS.sec
}
pub fn csc() -> Symbol {
GEOMETRICS.csc
}
pub fn asin() -> Symbol {
GEOMETRICS.asin
}
pub fn acos() -> Symbol {
GEOMETRICS.acos
}
pub fn atan() -> Symbol {
GEOMETRICS.atan
}
pub fn acot() -> Symbol {
GEOMETRICS.acot
}
pub fn asec() -> Symbol {
GEOMETRICS.asec
}
pub fn acsc() -> Symbol {
GEOMETRICS.acsc
}
pub fn sinh() -> Symbol {
GEOMETRICS.sinh
}
pub fn cosh() -> Symbol {
GEOMETRICS.cosh
}
pub fn tanh() -> Symbol {
GEOMETRICS.tanh
}
pub fn coth() -> Symbol {
GEOMETRICS.coth
}
pub fn sech() -> Symbol {
GEOMETRICS.sech
}
pub fn csch() -> Symbol {
GEOMETRICS.csch
}
pub fn asinh() -> Symbol {
GEOMETRICS.asinh
}
pub fn acosh() -> Symbol {
GEOMETRICS.acosh
}
pub fn atanh() -> Symbol {
GEOMETRICS.atanh
}
pub fn acoth() -> Symbol {
GEOMETRICS.acoth
}
pub fn asech() -> Symbol {
GEOMETRICS.asech
}
pub fn acsch() -> Symbol {
GEOMETRICS.acsch
}
pub fn bessel_j() -> Symbol {
BESSELS.bessel_j
}
pub fn bessel_y() -> Symbol {
BESSELS.bessel_y
}
pub fn bessel_i() -> Symbol {
BESSELS.bessel_i
}
pub fn bessel_k() -> Symbol {
BESSELS.bessel_k
}
macro_rules! unary_transcendental_methods {
($(
$(#[$meta:meta])*
$method:ident => $symbol:ident;
)*) => {
$(
$(#[$meta])*
fn $method(&self) -> Self::Output {
unary_transcendental_function(self, crate::transcendental::$symbol())
}
)*
};
}
pub trait TranscendentalFunctions: AtomCore {
unary_transcendental_methods! {
tan => tan;
cot => cot;
sec => sec;
csc => csc;
asin => asin;
acos => acos;
atan => atan;
acot => acot;
asec => asec;
acsc => acsc;
sinh => sinh;
cosh => cosh;
tanh => tanh;
coth => coth;
sech => sech;
csch => csch;
asinh => asinh;
acosh => acosh;
atanh => atanh;
acoth => acoth;
asech => asech;
acsch => acsch;
gamma => gamma;
erf => erf;
zeta => zeta;
}
fn atan2<'a, T: Into<AtomOrView<'a>>>(&self, y: T) -> Self::Output {
receiver_first_transcendental_function(self, crate::transcendental::atan(), y)
}
fn polygamma<'a, T: Into<AtomOrView<'a>>>(&self, n: T) -> Self::Output {
receiver_last_transcendental_function(self, crate::transcendental::polygamma(), n)
}
fn polylog<'a, T: Into<AtomOrView<'a>>>(&self, s: T) -> Self::Output {
receiver_last_transcendental_function(self, crate::transcendental::polylog(), s)
}
fn bessel_j<'a, T: Into<AtomOrView<'a>>>(&self, nu: T) -> Self::Output {
receiver_last_transcendental_function(self, crate::transcendental::bessel_j(), nu)
}
fn bessel_y<'a, T: Into<AtomOrView<'a>>>(&self, nu: T) -> Self::Output {
receiver_last_transcendental_function(self, crate::transcendental::bessel_y(), nu)
}
fn bessel_i<'a, T: Into<AtomOrView<'a>>>(&self, nu: T) -> Self::Output {
receiver_last_transcendental_function(self, crate::transcendental::bessel_i(), nu)
}
fn bessel_k<'a, T: Into<AtomOrView<'a>>>(&self, nu: T) -> Self::Output {
receiver_last_transcendental_function(self, crate::transcendental::bessel_k(), nu)
}
}
impl<T: AtomCore> TranscendentalFunctions for T {}
fn unary_transcendental_function<T: AtomCore>(arg: &T, symbol: Symbol) -> T::Output {
arg.atom_to_output(
FunctionBuilder::new(symbol)
.add_arg(arg.as_atom_view())
.finish(),
)
}
fn receiver_last_transcendental_function<'a, T, U>(
arg: &T,
symbol: Symbol,
first_arg: U,
) -> T::Output
where
T: AtomCore,
U: Into<AtomOrView<'a>>,
{
arg.atom_to_output(
FunctionBuilder::new(symbol)
.add_arg(first_arg)
.add_arg(arg.as_atom_view())
.finish(),
)
}
fn receiver_first_transcendental_function<'a, T, U>(
arg: &T,
symbol: Symbol,
second_arg: U,
) -> T::Output
where
T: AtomCore,
U: Into<AtomOrView<'a>>,
{
arg.atom_to_output(
FunctionBuilder::new(symbol)
.add_arg(arg.as_atom_view())
.add_arg(second_arg)
.finish(),
)
}
fn unary_eval_complex_float<F>(args: &[Complex<Float>], evaluator: F) -> Complex<Float>
where
F: FnOnce(&Complex<Float>, u32) -> Complex<Float>,
{
let prec = args
.first()
.map(|x| x.re.prec().max(x.im.prec()))
.unwrap_or(53);
let [arg] = args else {
return Complex::new(Float::with_val(53, f64::NAN), Float::with_val(53, f64::NAN));
};
evaluator(arg, prec)
}
fn tagged_unary_eval_to_float(
name: &str,
tags: &[AtomView],
args: &[Complex<Float>],
prec: u32,
evaluator: fn(&Complex<Float>, &Complex<Float>, u32) -> Option<Complex<Float>>,
) -> Result<Complex<Float>, String> {
let order = complex_float_tag(name, tags, prec)?;
let [arg] = args else {
return Err(format!(
"{name} expects exactly one argument after its order tag, got {}",
args.len()
));
};
evaluator(&order, arg, prec).ok_or_else(|| format!("{name} numeric evaluation failed"))
}
fn tagged_unary_eval_complex_float(
name: &str,
tags: &[Atom],
args: &[Complex<Float>],
evaluator: fn(&Complex<Float>, &Complex<Float>, u32) -> Option<Complex<Float>>,
) -> Complex<Float> {
let prec = args
.first()
.map(|x| x.re.prec().max(x.im.prec()))
.unwrap_or(53);
let tag_views = tags.iter().map(|x| x.as_view()).collect::<Vec<_>>();
tagged_unary_eval_to_float(name, &tag_views, args, prec, evaluator).unwrap_or_else(|_| {
Complex::new(Float::with_val(53, f64::NAN), Float::with_val(53, f64::NAN))
})
}
fn nonnegative_integer_tag(name: &str, tags: &[AtomView]) -> Result<u32, String> {
let [tag] = tags else {
return Err(format!(
"{name} expects exactly one order tag, got {}",
tags.len()
));
};
let Some(order) = atom_to_integer(*tag) else {
return Err(format!("{name} order tag must be an integer"));
};
if order < 0 || order > u32::MAX as i64 {
return Err(format!("{name} order tag must be a nonnegative integer"));
}
Ok(order as u32)
}
fn complex_float_tag(name: &str, tags: &[AtomView], prec: u32) -> Result<Complex<Float>, String> {
let [tag] = tags else {
return Err(format!(
"{name} expects exactly one order tag, got {}",
tags.len()
));
};
if let Ok(value) = Complex::<Float>::try_from(*tag) {
return Ok(value);
}
if let Ok(value) = Complex::<Rational>::try_from(*tag) {
return Ok(Complex::new(
value.re.to_multi_prec_float(prec),
value.im.to_multi_prec_float(prec),
));
}
Err(format!("{name} order tag must be numeric"))
}
fn unary_eval_real_f64(args: &[f64], evaluator: fn(f64) -> f64) -> f64 {
let [arg] = args else {
return f64::NAN;
};
evaluator(*arg)
}
fn unary_eval_real<T: Real, F: FnOnce(&T) -> T>(args: &[T], evaluator: F) -> T {
let [arg] = args else {
panic!(
"Wrong number of arguments for unary function: {}",
args.len()
);
};
evaluator(arg)
}
fn atan_eval_real<T: Real>(args: &[T]) -> T {
match args {
[z] => z.atan2(&z.one()),
[x, y] => y.atan2(x),
_ => panic!("Wrong number of arguments for atan: {}", args.len()),
}
}
fn unary_eval_f64(args: &[f64], evaluator: fn(&Complex<Float>, u32) -> Complex<Float>) -> f64 {
let [arg] = args else {
return f64::NAN;
};
evaluator(&Complex::new(Float::with_val(53, *arg), Float::new(53)), 53)
.re
.to_f64()
}
fn unary_eval_complex_f64(
args: &[Complex<f64>],
evaluator: fn(&Complex<Float>, u32) -> Complex<Float>,
) -> Complex<f64> {
let [arg] = args else {
return Complex::new(f64::NAN, f64::NAN);
};
evaluator(&complex_f64_to_float(arg, 53), 53).to_f64()
}
fn unary_eval_complex_real_f64<F>(args: &[Complex<f64>], evaluator: F) -> Complex<f64>
where
F: FnOnce(Complex<f64>) -> Complex<f64>,
{
let [arg] = args else {
return Complex::new(f64::NAN, f64::NAN);
};
evaluator(*arg)
}
fn tagged_unary_eval_real_f64(
order: Option<Complex<Float>>,
args: &[f64],
evaluator: fn(&Complex<Float>, &Complex<Float>, u32) -> Option<Complex<Float>>,
) -> f64 {
let Some(order) = order else {
return f64::NAN;
};
let [arg] = args else {
return f64::NAN;
};
evaluator(
&order,
&Complex::new(Float::with_val(53, *arg), Float::new(53)),
53,
)
.map(|z| z.re.to_f64())
.unwrap_or(f64::NAN)
}
fn tagged_unary_eval_complex_f64(
order: Option<Complex<Float>>,
args: &[Complex<f64>],
evaluator: fn(&Complex<Float>, &Complex<Float>, u32) -> Option<Complex<Float>>,
) -> Complex<f64> {
let Some(order) = order else {
return Complex::new(f64::NAN, f64::NAN);
};
let [arg] = args else {
return Complex::new(f64::NAN, f64::NAN);
};
evaluator(&order, &complex_f64_to_float(arg, 53), 53)
.map(|z| z.to_f64())
.unwrap_or_else(|| Complex::new(f64::NAN, f64::NAN))
}
fn atan_numeric_eval(z: &Complex<Float>, binary_prec: u32) -> Complex<Float> {
let i = z.i();
let one = complex_one(binary_prec);
let two = Complex::new(Float::with_val(binary_prec, 2), Float::new(binary_prec));
let iz = i.clone() * z;
(i / two) * ((one.clone() - &iz).log() - (one + iz).log())
}
fn atan2_numeric_eval(x: &Complex<Float>, y: &Complex<Float>, binary_prec: u32) -> Complex<Float> {
if y.im.is_zero() && x.im.is_zero() {
Complex::new(y.re.atan2(&x.re), Float::new(binary_prec))
} else {
y.atan2(x)
}
}
fn complex_one(prec: u32) -> Complex<Float> {
Complex::new(Float::with_val(prec, 1), Float::new(prec))
}
fn atan_eval_complex_f64(z: Complex<f64>) -> Complex<f64> {
let i = z.i();
let one = Complex::new(1.0, 0.0);
let two = Complex::new(2.0, 0.0);
let iz = i * z;
(i / two) * ((one - iz).log() - (one + iz).log())
}
fn bessel_j_numeric_eval(
order: &Complex<Float>,
z: &Complex<Float>,
binary_prec: u32,
) -> Option<Complex<Float>> {
if z.is_zero()
&& let Some(n) = complex_float_to_integer(order)
&& n >= 0
{
return Some(if n == 0 {
complex_one(binary_prec)
} else {
Complex::new(Float::new(binary_prec), Float::new(binary_prec))
});
}
if let Some(n) = complex_float_to_integer(order)
&& n < 0
{
let positive = Complex::new(Float::with_val(binary_prec, -n), Float::new(binary_prec));
let value = bessel_j_numeric_eval(&positive, z, binary_prec)?;
return Some(if n % 2 == 0 { value } else { -value });
}
bessel_series_eval(order, z, binary_prec, true)
}
fn bessel_i_numeric_eval(
order: &Complex<Float>,
z: &Complex<Float>,
binary_prec: u32,
) -> Option<Complex<Float>> {
if z.is_zero()
&& let Some(n) = complex_float_to_integer(order)
&& n >= 0
{
return Some(if n == 0 {
complex_one(binary_prec)
} else {
Complex::new(Float::new(binary_prec), Float::new(binary_prec))
});
}
if let Some(n) = complex_float_to_integer(order)
&& n < 0
{
let positive = Complex::new(Float::with_val(binary_prec, -n), Float::new(binary_prec));
return bessel_i_numeric_eval(&positive, z, binary_prec);
}
bessel_series_eval(order, z, binary_prec, false)
}
fn bessel_y_numeric_eval(
order: &Complex<Float>,
z: &Complex<Float>,
binary_prec: u32,
) -> Option<Complex<Float>> {
if z.is_zero() {
return Some(Complex::new(
Float::with_val(binary_prec, f64::INFINITY),
Float::new(binary_prec),
));
}
if let Some(n) = complex_float_to_integer(order)
&& n < 0
{
let positive = Complex::new(Float::with_val(binary_prec, -n), Float::new(binary_prec));
let value = bessel_y_numeric_eval(&positive, z, binary_prec)?;
return Some(if n % 2 == 0 { value } else { -value });
}
let order = bessel_regularized_order(order, binary_prec);
let pi_order = complex_pi(binary_prec) * order.clone();
let j_pos = bessel_j_numeric_eval(&order, z, binary_prec)?;
let j_neg = bessel_j_numeric_eval(&(-order.clone()), z, binary_prec)?;
Some((j_pos * pi_order.clone().cos() - j_neg) / pi_order.sin())
}
fn bessel_k_numeric_eval(
order: &Complex<Float>,
z: &Complex<Float>,
binary_prec: u32,
) -> Option<Complex<Float>> {
if z.is_zero() {
return Some(Complex::new(
Float::with_val(binary_prec, f64::INFINITY),
Float::new(binary_prec),
));
}
let order = if let Some(n) = complex_float_to_integer(order) {
Complex::new(
Float::with_val(binary_prec, n.abs()),
Float::new(binary_prec),
)
} else {
order.clone()
};
let order = bessel_regularized_order(&order, binary_prec);
let pi_order = complex_pi(binary_prec) * order.clone();
let i_neg = bessel_i_numeric_eval(&(-order.clone()), z, binary_prec)?;
let i_pos = bessel_i_numeric_eval(&order, z, binary_prec)?;
let pref = Complex::new(
Float::with_val(binary_prec, Constant::Pi) / Float::with_val(binary_prec, 2),
Float::new(binary_prec),
);
Some(pref * (i_neg - i_pos) / pi_order.sin())
}
fn bessel_series_eval(
order: &Complex<Float>,
z: &Complex<Float>,
binary_prec: u32,
alternating: bool,
) -> Option<Complex<Float>> {
let zero = Float::new(binary_prec);
let z_half = z.clone() / Complex::new(Float::with_val(binary_prec, 2), zero.clone());
let one = complex_one(binary_prec);
let gamma = gamma_numeric_eval(&(order.clone() + one), binary_prec);
let mut term = z_half.clone().powf(order) / gamma;
let mut sum = term.clone();
let z_half_sq = z_half.clone() * z_half;
let threshold = 2f64.powi(-(binary_prec.min(900) as i32));
for k in 1..(16 * binary_prec.max(16)) {
let kf = Float::with_val(binary_prec, k);
let denom = Complex::new(kf.clone(), zero.clone())
* (order.clone() + Complex::new(kf, zero.clone()));
let factor = if alternating {
-z_half_sq.clone()
} else {
z_half_sq.clone()
};
term = term * factor / denom;
let term_size = term.norm().re.to_f64().abs();
sum += term.clone();
if k > 16 && (term_size == 0.0 || term_size < threshold) {
return Some(sum);
}
}
Some(sum)
}
fn bessel_regularized_order(order: &Complex<Float>, binary_prec: u32) -> Complex<Float> {
if let Some(n) = complex_float_to_integer(order) {
let eps = bessel_order_epsilon(binary_prec);
Complex::new(
Float::with_val(binary_prec, n) + eps,
Float::new(binary_prec),
)
} else {
order.clone()
}
}
fn bessel_order_epsilon(binary_prec: u32) -> Float {
let eps = 2f64.powi(-((binary_prec.min(200) as i32) / 3));
Float::with_val(binary_prec, eps.max(1e-8))
}
fn complex_pi(prec: u32) -> Complex<Float> {
Complex::new(Float::with_val(prec, Constant::Pi), Float::new(prec))
}
fn polygamma_eval_f64(order: Option<u32>, args: &[f64]) -> f64 {
let Some(order) = order else {
return f64::NAN;
};
let [arg] = args else {
return f64::NAN;
};
polygamma_numeric_eval(
order,
&Complex::new(Float::with_val(53, *arg), Float::new(53)),
53,
)
.re
.to_f64()
}
fn polygamma_eval_complex_f64(order: Option<u32>, args: &[Complex<f64>]) -> Complex<f64> {
let Some(order) = order else {
return Complex::new(f64::NAN, f64::NAN);
};
let [arg] = args else {
return Complex::new(f64::NAN, f64::NAN);
};
polygamma_numeric_eval(order, &complex_f64_to_float(arg, 53), 53).to_f64()
}
fn polylog_eval_f64(order: Option<Complex<Float>>, args: &[f64]) -> f64 {
let Some(order) = order else {
return f64::NAN;
};
let [arg] = args else {
return f64::NAN;
};
let arg = Complex::new(Float::with_val(53, *arg), Float::new(53));
polylog_numeric_eval(&order, &arg, 53)
.map(|x| x.re.to_f64())
.unwrap_or(f64::NAN)
}
fn polylog_eval_complex_f64(order: Option<Complex<Float>>, args: &[Complex<f64>]) -> Complex<f64> {
let Some(order) = order else {
return Complex::new(f64::NAN, f64::NAN);
};
let [arg] = args else {
return Complex::new(f64::NAN, f64::NAN);
};
polylog_numeric_eval(&order, &complex_f64_to_float(arg, 53), 53)
.map(|x| x.to_f64())
.unwrap_or_else(|| Complex::new(f64::NAN, f64::NAN))
}
fn zeta_eval_f64(args: &[f64]) -> f64 {
let [arg] = args else {
return f64::NAN;
};
zeta_numeric_eval(&Complex::new(Float::with_val(53, *arg), Float::new(53)), 53)
.re
.to_f64()
}
fn zeta_eval_complex_f64(args: &[Complex<f64>]) -> Complex<f64> {
let [arg] = args else {
return Complex::new(f64::NAN, f64::NAN);
};
zeta_numeric_eval(&complex_f64_to_float(arg, 53), 53).to_f64()
}
fn complex_f64_to_float(value: &Complex<f64>, prec: u32) -> Complex<Float> {
Complex::new(
Float::with_val(prec, value.re),
Float::with_val(prec, value.im),
)
}
fn complex_float_to_integer(value: &Complex<Float>) -> Option<i64> {
if value.im.to_f64().abs() > 1e-12 {
return None;
}
let re = value.re.to_f64();
if !re.is_finite() {
return None;
}
let rounded = re.round();
if (re - rounded).abs() > 1e-12 || rounded < i64::MIN as f64 || rounded > i64::MAX as f64 {
return None;
}
Some(rounded as i64)
}
fn function_arguments<'a, const N: usize>(view: AtomView<'a>) -> Option<[AtomView<'a>; N]> {
let AtomView::Fun(fun) = view else {
return None;
};
fun.iter().try_into().ok()
}
fn gamma_exact_rational(rat: &Rational) -> Option<Atom> {
let denominator = rat.denominator();
if denominator == 1 {
let numerator = rat.numerator();
if numerator <= 0 {
return Some(Atom::num(Coefficient::complex_infinity()));
}
let n = numerator.to_i64()?;
if n > u32::MAX as i64 {
return None;
}
return Some(Atom::num(Integer::factorial((n - 1) as u32)));
}
if denominator == 2 {
let numerator = rat.numerator().to_i64()?;
if numerator % 2 == 0 {
return None;
}
let prefactor = gamma_half_integer_prefactor(numerator);
return Some(Atom::num(prefactor) * Atom::from(State::PI).pow(Atom::num((1, 2))));
}
None
}
fn zeta_exact_rational(rat: &Rational) -> Option<Atom> {
if rat.denominator() != 1 {
return None;
}
let n = rat.numerator().to_i64()?;
match n {
0 => return Some(Atom::num((-1, 2))),
1 => return Some(Atom::num(Coefficient::complex_infinity())),
-1 => return Some(Atom::num((-1, 12))),
_ => {}
}
if n < 0 {
let Some(bernoulli_index) = n.checked_neg().and_then(|x| x.checked_add(1)) else {
return None;
};
if bernoulli_index > u32::MAX as i64 {
return None;
}
let bernoulli_index = bernoulli_index as u32;
let value = bernoulli_number(bernoulli_index) / Rational::from((n - 1, 1));
if value.is_zero() {
return Some(Atom::new());
}
return Some(Atom::num(value));
}
if n % 2 != 0 {
return None;
}
if n > u32::MAX as i64 {
return None;
}
let half_order = n as u32 / 2;
let mut bernoulli = bernoulli_number(n as u32);
if half_order.is_multiple_of(2) {
bernoulli = -bernoulli;
}
let pi_power = (Atom::num(2) * Atom::from(State::PI)).pow(Atom::num(n));
Some(Atom::num(bernoulli) * pi_power / (Atom::num(2) * Atom::num(Integer::factorial(n as u32))))
}
fn bernoulli_number(n: u32) -> Rational {
let mut a = vec![Rational::zero(); n as usize + 1];
for m in 0..=n as usize {
a[m] = Rational::from((1, m as i64 + 1));
for j in (1..=m).rev() {
a[j - 1] = (a[j - 1].clone() - a[j].clone()) * Rational::from((j as i64, 1));
}
}
a[0].clone()
}
fn polygamma_order_zero_exact_rational(rat: &Rational) -> Option<Atom> {
let denominator = rat.denominator();
if denominator == 1 {
let numerator = rat.numerator();
if numerator <= 0 {
return Some(Atom::num(Coefficient::complex_infinity()));
}
let n = numerator.to_i64()?;
if n > u32::MAX as i64 {
return None;
}
return Some(Atom::num(harmonic_rational((n - 1) as u32)) - Atom::from(euler_gamma()));
}
if denominator == 2 {
let numerator = rat.numerator().to_i64()?;
if numerator % 2 == 0 {
return None;
}
return Some(
polygamma_order_zero_half_integer_base()
+ Atom::num(polygamma_order_zero_half_integer_shift(numerator)),
);
}
None
}
fn polygamma_exact_rational(order: u32, rat: &Rational) -> Option<Atom> {
if order == 0 {
return polygamma_order_zero_exact_rational(rat);
}
if *rat != Rational::one() {
return None;
}
match order {
1 => Some(Atom::from(State::PI).pow(2) / 6),
3 => Some(Atom::from(State::PI).pow(4) / 15),
_ => None,
}
}
fn gamma_half_integer_prefactor(numerator: i64) -> Rational {
let one = Rational::from((1, 1));
let target = Rational::from((numerator, 2));
let mut current = Rational::from((1, 2));
let mut prefactor = Rational::from((1, 1));
while current < target {
prefactor *= ¤t;
current += &one;
}
while current > target {
current -= &one;
prefactor /= ¤t;
}
prefactor
}
fn harmonic_rational(n: u32) -> Rational {
let mut sum = Rational::from((0, 1));
let one = Rational::from((1, 1));
for k in 1..=n {
sum += &(one.clone() / Rational::from((k as i64, 1)));
}
sum
}
fn polygamma_order_zero_half_integer_base() -> Atom {
-Atom::from(euler_gamma()) - Atom::num(2) * function!(State::LOG, Atom::num(2))
}
fn polygamma_order_zero_half_integer_shift(numerator: i64) -> Rational {
let one = Rational::from((1, 1));
let target = Rational::from((numerator, 2));
let mut current = Rational::from((1, 2));
let mut correction = Rational::from((0, 1));
while current < target {
correction += &(one.clone() / current.clone());
current += &one;
}
while current > target {
current -= &one;
correction -= &(one.clone() / current.clone());
}
correction
}
fn is_tan_pole(point: AtomView) -> bool {
let Some(ratio) = rational_multiple_of_pi(point) else {
return false;
};
ratio.denominator() == 2
}
fn is_cot_csc_pole(point: AtomView) -> bool {
rational_multiple_of_pi(point).is_some_and(|r| r.denominator() == 1)
}
fn is_tanh_sech_pole(point: AtomView) -> bool {
let Some(ratio) = rational_multiple_of_i_pi(point) else {
return false;
};
ratio.denominator() == 2
}
fn is_coth_csch_pole(point: AtomView) -> bool {
rational_multiple_of_i_pi(point).is_some_and(|r| r.denominator() == 1)
}
fn sec_residue(point: AtomView) -> Option<Atom> {
let ratio = rational_multiple_of_pi(point)?;
let sin = half_integer_sin_sign(&ratio)?;
Some(Atom::num(-sin))
}
fn csc_residue(point: AtomView) -> Option<Atom> {
let ratio = rational_multiple_of_pi(point)?;
let cos = integer_cos_sign(&ratio)?;
Some(Atom::num(cos))
}
fn sech_residue(point: AtomView) -> Option<Atom> {
let ratio = rational_multiple_of_i_pi(point)?;
let sin = half_integer_sin_sign(&ratio)?;
Some(Atom::num(-sin) * Atom::i())
}
fn csch_residue(point: AtomView) -> Option<Atom> {
let ratio = rational_multiple_of_i_pi(point)?;
let cos = integer_cos_sign(&ratio)?;
Some(Atom::num(cos))
}
fn half_integer_sin_sign(ratio: &Rational) -> Option<i64> {
if ratio.denominator() != 2 {
return None;
}
let rem = ratio.numerator() % 4;
if rem == 1 {
Some(1)
} else if rem == 3 {
Some(-1)
} else {
None
}
}
fn integer_cos_sign(ratio: &Rational) -> Option<i64> {
if ratio.denominator() != 1 {
return None;
}
if ratio.numerator() % 2 == 0 {
Some(1)
} else {
Some(-1)
}
}
fn rational_multiple_of_i_pi(point: AtomView) -> Option<Rational> {
let ratio = complex_rational_multiple_of_pi(point)?;
if ratio.re.is_zero() {
Some(ratio.im)
} else {
None
}
}
fn rational_multiple_of_pi(point: AtomView) -> Option<Rational> {
let ratio = complex_rational_multiple_of_pi(point)?;
if ratio.im.is_zero() {
Some(ratio.re)
} else {
None
}
}
fn complex_rational_multiple_of_pi(point: AtomView) -> Option<Complex<Rational>> {
match point {
AtomView::Num(_) => {
let r = Complex::<Rational>::try_from(point).ok()?;
if r.re.is_zero() && r.im.is_zero() {
Some(Complex::new(Rational::zero(), Rational::zero()))
} else {
None
}
}
AtomView::Var(v) => {
if v.get_symbol() == State::PI {
Some(Complex::new(Rational::one(), Rational::zero()))
} else {
None
}
}
AtomView::Mul(m) => {
let mut coeff = Complex::new(Rational::one(), Rational::zero());
let mut has_pi = false;
for factor in m {
if let Ok(r) = Complex::<Rational>::try_from(factor) {
coeff *= &r;
continue;
}
if let AtomView::Var(v) = factor
&& v.get_symbol() == State::PI
&& !has_pi
{
has_pi = true;
continue;
}
return None;
}
if has_pi { Some(coeff) } else { None }
}
AtomView::Add(a) => {
let mut sum = Complex::new(Rational::zero(), Rational::zero());
for term in a {
sum += &complex_rational_multiple_of_pi(term)?;
}
Some(sum)
}
_ => None,
}
}
fn is_negative_atom(arg: AtomView) -> bool {
match arg {
AtomView::Num(n) => n.get_coeff_view().to_owned().is_negative(),
AtomView::Mul(m) => m.into_iter().any(is_negative_atom),
_ => false,
}
}
fn atom_float_precision(value: AtomView) -> Option<u32> {
let AtomView::Num(number) = value else {
return None;
};
match number.get_coeff_view() {
CoefficientView::Float(r, i) => Some(r.to_float().prec().max(i.to_float().prec())),
_ => None,
}
}
fn atom_to_integer(value: AtomView) -> Option<i64> {
let rat = Rational::try_from(value).ok()?;
if !rat.is_integer() {
return None;
}
rat.numerator().to_i64()
}
fn maybe_eval_unary_float_in_norm(
arg: AtomView,
out: &mut Settable<Atom>,
evaluator: fn(&Complex<Float>, u32) -> Complex<Float>,
) -> bool {
let Some(prec) = atom_float_precision(arg) else {
return false;
};
let Some(z) = atom_to_complex_float(arg, prec) else {
return false;
};
out.to_num(evaluator(&z, prec));
true
}
fn maybe_eval_polygamma_float_in_norm(order: u32, arg: AtomView, out: &mut Settable<Atom>) -> bool {
let Some(prec) = atom_float_precision(arg) else {
return false;
};
let Some(z) = atom_to_complex_float(arg, prec) else {
return false;
};
out.to_num(polygamma_numeric_eval(order, &z, prec));
true
}
fn maybe_eval_binary_float_in_norm(
lhs: AtomView,
rhs: AtomView,
out: &mut Settable<Atom>,
evaluator: fn(&Complex<Float>, &Complex<Float>, u32) -> Option<Complex<Float>>,
) -> bool {
let precision = atom_float_precision(lhs).or_else(|| atom_float_precision(rhs));
let Some(prec) = precision else {
return false;
};
let Some(lhs_float) = atom_to_complex_float(lhs, prec) else {
return false;
};
let Some(rhs_float) = atom_to_complex_float(rhs, prec) else {
return false;
};
let Some(value) = evaluator(&lhs_float, &rhs_float, prec) else {
return false;
};
out.to_num(value);
true
}
fn atom_to_complex_float(value: AtomView, binary_prec: u32) -> Option<Complex<Float>> {
let AtomView::Num(number) = value else {
return None;
};
match number.get_coeff_view() {
CoefficientView::Natural(nr, dr, ni, di) => Some(Complex::new(
Float::with_val(binary_prec, nr) / Float::with_val(binary_prec, dr),
Float::with_val(binary_prec, ni) / Float::with_val(binary_prec, di),
)),
CoefficientView::Large(r, i) => Some(Complex::new(
r.to_rat().to_multi_prec_float(binary_prec),
i.to_rat().to_multi_prec_float(binary_prec),
)),
CoefficientView::Float(r, i) => {
let mut re = r.to_float();
let mut im = i.to_float();
if re.prec() > binary_prec {
re.set_prec(binary_prec);
}
if im.prec() > binary_prec {
im.set_prec(binary_prec);
}
Some(Complex::new(re, im))
}
_ => None,
}
}
fn gamma_numeric_eval(z: &Complex<Float>, binary_prec: u32) -> Complex<Float> {
#[cfg(feature = "gmp")]
{
if z.im.to_f64() == 0.0 {
return Complex::new(
z.re.clone().into_inner().gamma().into(),
Float::new(binary_prec),
);
}
}
gamma_complex_spouge(z, binary_prec)
}
fn erf_numeric_eval(z: &Complex<Float>, binary_prec: u32) -> Complex<Float> {
if z.is_zero() {
return Complex::new(Float::new(binary_prec), Float::new(binary_prec));
}
if z.re.to_f64() < 0.0 {
return -erf_numeric_eval(&(-z.clone()), binary_prec);
}
let abs_z = z.norm().re.to_f64().abs();
if abs_z <= 4.0 || z.re.to_f64().abs() < 1.0 {
erf_series_eval(z, binary_prec)
} else {
erf_asymptotic_eval(z, binary_prec)
}
}
fn polygamma_order_zero_numeric_eval(z: &Complex<Float>, binary_prec: u32) -> Complex<Float> {
#[cfg(feature = "gmp")]
{
if z.im.to_f64() == 0.0 {
return Complex::new(
z.re.clone().into_inner().digamma().into(),
Float::new(binary_prec),
);
}
}
polygamma_order_zero_complex(z, binary_prec)
}
fn polygamma_numeric_eval(order: u32, z: &Complex<Float>, binary_prec: u32) -> Complex<Float> {
if order == 0 {
return polygamma_order_zero_numeric_eval(z, binary_prec);
}
let zero = Float::new(binary_prec);
let one = Float::with_val(binary_prec, 1);
let exponent = Complex::new(Float::with_val(binary_prec, order + 1), zero.clone());
let factorial = Float::with_val(binary_prec, Integer::factorial(order).to_multi_prec());
let sign = if order.is_multiple_of(2) {
-factorial
} else {
factorial
};
let prefactor = Complex::new(sign, zero.clone());
let mut shifted = z.clone();
let mut correction = Complex::new(zero.clone(), zero.clone());
while shifted.re.to_f64() < 8.0 {
let denom = shifted.clone().powf(&exponent);
correction += prefactor.clone() / denom;
shifted += Complex::new(one.clone(), zero.clone());
}
let threshold = 2f64.powi(-(binary_prec.min(900) as i32));
let mut sum = Complex::new(zero.clone(), zero.clone());
for n in 0..(8 * binary_prec.max(16)) {
let shift = Complex::new(Float::with_val(binary_prec, n), zero.clone());
let term = prefactor.clone() / (shifted.clone() + shift).powf(&exponent);
let term_size = term.norm().re.to_f64().abs();
sum += term;
if n > 16 && (term_size == 0.0 || term_size < threshold) {
break;
}
}
correction + sum
}
fn polylog_exact(s: AtomView, z: AtomView) -> Option<Atom> {
if z.is_zero() {
return Some(Atom::new());
}
if let Some(order) = atom_to_integer(s) {
if order < 0 {
return polylog_negative_integer_exact(order, z);
}
if order == 0 {
return Some(z.to_owned() / (Atom::num(1) - z));
}
if order == 1 {
return Some(-function!(State::LOG, Atom::num(1) - z));
}
}
if z.is_one() {
if let Some(order) = atom_to_integer(s)
&& order == 1
{
return Some(Atom::num(Coefficient::complex_infinity()));
}
if let Ok(rat) = Rational::try_from(s)
&& let Some(exact) = zeta_exact_rational(&rat)
{
return Some(exact);
}
return Some(function!(SPECIALS.zeta, s.to_owned()));
}
if let Ok(rat) = Rational::try_from(z)
&& rat == (-1, 1)
{
if let Some(order) = atom_to_integer(s)
&& order == 1
{
return Some(-function!(State::LOG, Atom::num(2)));
}
if let Ok(s_rat) = Rational::try_from(s)
&& let Some(exact) = polylog_minus_one_exact_rational(&s_rat)
{
return Some(exact);
}
return Some(
-(Atom::num(1) - Atom::num(2).pow(Atom::num(1) - s.to_owned()))
* function!(SPECIALS.zeta, s.to_owned()),
);
}
None
}
fn polylog_numeric_eval(
s: &Complex<Float>,
z: &Complex<Float>,
binary_prec: u32,
) -> Option<Complex<Float>> {
let zero = Float::new(binary_prec);
let one = Float::with_val(binary_prec, 1);
if z.norm().re.to_f64() == 0.0 {
return Some(Complex::new(zero.clone(), zero.clone()));
}
if s.im.to_f64() == 0.0 && s.re.to_f64() == 0.0 {
return Some(z.clone() / (Complex::new(one.clone(), zero.clone()) - z.clone()));
}
if s.im.to_f64() == 0.0 && s.re.to_f64() == 1.0 {
return Some(-(Complex::new(one.clone(), zero.clone()) - z.clone()).log());
}
if s.im.to_f64() == 0.0 {
let rounded = s.re.to_f64().round();
if (s.re.to_f64() - rounded).abs() < 1e-14
&& rounded >= i64::MIN as f64
&& rounded <= i64::MAX as f64
{
return polylog_integer_numeric_eval(rounded as i64, z, binary_prec, 0);
}
}
if z.norm().re.to_f64() >= 0.95 {
return None;
}
let threshold = 2f64.powi(-(binary_prec.min(900) as i32));
let mut z_pow = z.clone();
let mut sum = Complex::new(zero.clone(), zero.clone());
for k in 1..(12 * binary_prec.max(16)) {
let base = Complex::new(Float::with_val(binary_prec, k), zero.clone());
let denom = base.powf(s);
let term = z_pow.clone() / denom;
let term_size = term.norm().re.to_f64().abs();
sum += term;
if k > 16 && (term_size == 0.0 || term_size < threshold) {
return Some(sum);
}
z_pow *= z.clone();
}
Some(sum)
}
fn polylog_minus_one_exact_rational(s: &Rational) -> Option<Atom> {
if *s == Rational::one() {
return Some(-function!(State::LOG, Atom::num(2)));
}
let zeta_term = if let Some(exact) = zeta_exact_rational(s) {
exact
} else {
function!(SPECIALS.zeta, Atom::num(s.clone()))
};
Some(-(Atom::num(1) - Atom::num(2).pow(Atom::num(1) - Atom::num(s.clone()))) * zeta_term)
}
fn polylog_negative_integer_exact(order: i64, z: AtomView) -> Option<Atom> {
let n = order.checked_neg()? as u32;
let coeffs = eulerian_row(n);
let mut numerator = Atom::new();
for (k, coeff) in coeffs.into_iter().enumerate() {
numerator += Atom::num(coeff) * z.to_owned().pow(Atom::num((k + 1) as i64));
}
Some(numerator / (Atom::num(1) - z).pow(Atom::num((n + 1) as i64)))
}
fn eulerian_row(n: u32) -> Vec<Integer> {
if n == 0 {
return vec![Integer::from(1)];
}
let mut row = vec![Integer::from(1)];
for m in 2..=n {
let mut next = vec![Integer::from(0); m as usize];
for k in 0..m as usize {
let mut value = Integer::from(0);
if k > 0 {
value += Integer::from(m - k as u32) * row[k - 1].clone();
}
if k < row.len() {
value += Integer::from(k as u32 + 1) * row[k].clone();
}
next[k] = value;
}
row = next;
}
row
}
fn polylog_integer_numeric_eval(
order: i64,
z: &Complex<Float>,
binary_prec: u32,
depth: u32,
) -> Option<Complex<Float>> {
if depth > 8 {
return None;
}
let zero = Float::new(binary_prec);
if order <= 0 {
return Some(polylog_negative_integer_numeric_eval(order, z, binary_prec));
}
if z.im.to_f64() == 0.0 && z.re.to_f64() == 1.0 {
return Some(zeta_integer_numeric_eval(order, binary_prec));
}
if z.im.to_f64() == 0.0 && z.re.to_f64() == -1.0 {
let zeta = zeta_integer_numeric_eval(order, binary_prec);
let exponent = Complex::new(Float::with_val(binary_prec, 1 - order), zero.clone());
let factor = complex_one(binary_prec)
- Complex::new(Float::with_val(binary_prec, 2), zero.clone()).powf(&exponent);
return Some(-factor * zeta);
}
if order == 1 {
return Some(-(complex_one(binary_prec) - z.clone()).log());
}
#[cfg(feature = "gmp")]
if order == 2 && z.im.to_f64() == 0.0 && z.re.to_f64() <= 1.0 {
return Some(Complex::new(
z.re.clone().into_inner().li2().into(),
zero.clone(),
));
}
let abs_z = z.norm().re.to_f64().abs();
if abs_z < 0.95 {
return polylog_series_integer(order as u32, z, binary_prec, 64);
}
let log_z = z.log();
if log_z.norm().re.to_f64().abs() < 1.0 {
return Some(polylog_real_branch_fix(
z,
polylog_jonquiere_integer(order as u32, &log_z, binary_prec),
));
}
if abs_z <= 1.0 + 1e-14 {
return polylog_series_integer(order as u32, z, binary_prec, 64);
}
let inv = z.inv();
if order == 2 {
let continued = polylog_integer_numeric_eval(order, &inv, binary_prec, depth + 1)?;
let log_minus_z = polylog_log_minus(z, binary_prec);
let pi_sq_over_six = Complex::new(
Float::with_val(binary_prec, Constant::Pi).pow(2) / Float::with_val(binary_prec, 6),
zero.clone(),
);
return Some(polylog_real_branch_fix(
z,
-continued
- pi_sq_over_six
- log_minus_z.clone() * log_minus_z
/ Complex::new(Float::with_val(binary_prec, 2), zero.clone()),
));
}
let continued = polylog_integer_numeric_eval(order, &inv, binary_prec, depth + 1)?;
let log_minus_z = polylog_log_minus(z, binary_prec);
let i = log_minus_z.i();
let bernoulli = bernoulli_polynomial(
order as u32,
&(Complex::new(Float::with_val(binary_prec, 1), zero.clone())
/ Complex::new(Float::with_val(binary_prec, 2), zero.clone())
+ log_minus_z.clone()
/ (Complex::new(Float::with_val(binary_prec, 2), zero.clone())
* Complex::new(Float::with_val(binary_prec, Constant::Pi), zero.clone())
* i.clone())),
binary_prec,
);
let two_pi_i = Complex::new(
Float::with_val(binary_prec, 2) * Float::with_val(binary_prec, Constant::Pi),
zero.clone(),
) * i;
let prefactor = -pow_complex_u32(&two_pi_i, order as u32)
/ Complex::new(factorial_float(order as u32, binary_prec), zero.clone());
if order % 2 == 0 {
Some(polylog_real_branch_fix(
z,
prefactor * bernoulli - continued,
))
} else {
Some(polylog_real_branch_fix(
z,
prefactor * bernoulli + continued,
))
}
}
fn polylog_negative_integer_numeric_eval(
order: i64,
z: &Complex<Float>,
binary_prec: u32,
) -> Complex<Float> {
let n = order.checked_neg().unwrap_or_default() as u32;
let coeffs = eulerian_row(n);
let mut poly = Complex::new(Float::new(binary_prec), Float::new(binary_prec));
for coeff in coeffs.iter().rev() {
poly = z.clone()
* (poly
+ Complex::new(
Float::with_val(binary_prec, coeff.clone().to_multi_prec()),
Float::new(binary_prec),
));
}
poly / pow_complex_u32(&(complex_one(binary_prec) - z.clone()), n + 1)
}
fn polylog_series_integer(
order: u32,
z: &Complex<Float>,
binary_prec: u32,
max_terms_factor: u32,
) -> Option<Complex<Float>> {
let zero = Float::new(binary_prec);
let threshold = 2f64.powi(-(binary_prec.min(900) as i32));
let mut z_pow = z.clone();
let mut sum = Complex::new(zero.clone(), zero.clone());
for k in 1..(max_terms_factor * binary_prec.max(16)) {
let denom = Float::with_val(binary_prec, k).pow(order as u64);
let term = z_pow.clone() / Complex::new(denom, zero.clone());
let term_size = term.norm().re.to_f64().abs();
sum += term;
if k > 16 && (term_size == 0.0 || term_size < threshold) {
return Some(sum);
}
z_pow *= z.clone();
}
Some(sum)
}
fn polylog_jonquiere_integer(order: u32, mu: &Complex<Float>, binary_prec: u32) -> Complex<Float> {
let zero = Float::new(binary_prec);
let threshold = 2f64.powi(-(binary_prec.min(900) as i32));
let mut sum = Complex::new(zero.clone(), zero.clone());
let mut mu_pow = Complex::new(Float::with_val(binary_prec, 1), zero.clone());
let mut factorial = Float::with_val(binary_prec, 1);
for k in 0u32..(32 * binary_prec.max(16)) {
let term = if k + 1 == order {
mu_pow.clone() / Complex::new(factorial.clone(), zero.clone())
* Complex::new(
harmonic_rational(order - 1).to_multi_prec_float(binary_prec),
zero.clone(),
)
- mu_pow.clone() / Complex::new(factorial.clone(), zero.clone())
* (-mu.clone()).log()
} else {
mu_pow.clone() / Complex::new(factorial.clone(), zero.clone())
* zeta_integer_numeric_eval(order as i64 - k as i64, binary_prec)
};
let term_size = term.norm().re.to_f64().abs();
sum += term;
if k > order + 8 && (term_size == 0.0 || term_size < threshold) {
return sum;
}
mu_pow *= mu.clone();
factorial *= Float::with_val(binary_prec, k + 1);
}
sum
}
fn zeta_integer_numeric_eval(order: i64, binary_prec: u32) -> Complex<Float> {
zeta_numeric_eval(
&Complex::new(Float::with_val(binary_prec, order), Float::new(binary_prec)),
binary_prec,
)
}
fn bernoulli_polynomial(n: u32, x: &Complex<Float>, binary_prec: u32) -> Complex<Float> {
let zero = Float::new(binary_prec);
let mut sum = Complex::new(zero.clone(), zero.clone());
for k in 0..=n {
let coeff = Float::with_val(
binary_prec,
Integer::binom(n.into(), k.into()).to_multi_prec(),
) * bernoulli_number(n - k).to_multi_prec_float(binary_prec);
sum += Complex::new(coeff, zero.clone()) * pow_complex_u32(x, k);
}
sum
}
fn polylog_log_minus(z: &Complex<Float>, binary_prec: u32) -> Complex<Float> {
if z.im.to_f64() == 0.0 && z.re.to_f64() > 0.0 {
return z.log()
+ Complex::new(
Float::new(binary_prec),
Float::with_val(binary_prec, Constant::Pi),
);
}
(-z.clone()).log()
}
fn polylog_real_branch_fix(z: &Complex<Float>, value: Complex<Float>) -> Complex<Float> {
if z.im.to_f64() == 0.0 && z.re.to_f64() > 1.0 && value.im.to_f64() > 0.0 {
return Complex::new(value.re, -value.im);
}
value
}
fn pow_complex_u32(z: &Complex<Float>, n: u32) -> Complex<Float> {
if n == 0 {
return complex_one(z.re.prec().max(z.im.prec()));
}
let prec = z.re.prec().max(z.im.prec());
let mut base = z.clone();
let mut exp = n;
let mut out = complex_one(prec);
while exp > 0 {
if exp % 2 == 1 {
out *= base.clone();
}
exp /= 2;
if exp > 0 {
base *= base.clone();
}
}
out
}
fn factorial_float(n: u32, binary_prec: u32) -> Float {
Float::with_val(binary_prec, Integer::factorial(n).to_multi_prec())
}
fn zeta_numeric_eval(s: &Complex<Float>, binary_prec: u32) -> Complex<Float> {
let zero = Float::new(binary_prec);
if s.is_real() {
if s.re.to_f64() == 1.0 {
return Complex::new(Float::with_val(binary_prec, f64::INFINITY), zero);
}
#[cfg(feature = "gmp")]
{
return Complex::new(s.re.clone().into_inner().zeta().into(), zero);
}
#[cfg(not(feature = "gmp"))]
{
return zeta_complex_hasse(s, binary_prec);
}
}
if (s.re.to_f64() - 1.0).abs() < 1e-14 && s.im.to_f64().abs() < 1e-14 {
return Complex::new(Float::with_val(binary_prec, f64::INFINITY), zero);
}
let work_prec = binary_prec
.saturating_mul(4)
.max(binary_prec.saturating_add(64));
let mut s_re = s.re.clone();
let mut s_im = s.im.clone();
s_re.set_prec(work_prec);
s_im.set_prec(work_prec);
let s = Complex::new(s_re, s_im);
let result = if s.re.to_f64() < 0.0 {
zeta_complex_reflection(&s, work_prec)
} else {
zeta_complex_hasse(&s, work_prec)
};
let mut re = result.re;
let mut im = result.im;
re.set_prec(binary_prec);
im.set_prec(binary_prec);
Complex::new(re, im)
}
fn zeta_complex_reflection(s: &Complex<Float>, binary_prec: u32) -> Complex<Float> {
let zero = Float::new(binary_prec);
let one = Float::with_val(binary_prec, 1);
let two = Float::with_val(binary_prec, 2);
let pi = Float::with_val(binary_prec, Constant::Pi);
let one_c = Complex::new(one.clone(), zero.clone());
let two_c = Complex::new(two, zero.clone());
let pi_c = Complex::new(pi.clone(), zero.clone());
let reflected = one_c.clone() - s.clone();
let two_power = two_c.powf(s);
let pi_power = pi_c.clone().powf(&(s.clone() - one_c));
let sine =
(pi_c * s.clone() / Complex::new(Float::with_val(binary_prec, 2), zero.clone())).sin();
let gamma = gamma_numeric_eval(&reflected, binary_prec);
#[cfg(feature = "gmp")]
let zeta = if reflected.im.to_f64() == 0.0 {
Complex::new(
reflected.re.clone().into_inner().zeta().into(),
zero.clone(),
)
} else {
zeta_complex_hasse(&reflected, binary_prec)
};
#[cfg(not(feature = "gmp"))]
let zeta = zeta_complex_hasse(&reflected, binary_prec);
two_power * pi_power * sine * gamma * zeta
}
fn zeta_complex_hasse(s: &Complex<Float>, binary_prec: u32) -> Complex<Float> {
let zero = Float::new(binary_prec);
let one = Float::with_val(binary_prec, 1);
let two = Float::with_val(binary_prec, 2);
let one_c = Complex::new(one.clone(), zero.clone());
let two_c = Complex::new(two.clone(), zero.clone());
let denominator = one_c.clone() - two_c.powf(&(one_c - s.clone()));
let threshold = 2f64.powi(-(binary_prec.min(900) as i32));
let mut outer = Complex::new(zero.clone(), zero.clone());
let mut two_power = Float::with_val(binary_prec, 2);
for n in 0..(8 * binary_prec.max(16)) {
let mut inner = Complex::new(zero.clone(), zero.clone());
let mut binom = Float::with_val(binary_prec, 1);
for k in 0..=n {
let base = Complex::new(Float::with_val(binary_prec, k + 1), zero.clone());
let mut term = Complex::new(binom.clone(), zero.clone()) / base.powf(s);
if k % 2 == 1 {
term = -term;
}
inner += term;
if k < n {
binom *= Float::with_val(binary_prec, n - k);
binom /= Float::with_val(binary_prec, k + 1);
}
}
let term = inner / Complex::new(two_power.clone(), zero.clone());
let term_size = term.norm().re.to_f64().abs();
outer += term;
if n > 16 && (term_size == 0.0 || term_size < threshold) {
break;
}
two_power *= Float::with_val(binary_prec, 2);
}
outer / denominator
}
fn gamma_complex_spouge(z: &Complex<Float>, binary_prec: u32) -> Complex<Float> {
let zero = Float::new(binary_prec);
let one = Float::with_val(binary_prec, 1);
let half = Float::with_val(binary_prec, 1) / Float::with_val(binary_prec, 2);
let pi = Float::with_val(binary_prec, Constant::Pi);
if z.re.to_f64() < 0.5 {
let numerator = Complex::new(pi.clone(), zero.clone());
let pi_z = Complex::new(pi, zero.clone()) * z.clone();
let reflected = Complex::new(one.clone(), zero.clone()) - z.clone();
return numerator / (pi_z.sin() * gamma_complex_spouge(&reflected, binary_prec));
}
let a = spouge_parameter(binary_prec);
let mut sum = Complex::new(
(Float::with_val(binary_prec, 2) * Float::with_val(binary_prec, Constant::Pi)).sqrt(),
zero.clone(),
);
for k in 1..a {
let coeff = spouge_coefficient(a, k, binary_prec);
let denom = z.clone() + Complex::new(Float::with_val(binary_prec, k - 1), zero.clone());
sum += Complex::new(coeff, zero.clone()) / denom;
}
let t = z.clone() + Complex::new(Float::with_val(binary_prec, a - 1), zero.clone());
let exponent = z.clone() - Complex::new(half, zero);
t.powf(&exponent) * (-t).exp() * sum
}
fn polygamma_order_zero_complex(z: &Complex<Float>, binary_prec: u32) -> Complex<Float> {
let zero = Float::new(binary_prec);
let one = Float::with_val(binary_prec, 1);
let pi = Float::with_val(binary_prec, Constant::Pi);
if z.re.to_f64() < 0.5 {
let pi_z = Complex::new(pi.clone(), zero.clone()) * z.clone();
let reflected = Complex::new(one.clone(), zero.clone()) - z.clone();
return polygamma_order_zero_complex(&reflected, binary_prec)
- Complex::new(pi, zero.clone()) / pi_z.tan();
}
let mut shifted = z.clone();
let mut correction = Complex::new(zero.clone(), zero.clone());
while shifted.re.to_f64() < 8.0 {
correction -= Complex::new(one.clone(), zero.clone()) / shifted.clone();
shifted += Complex::new(one.clone(), zero.clone());
}
let mut sum = Complex::new(-Float::with_val(binary_prec, Constant::Euler), zero.clone());
let threshold = 2f64.powi(-(binary_prec.min(900) as i32));
for n in 0..(8 * binary_prec.max(16)) {
let n1 = Float::with_val(binary_prec, n + 1);
let term = Complex::new(one.clone() / &n1, zero.clone())
- Complex::new(one.clone(), zero.clone())
/ (shifted.clone() + Complex::new(n1, zero.clone()));
let term_size = term.norm().re.to_f64().abs();
sum += term;
if n > 16 && (term_size == 0.0 || term_size < threshold) {
break;
}
}
correction + sum
}
fn erf_series_eval(z: &Complex<Float>, binary_prec: u32) -> Complex<Float> {
let work_prec = binary_prec
.saturating_mul(2)
.max(binary_prec.saturating_add(64));
let mut z = z.clone();
z.re.set_prec(work_prec);
z.im.set_prec(work_prec);
let zero = Float::new(work_prec);
let z_squared = z.clone() * z.clone();
let mut term = z.clone();
let mut sum = term.clone();
let threshold = 2f64.powi(-(binary_prec.min(900) as i32));
for n in 0..(64 * binary_prec.max(16)) {
let numerator = Complex::new(Float::with_val(work_prec, 2 * n + 1), zero.clone());
let denominator = Complex::new(
Float::with_val(work_prec, n + 1) * Float::with_val(work_prec, 2 * n + 3),
zero.clone(),
);
term = -term * z_squared.clone() * numerator / denominator;
let term_size = term.norm().re.to_f64().abs();
sum += term.clone();
if n > 16 && (term_size == 0.0 || term_size < threshold) {
break;
}
}
let prefactor = Complex::new(
Float::with_val(work_prec, 2) / Float::with_val(work_prec, Constant::Pi).sqrt(),
zero,
);
let mut result = prefactor * sum;
result.re.set_prec(binary_prec);
result.im.set_prec(binary_prec);
result
}
fn erf_asymptotic_eval(z: &Complex<Float>, binary_prec: u32) -> Complex<Float> {
let work_prec = binary_prec
.saturating_mul(2)
.max(binary_prec.saturating_add(64));
let mut z = z.clone();
z.re.set_prec(work_prec);
z.im.set_prec(work_prec);
let zero = Float::new(work_prec);
let one = Complex::new(Float::with_val(work_prec, 1), zero.clone());
let two = Complex::new(Float::with_val(work_prec, 2), zero.clone());
let z_squared = z.clone() * z.clone();
let mut term = one.clone();
let mut sum = term.clone();
let threshold = 2f64.powi(-(binary_prec.min(900) as i32));
for n in 0..(16 * binary_prec.max(16)) {
let factor = Complex::new(Float::with_val(work_prec, 2 * n + 1), zero.clone())
/ (two.clone() * z_squared.clone());
term = -term * factor;
let term_size = term.norm().re.to_f64().abs();
if n > 0 && term_size > sum.norm().re.to_f64().abs() {
break;
}
sum += term.clone();
if n > 4 && (term_size == 0.0 || term_size < threshold) {
break;
}
}
let sqrt_pi = Float::with_val(work_prec, Constant::Pi).sqrt();
let erfc = (-z_squared).exp() * sum / (Complex::new(sqrt_pi, zero) * z);
let mut result = one - erfc;
result.re.set_prec(binary_prec);
result.im.set_prec(binary_prec);
result
}
fn spouge_parameter(binary_prec: u32) -> u32 {
((binary_prec as f64 / (2.0 * std::f64::consts::PI).log2()).ceil() as u32).max(12) + 2
}
fn spouge_coefficient(a: u32, k: u32, binary_prec: u32) -> Float {
let half = Float::with_val(binary_prec, 1) / Float::with_val(binary_prec, 2);
let a_minus_k = Float::with_val(binary_prec, a - k);
let exponent = Float::with_val(binary_prec, k) - half;
let factorial = Float::with_val(binary_prec, Integer::factorial(k - 1).to_multi_prec());
let mut coeff = a_minus_k.powf(&exponent) * a_minus_k.exp() / factorial;
if k.is_multiple_of(2) {
coeff = -coeff;
}
coeff
}
#[cfg(test)]
mod tests {
use std::{
io::Write,
process::{Command, Stdio},
};
use crate::{
atom::{Atom, AtomCore},
coefficient::Coefficient,
domains::float::{Complex, Float, Real, RealLike},
parse, symbol,
};
#[test]
fn gamma_exact_normalization() {
assert_eq!(parse!("gamma(5)"), Atom::num(24));
assert_eq!(parse!("gamma(1/2)"), parse!("pi^(1/2)"));
assert_eq!(parse!("gamma(-5/2)"), parse!("-8/15*pi^(1/2)"));
assert_eq!(parse!("gamma(0.3)"), parse!("2.991568987687591"));
assert_eq!(
parse!("gamma(0)"),
Atom::num(Coefficient::complex_infinity())
);
}
#[test]
fn gamma_derivative() {
assert_eq!(
parse!("gamma(x)").derivative(symbol!("x")),
parse!("gamma(x)*polygamma(0,x)")
);
}
#[test]
fn erf_basic_normalization() {
assert_eq!(parse!("erf(0)"), Atom::new());
assert_eq!(parse!("erf(-x)"), parse!("-erf(x)"));
let actual = Complex::<Float>::try_from(parse!("erf(1)").to_float(80)).unwrap();
let expected = Complex::new(
Float::parse("0.84270079294971486934", Some(80)).unwrap(),
Float::new(80),
);
assert_close_complex(&actual, &expected, "1e-20");
}
#[test]
fn erf_derivative() {
assert_eq!(
parse!("erf(x)").derivative(symbol!("x")),
parse!("2*exp(-x^2)/sqrt(pi)")
);
}
#[test]
fn gamma_laurent_series() {
let x = symbol!("x");
assert_eq!(
parse!("gamma(x)").series(x, 0, 3).unwrap().to_atom(),
parse!(
"x^-1-euler_gamma+1/2*(euler_gamma^2+1/6*pi^2)*x+1/6*(-euler_gamma^3-1/2*euler_gamma*pi^2+polygamma(2,1))*x^2+1/24*(euler_gamma^4+euler_gamma^2*pi^2-4*euler_gamma*polygamma(2,1)+3/20*pi^4)*x^3"
)
);
assert_eq!(
parse!("gamma(x)").series(x, 0, 0).unwrap().to_atom(),
parse!("x^-1-euler_gamma")
);
assert_eq!(
parse!("gamma(x-1)").series(x, 1, 0).unwrap().to_atom(),
parse!("(x-1)^-1-euler_gamma")
);
}
#[test]
fn gamma_matches_ginac() {
if Command::new("ginsh").arg("--help").output().is_err() {
return;
}
let real_expected = ginsh_gamma("1.25", 50);
let real_actual = Complex::<Float>::try_from(parse!("gamma(5/4)").to_float(50)).unwrap();
assert_close_complex(&real_actual, &real_expected, "1e-40");
let complex_expected = ginsh_gamma("0.5+0.25*I", 50);
let complex_actual =
Complex::<Float>::try_from(parse!("gamma(1/2+1i/4)").to_float(50)).unwrap();
assert_close_complex(&complex_actual, &complex_expected, "1e-36");
}
#[test]
fn polygamma_order_zero_exact_normalization() {
assert_eq!(parse!("polygamma(0,1)"), parse!("-euler_gamma"));
assert_eq!(parse!("polygamma(0,1/2)"), parse!("-euler_gamma-2*log(2)"));
assert_eq!(
parse!("polygamma(0,5/2)"),
parse!("8/3-euler_gamma-2*log(2)")
);
assert_eq!(
parse!("polygamma(0,0)"),
Atom::num(Coefficient::complex_infinity())
);
}
#[test]
fn polygamma_order_zero_numeric_recurrence() {
let lhs = Complex::<Float>::try_from(
(parse!("polygamma(0,3/2+1i/4)") - parse!("polygamma(0,1/2+1i/4)")).to_float(50),
)
.unwrap();
let rhs = Complex::<Float>::try_from(parse!("1/(1/2+1i/4)").to_float(50)).unwrap();
assert_close_complex(&lhs, &rhs, "1e-30");
}
#[test]
fn special_float_inputs_normalize_immediately() {
assert!(Complex::<Float>::try_from(parse!("gamma(1.25)")).is_ok());
assert!(Complex::<Float>::try_from(parse!("erf(1.25)")).is_ok());
assert!(Complex::<Float>::try_from(parse!("polygamma(0,1.25)")).is_ok());
assert!(Complex::<Float>::try_from(parse!("polygamma(1,1.25)")).is_ok());
assert!(Complex::<Float>::try_from(parse!("zeta(1.25)")).is_ok());
}
#[test]
fn polygamma_basic_normalization() {
assert_eq!(
parse!("polygamma(1,0)"),
Atom::num(Coefficient::complex_infinity())
);
assert_eq!(
parse!("polygamma(0,x)").derivative(symbol!("x")),
parse!("polygamma(1,x)")
);
assert_eq!(
parse!("polygamma(3,x)").derivative(symbol!("x")),
parse!("polygamma(4,x)")
);
}
#[test]
fn polygamma_numeric_recurrence() {
let lhs = Complex::<Float>::try_from(
(parse!("polygamma(1,3/2+1i/4)") - parse!("polygamma(1,1/2+1i/4)")).to_float(50),
)
.unwrap();
let rhs = Complex::<Float>::try_from(parse!("-1/(1/2+1i/4)^2").to_float(50)).unwrap();
assert_close_complex(&lhs, &rhs, "1e-28");
}
#[test]
fn polygamma_laurent_series() {
let x = symbol!("x");
assert_eq!(
parse!("polygamma(0,x)").series(x, 0, 2).unwrap().to_atom(),
parse!("-x^-1-euler_gamma+1/6*pi^2*x+1/2*polygamma(2,1)*x^2")
);
assert_eq!(
parse!("polygamma(1,x)").series(x, 0, 2).unwrap().to_atom(),
parse!("x^-2+1/6*pi^2+polygamma(2,1)*x+1/30*pi^4*x^2")
);
}
#[test]
fn polylog_basic_normalization() {
assert_eq!(parse!("polylog(s,0)"), Atom::new());
assert_eq!(parse!("polylog(0,x)"), parse!("x/(1-x)"));
assert_eq!(parse!("polylog(-1,x)"), parse!("x/(1-x)^2"));
assert_eq!(parse!("polylog(-2,x)"), parse!("(x+x^2)/(1-x)^3"));
assert_eq!(parse!("polylog(1,x)"), parse!("-log(1-x)"));
assert_eq!(parse!("polylog(2,1)"), parse!("1/6*pi^2"));
assert_eq!(parse!("polylog(3,1)"), parse!("zeta(3)"));
assert_eq!(parse!("polylog(1,-1)"), parse!("-log(2)"));
assert_eq!(parse!("polylog(2,-1)"), parse!("-1/12*pi^2"));
assert_eq!(
parse!("polylog(3,x)").derivative(symbol!("x")),
parse!("polylog(2,x)/x")
);
}
#[test]
fn polylog_float_normalization() {
assert!(Complex::<Float>::try_from(parse!("polylog(2,0.5)")).is_ok());
}
#[test]
fn polylog_matches_ginac_for_dilog() {
if Command::new("ginsh").arg("--help").output().is_err() {
return;
}
let expected = ginsh_polylog("2", "0.5", 50);
let actual = Complex::<Float>::try_from(parse!("polylog(2,1/2)").to_float(50)).unwrap();
assert_close_complex(&actual, &expected, "1e-40");
}
#[test]
fn polylog_matches_ginac_off_series_region() {
if Command::new("ginsh").arg("--help").output().is_err() {
return;
}
let expected = ginsh_polylog("2", "2", 50);
let actual = Complex::<Float>::try_from(parse!("polylog(2,2)").to_float(50)).unwrap();
assert_close_complex(&actual, &expected, "1e-12");
let expected = ginsh_polylog("3", "11/10", 50);
let actual = Complex::<Float>::try_from(parse!("polylog(3,11/10)").to_float(50)).unwrap();
assert_close_complex(&actual, &expected, "1e-24");
}
#[test]
fn zeta_exact_normalization() {
assert_eq!(parse!("zeta(0)"), parse!("-1/2"));
assert_eq!(parse!("zeta(-1)"), parse!("-1/12"));
assert_eq!(parse!("zeta(-2)"), Atom::new());
assert_eq!(parse!("zeta(2)"), parse!("1/6*pi^2"));
assert_eq!(parse!("zeta(4)"), parse!("1/90*pi^4"));
assert_eq!(
parse!("zeta(1)"),
Atom::num(Coefficient::complex_infinity())
);
}
#[test]
fn zeta_float_normalization() {
assert!(Complex::<Float>::try_from(parse!("zeta(1.25)")).is_ok());
assert!(Complex::<Float>::try_from(parse!("zeta(1/2+1i/4)").to_float(50)).is_ok());
}
#[test]
fn special_functions_register_eval_info() {
let mut evaluator = parse!(
"gamma(x)+erf(x)+polygamma(0,x)+polygamma(1,x)+polylog(2,x)+zeta(x)+euler_gamma"
)
.evaluator(&[parse!("x")])
.build()
.unwrap()
.map_coeff(&|x| x.re.to_f64());
let mut out = [0.0];
evaluator.evaluate(&[0.25], &mut out);
assert!(out[0].is_finite());
}
#[test]
fn geometric_float_inputs_normalize_immediately() {
assert!(Complex::<Float>::try_from(parse!("tan(0.25)")).is_ok());
assert!(Complex::<Float>::try_from(parse!("atan(0.25)")).is_ok());
assert!(Complex::<Float>::try_from(parse!("atan(-1.0,1.0)")).is_ok());
assert!(Complex::<Float>::try_from(parse!("sinh(0.25)")).is_ok());
assert!(Complex::<Float>::try_from(parse!("atanh(0.25)")).is_ok());
assert!(Complex::<Float>::try_from(parse!("sec(0.25)")).is_ok());
assert!(Complex::<Float>::try_from(parse!("coth(1.25)")).is_ok());
}
#[test]
fn atan_two_argument_form_uses_quadrant() {
let actual = Complex::<Float>::try_from(parse!("atan(-1.0,1.0)")).unwrap();
let expected = Complex::new(Float::with_val(53, 1.0_f64.atan2(-1.0)), Float::new(53));
assert_close_complex(&actual, &expected, "1e-15");
let exact_actual = Complex::<Float>::try_from(parse!("atan(-1,1)").to_float(80)).unwrap();
let one = Float::with_val(80, 1);
let minus_one = Float::with_val(80, -1);
let exact_expected = Complex::new(one.atan2(&minus_one), Float::new(80));
assert_close_complex(&exact_actual, &exact_expected, "1e-23");
}
#[test]
fn transcendental_functions_extension_trait() {
use super::TranscendentalFunctions;
let x = parse!("x");
assert_eq!(x.atan(), parse!("atan(x)"));
assert_eq!(x.atan2(parse!("y")), parse!("atan(x,y)"));
assert_eq!(x.gamma(), parse!("gamma(x)"));
assert_eq!(x.erf(), parse!("erf(x)"));
assert_eq!(x.zeta(), parse!("zeta(x)"));
assert_eq!(x.polylog(2), parse!("polylog(2,x)"));
assert_eq!(x.polygamma(1), parse!("polygamma(1,x)"));
assert_eq!(x.bessel_j(parse!("nu")), parse!("bessel_j(nu,x)"));
}
#[test]
fn geometric_derivatives() {
assert_eq!(
parse!("tan(x)").derivative(symbol!("x")),
parse!("sec(x)^2")
);
assert_eq!(
parse!("atan(x)").derivative(symbol!("x")),
parse!("1/(1+x^2)")
);
assert_eq!(
parse!("atan(x,y)").derivative(symbol!("x")),
parse!("-y/(x^2+y^2)")
);
assert_eq!(
parse!("atan(x,y)").derivative(symbol!("y")),
parse!("x/(x^2+y^2)")
);
assert_eq!(
parse!("atan(y,x)").derivative(symbol!("x")),
parse!("y/(x^2+y^2)")
);
assert_eq!(
parse!("sech(x)").derivative(symbol!("x")),
parse!("-sech(x)*tanh(x)")
);
assert_eq!(
parse!("atanh(x)").derivative(symbol!("x")),
parse!("1/(1-x^2)")
);
}
#[test]
fn tan_laurent_series() {
let x = symbol!("x");
assert_eq!(
parse!("tan(x)")
.series(x, parse!("pi/2"), 0)
.unwrap()
.to_atom(),
parse!("-(x-pi/2)^-1")
);
assert_eq!(
parse!("sec(x)")
.series(x, parse!("pi/2"), 0)
.unwrap()
.to_atom(),
parse!("-(x-pi/2)^-1")
);
assert_eq!(
parse!("sec(x)")
.series(x, parse!("3*pi/2"), 0)
.unwrap()
.to_atom(),
parse!("(x-3/2*pi)^-1")
);
assert_eq!(
parse!("cot(x)").series(x, 0, 0).unwrap().to_atom(),
parse!("x^-1")
);
assert_eq!(
parse!("csc(x)").series(x, 0, 0).unwrap().to_atom(),
parse!("x^-1")
);
assert_eq!(
parse!("csc(x)")
.series(x, parse!("pi"), 0)
.unwrap()
.to_atom(),
parse!("-(x-pi)^-1")
);
}
#[test]
fn hyperbolic_laurent_series() {
let x = symbol!("x");
assert_eq!(
parse!("tanh(x)")
.series(x, parse!("1i*pi/2"), 0)
.unwrap()
.to_atom(),
parse!("(x-1i*pi/2)^-1")
);
assert_eq!(
parse!("coth(x)").series(x, 0, 0).unwrap().to_atom(),
parse!("x^-1")
);
assert_eq!(
parse!("sech(x)")
.series(x, parse!("1i*pi/2"), 0)
.unwrap()
.to_atom(),
parse!("-1i*(x-1i*pi/2)^-1")
);
assert_eq!(
parse!("sech(x)")
.series(x, parse!("3i*pi/2"), 0)
.unwrap()
.to_atom(),
parse!("1i*(x-3i*pi/2)^-1")
);
assert_eq!(
parse!("csch(x)").series(x, 0, 0).unwrap().to_atom(),
parse!("x^-1")
);
assert_eq!(
parse!("csch(x)")
.series(x, parse!("1i*pi"), 0)
.unwrap()
.to_atom(),
parse!("-(x-1i*pi)^-1")
);
}
#[test]
fn geometric_functions_register_eval_info() {
let mut evaluator = parse!("tan(x)+atan(x)+atan(x,-1)+sinh(x)+atanh(x/2)+sec(x)+coth(x+2)")
.evaluator(&[parse!("x")])
.build()
.unwrap()
.map_coeff(&|x| x.re.to_f64());
let mut out = [0.0];
evaluator.evaluate(&[0.25], &mut out);
assert!(out[0].is_finite());
}
#[test]
fn bessel_float_inputs_normalize_immediately() {
assert!(Complex::<Float>::try_from(parse!("bessel_j(1/2,2.0)")).is_ok());
assert!(Complex::<Float>::try_from(parse!("bessel_y(1/2,2.0)")).is_ok());
assert!(Complex::<Float>::try_from(parse!("bessel_i(1/2,2.0)")).is_ok());
assert!(Complex::<Float>::try_from(parse!("bessel_k(1/2,2.0)")).is_ok());
}
#[test]
fn bessel_known_half_integer_values() {
let j = Complex::<Float>::try_from(parse!("bessel_j(1/2,2)").to_float(80)).unwrap();
let y = Complex::<Float>::try_from(parse!("bessel_y(1/2,2)").to_float(80)).unwrap();
let i = Complex::<Float>::try_from(parse!("bessel_i(1/2,2)").to_float(80)).unwrap();
let k = Complex::<Float>::try_from(parse!("bessel_k(1/2,2)").to_float(80)).unwrap();
let j_expected =
Complex::<Float>::try_from(parse!("pi^(-1/2)*sin(2)").to_float(80)).unwrap();
let y_expected =
Complex::<Float>::try_from(parse!("-pi^(-1/2)*cos(2)").to_float(80)).unwrap();
let i_expected =
Complex::<Float>::try_from(parse!("pi^(-1/2)*sinh(2)").to_float(80)).unwrap();
let k_expected =
Complex::<Float>::try_from(parse!("1/2*pi^(1/2)*exp(-2)").to_float(80)).unwrap();
assert_close_complex(&j, &j_expected, "1e-20");
assert_close_complex(&y, &y_expected, "1e-20");
assert_close_complex(&i, &i_expected, "1e-20");
assert_close_complex(&k, &k_expected, "1e-18");
}
#[test]
fn bessel_derivatives() {
assert_eq!(
parse!("bessel_j(nu,x)").derivative(symbol!("x")),
parse!("(bessel_j(nu-1,x)-bessel_j(nu+1,x))/2")
);
assert_eq!(
parse!("bessel_i(nu,x)").derivative(symbol!("x")),
parse!("(bessel_i(nu-1,x)+bessel_i(nu+1,x))/2")
);
}
#[test]
fn bessel_functions_register_eval_info() {
let mut evaluator =
parse!("bessel_j(1/2,x)+bessel_y(1/2,x)+bessel_i(1/2,x)+bessel_k(1/2,x)")
.evaluator(&[parse!("x")])
.build()
.unwrap()
.map_coeff(&|x| x.re.to_f64());
let mut out = [0.0];
evaluator.evaluate(&[2.0], &mut out);
assert!(out[0].is_finite());
}
fn ginsh_gamma(argument: &str, digits: u32) -> Complex<Float> {
let script = format!("Digits={digits}:\nevalf(tgamma({argument}));\nquit;\n");
let mut child = Command::new("ginsh")
.stdin(Stdio::piped())
.stdout(Stdio::piped())
.spawn()
.expect("ginsh must be available");
child
.stdin
.as_mut()
.unwrap()
.write_all(script.as_bytes())
.expect("ginsh stdin must be writable");
let output = child.wait_with_output().expect("ginsh must complete");
assert!(
output.status.success(),
"ginsh failed: {}",
String::from_utf8_lossy(&output.stderr)
);
let stdout = String::from_utf8(output.stdout).unwrap();
let value = stdout
.lines()
.map(str::trim)
.find(|line| !line.is_empty())
.unwrap();
parse_ginsh_complex(value)
}
fn ginsh_polylog(order: &str, argument: &str, digits: u32) -> Complex<Float> {
let script = format!("Digits={digits}:\nevalf(Li({order},{argument}));\nquit;\n");
let mut child = Command::new("ginsh")
.stdin(Stdio::piped())
.stdout(Stdio::piped())
.spawn()
.expect("ginsh must be available");
child
.stdin
.as_mut()
.unwrap()
.write_all(script.as_bytes())
.expect("ginsh stdin must be writable");
let output = child.wait_with_output().expect("ginsh must complete");
assert!(
output.status.success(),
"ginsh failed: {}",
String::from_utf8_lossy(&output.stderr)
);
let stdout = String::from_utf8(output.stdout).unwrap();
let value = stdout
.lines()
.map(str::trim)
.find(|line| !line.is_empty())
.unwrap();
parse_ginsh_complex(value)
}
fn assert_close_complex(actual: &Complex<Float>, expected: &Complex<Float>, tolerance: &str) {
let tol = Float::parse(tolerance, Some(256)).unwrap();
let diff = (actual.clone() - expected.clone()).norm().re;
let scale = expected.norm().re;
let limit = if scale.to_f64() == 0.0 {
tol
} else {
tol * scale
};
assert!(
diff <= limit,
"difference too large: actual={}, expected={}, diff={}, limit={}",
actual,
expected,
diff,
limit
);
}
fn parse_ginsh_complex(value: &str) -> Complex<Float> {
let value = value.trim();
if !value.contains("*I") {
return Complex::new(Float::parse(value, Some(256)).unwrap(), Float::new(256));
}
let imag_marker = value
.rfind("*I")
.expect("GiNaC complex output must contain *I");
let core = &value[..imag_marker];
let split = core
.char_indices()
.skip(1)
.filter(|(_, ch)| *ch == '+' || *ch == '-')
.filter(|(idx, _)| !matches!(core.as_bytes()[idx - 1], b'e' | b'E'))
.map(|(idx, _)| idx)
.last()
.expect("GiNaC complex output must contain a real and imaginary part");
let re = Float::parse(&core[..split], Some(256)).unwrap();
let im = Float::parse(&core[split..], Some(256)).unwrap();
Complex::new(re, im)
}
}