use std::sync::LazyLock;

use rug::float::Constant;

use crate::{
    atom::{Atom, AtomCore, AtomOrView, AtomView, EvaluationInfo, FunctionBuilder, Symbol},
    coefficient::{Coefficient, CoefficientView},
    domains::{
        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"),
    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,
    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: &[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, rug::float::Special::Nan),
                                    Float::with_val(53, rug::float::Special::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: &[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)
                        && let Some(order) = atom_to_integer(s)
                    {
                        if order == 1 {
                            **out = Atom::num(1) / (Atom::num(1) - z);
                        } else {
                            **out = function!(polylog_symbol, Atom::num(order - 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, rug::float::Special::Nan),
                                    Float::with_val(53, rug::float::Special::Nan),
                                );
                            };
                            let [arg] = args else {
                                return Complex::new(
                                    Float::with_val(53, rug::float::Special::Nan),
                                    Float::with_val(53, rug::float::Special::Nan),
                                );
                            };
                            polylog_numeric_eval(&order, arg, prec).unwrap_or_else(|| {
                                Complex::new(
                                    Float::with_val(53, rug::float::Special::Nan),
                                    Float::with_val(53, rug::float::Special::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;
                        }

                        // TODO: add Stieltjes constants
                        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: &[f64]| zeta_eval_f64(args))
                    .register(|args: &[Complex<f64>]| zeta_eval_complex_f64(args)),
                )
            }
        );

        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,
            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: &[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: &[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: &[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: &[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: &[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: &[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: &[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: &[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: &[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: &[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: &[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: &[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, rug::float::Special::Nan),
                        Float::with_val(53, rug::float::Special::Nan),
                    ),
                })
                .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: &[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: &[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: &[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: &[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: &[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: &[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: &[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: &[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: &[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 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 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 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 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"])
}

/// Return the built-in Euler gamma function symbol `gamma`.
///
/// `gamma(z)` is meromorphic with simple poles at the non-positive integers and no branch cuts.
pub fn gamma() -> Symbol {
    SPECIALS.gamma
}

/// Return the built-in Euler-Mascheroni constant symbol `euler_gamma`.
///
/// This is the constant `gamma_E`, approximately `0.57721`, used in expansions of `gamma`,
/// `polygamma`, and related special functions.
pub fn euler_gamma() -> Symbol {
    SPECIALS.euler_gamma
}

/// Return the built-in polygamma function symbol `polygamma`.
///
/// The first argument is the non-negative integer order. For fixed order, `polygamma(n, z)` is
/// meromorphic in `z` with poles at the non-positive integers and no branch cuts.
pub fn polygamma() -> Symbol {
    SPECIALS.polygamma
}

/// Return the built-in polylogarithm symbol `polylog`.
///
/// `polylog(s, z)` uses the principal branch in `z`; for generic `s` it has the standard branch cut
/// along the real interval `[1, +infinity)`.
pub fn polylog() -> Symbol {
    SPECIALS.polylog
}

/// Return the built-in Riemann zeta function symbol `zeta`.
///
/// `zeta(s)` is meromorphic with a simple pole at `s = 1` and no branch cuts.
pub fn zeta() -> Symbol {
    SPECIALS.zeta
}

/// Return the built-in tangent function symbol `tan`.
///
/// `tan(z)` is meromorphic with simple poles at `pi/2 + k pi`.
pub fn tan() -> Symbol {
    GEOMETRICS.tan
}

/// Return the built-in cotangent function symbol `cot`.
///
/// `cot(z)` is meromorphic with simple poles at `k pi`.
pub fn cot() -> Symbol {
    GEOMETRICS.cot
}

/// Return the built-in secant function symbol `sec`.
///
/// `sec(z)` is meromorphic with simple poles at `pi/2 + k pi`.
pub fn sec() -> Symbol {
    GEOMETRICS.sec
}

/// Return the built-in cosecant function symbol `csc`.
///
/// `csc(z)` is meromorphic with simple poles at `k pi`.
pub fn csc() -> Symbol {
    GEOMETRICS.csc
}

/// Return the built-in inverse sine symbol `asin`.
///
/// `asin(z)` uses the principal branch with branch cuts on `(-infinity, -1]` and `[1, +infinity)`.
pub fn asin() -> Symbol {
    GEOMETRICS.asin
}

/// Return the built-in inverse cosine symbol `acos`.
///
/// `acos(z)` uses the principal branch with branch cuts on `(-infinity, -1]` and `[1, +infinity)`.
pub fn acos() -> Symbol {
    GEOMETRICS.acos
}

/// Return the built-in inverse tangent symbol `atan`.
///
/// `atan(z)` uses the principal branch with branch cuts on the imaginary axis:
/// `(-i infinity, -i]` and `[i, i infinity)`.
///
/// `atan(x, y)` is the quadrant-aware inverse tangent for real numeric inputs,
/// equivalent to `atan2(y, x)`.
pub fn atan() -> Symbol {
    GEOMETRICS.atan
}

/// Return the built-in inverse cotangent symbol `acot`.
///
/// `acot(z)` uses the principal branch with branch cuts on the imaginary axis:
/// `(-i infinity, -i]` and `[i, i infinity)`.
pub fn acot() -> Symbol {
    GEOMETRICS.acot
}

/// Return the built-in inverse secant symbol `asec`.
///
/// `asec(z)` uses the principal branch with branch cuts on `[-1, 1]`.
pub fn asec() -> Symbol {
    GEOMETRICS.asec
}

/// Return the built-in inverse cosecant symbol `acsc`.
///
/// `acsc(z)` uses the principal branch with branch cuts on `[-1, 1]`.
pub fn acsc() -> Symbol {
    GEOMETRICS.acsc
}

/// Return the built-in hyperbolic sine symbol `sinh`.
///
/// `sinh(z)` is entire.
pub fn sinh() -> Symbol {
    GEOMETRICS.sinh
}

/// Return the built-in hyperbolic cosine symbol `cosh`.
///
/// `cosh(z)` is entire.
pub fn cosh() -> Symbol {
    GEOMETRICS.cosh
}

/// Return the built-in hyperbolic tangent symbol `tanh`.
///
/// `tanh(z)` is meromorphic with simple poles at `i (pi/2 + k pi)`.
pub fn tanh() -> Symbol {
    GEOMETRICS.tanh
}

/// Return the built-in hyperbolic cotangent symbol `coth`.
///
/// `coth(z)` is meromorphic with simple poles at `i k pi`.
pub fn coth() -> Symbol {
    GEOMETRICS.coth
}

/// Return the built-in hyperbolic secant symbol `sech`.
///
/// `sech(z)` is meromorphic with simple poles at `i (pi/2 + k pi)`.
pub fn sech() -> Symbol {
    GEOMETRICS.sech
}

/// Return the built-in hyperbolic cosecant symbol `csch`.
///
/// `csch(z)` is meromorphic with simple poles at `i k pi`.
pub fn csch() -> Symbol {
    GEOMETRICS.csch
}

/// Return the built-in inverse hyperbolic sine symbol `asinh`.
///
/// `asinh(z)` uses the principal branch with branch cuts on the imaginary axis:
/// `(-i infinity, -i]` and `[i, i infinity)`.
pub fn asinh() -> Symbol {
    GEOMETRICS.asinh
}

/// Return the built-in inverse hyperbolic cosine symbol `acosh`.
///
/// `acosh(z)` uses the principal branch with branch cut `(-infinity, 1]`.
pub fn acosh() -> Symbol {
    GEOMETRICS.acosh
}

/// Return the built-in inverse hyperbolic tangent symbol `atanh`.
///
/// `atanh(z)` uses the principal branch with branch cuts on `(-infinity, -1]` and `[1, +infinity)`.
pub fn atanh() -> Symbol {
    GEOMETRICS.atanh
}

/// Return the built-in inverse hyperbolic cotangent symbol `acoth`.
///
/// `acoth(z)` uses the principal branch with branch cut `[-1, 1]`.
pub fn acoth() -> Symbol {
    GEOMETRICS.acoth
}

/// Return the built-in inverse hyperbolic secant symbol `asech`.
///
/// `asech(z)` uses the principal branch with branch cuts on `(-infinity, 0]` and `[1, +infinity)`.
pub fn asech() -> Symbol {
    GEOMETRICS.asech
}

/// Return the built-in inverse hyperbolic cosecant symbol `acsch`.
///
/// `acsch(z)` uses the principal branch with branch cuts on the imaginary interval `[-i, i]`.
pub fn acsch() -> Symbol {
    GEOMETRICS.acsch
}

/// Return the built-in cylindrical Bessel function of the first kind symbol `bessel_j`.
///
/// For fixed order `nu`, `bessel_j(nu, z)` is entire in `z`.
pub fn bessel_j() -> Symbol {
    BESSELS.bessel_j
}

/// Return the built-in cylindrical Bessel function of the second kind symbol `bessel_y`.
///
/// For fixed order `nu`, `bessel_y(nu, z)` uses the principal branch in `z` with branch cut
/// on `(-infinity, 0]`.
pub fn bessel_y() -> Symbol {
    BESSELS.bessel_y
}

/// Return the built-in modified Bessel function of the first kind symbol `bessel_i`.
///
/// For fixed order `nu`, `bessel_i(nu, z)` is entire in `z`.
pub fn bessel_i() -> Symbol {
    BESSELS.bessel_i
}

/// Return the built-in modified Bessel function of the second kind symbol `bessel_k`.
///
/// For fixed order `nu`, `bessel_k(nu, z)` uses the principal branch in `z` with branch cut
/// on `(-infinity, 0]`.
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())
            }
        )*
    };
}

/// Convenience methods for transcendental functions on expressions.
///
/// This extension trait is implemented for every type that implements [`AtomCore`].
/// Import it to write `x.atan()`, `x.gamma()`, `x.polylog(2)`, and similar method
/// calls without constructing the function symbols manually.
///
/// Functions already provided by [`AtomCore`] itself, such as `exp`, `log`, `sin`,
/// `cos`, and `sqrt`, are not repeated here to avoid method-resolution ambiguity.
pub trait TranscendentalFunctions: AtomCore {
    unary_transcendental_methods! {
        /// Take the tangent.
        tan => tan;
        /// Take the cotangent.
        cot => cot;
        /// Take the secant.
        sec => sec;
        /// Take the cosecant.
        csc => csc;
        /// Take the inverse sine.
        asin => asin;
        /// Take the inverse cosine.
        acos => acos;
        /// Take the inverse tangent.
        atan => atan;
        /// Take the inverse cotangent.
        acot => acot;
        /// Take the inverse secant.
        asec => asec;
        /// Take the inverse cosecant.
        acsc => acsc;
        /// Take the hyperbolic sine.
        sinh => sinh;
        /// Take the hyperbolic cosine.
        cosh => cosh;
        /// Take the hyperbolic tangent.
        tanh => tanh;
        /// Take the hyperbolic cotangent.
        coth => coth;
        /// Take the hyperbolic secant.
        sech => sech;
        /// Take the hyperbolic cosecant.
        csch => csch;
        /// Take the inverse hyperbolic sine.
        asinh => asinh;
        /// Take the inverse hyperbolic cosine.
        acosh => acosh;
        /// Take the inverse hyperbolic tangent.
        atanh => atanh;
        /// Take the inverse hyperbolic cotangent.
        acoth => acoth;
        /// Take the inverse hyperbolic secant.
        asech => asech;
        /// Take the inverse hyperbolic cosecant.
        acsch => acsch;
        /// Apply the gamma function.
        gamma => gamma;
        /// Apply the Riemann zeta function.
        zeta => zeta;
    }

    /// Compute the quadrant-aware inverse tangent `atan(self, y)`.
    ///
    /// For real numeric inputs, this is equivalent to `atan2(y, self)`.
    fn atan2<'a, T: Into<AtomOrView<'a>>>(&self, y: T) -> Self::Output {
        receiver_first_transcendental_function(self, crate::transcendental::atan(), y)
    }

    /// Compute the polygamma function of order `n` at `self`.
    fn polygamma<'a, T: Into<AtomOrView<'a>>>(&self, n: T) -> Self::Output {
        receiver_last_transcendental_function(self, crate::transcendental::polygamma(), n)
    }

    /// Compute the polylogarithm of order `s` at `self`.
    fn polylog<'a, T: Into<AtomOrView<'a>>>(&self, s: T) -> Self::Output {
        receiver_last_transcendental_function(self, crate::transcendental::polylog(), s)
    }

    /// Compute the cylindrical Bessel function of the first kind of order `nu` at `self`.
    fn bessel_j<'a, T: Into<AtomOrView<'a>>>(&self, nu: T) -> Self::Output {
        receiver_last_transcendental_function(self, crate::transcendental::bessel_j(), nu)
    }

    /// Compute the cylindrical Bessel function of the second kind of order `nu` at `self`.
    fn bessel_y<'a, T: Into<AtomOrView<'a>>>(&self, nu: T) -> Self::Output {
        receiver_last_transcendental_function(self, crate::transcendental::bessel_y(), nu)
    }

    /// Compute the modified Bessel function of the first kind of order `nu` at `self`.
    fn bessel_i<'a, T: Into<AtomOrView<'a>>>(&self, nu: T) -> Self::Output {
        receiver_last_transcendental_function(self, crate::transcendental::bessel_i(), nu)
    }

    /// Compute the modified Bessel function of the second kind of order `nu` at `self`.
    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, rug::float::Special::Nan),
            Float::with_val(53, rug::float::Special::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, rug::float::Special::Nan),
            Float::with_val(53, rug::float::Special::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 *= &current;
        current += &one;
    }

    while current > target {
        current -= &one;
        prefactor /= &current;
    }

    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 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> {
    if z.im.to_f64() == 0.0 {
        Complex::new(
            z.re.clone().into_inner().gamma().into(),
            Float::new(binary_prec),
        )
    } else {
        gamma_complex_spouge(z, binary_prec)
    }
}

fn polygamma_order_zero_numeric_eval(z: &Complex<Float>, binary_prec: u32) -> Complex<Float> {
    if z.im.to_f64() == 0.0 {
        Complex::new(
            z.re.clone().into_inner().digamma().into(),
            Float::new(binary_prec),
        )
    } else {
        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());
    }

    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, rug::float::Special::Infinity),
                zero,
            );
        }

        return Complex::new(s.re.clone().into_inner().zeta().into(), zero);
    }

    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, rug::float::Special::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);
    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)
    };

    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 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 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!("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)+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.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)
    }
}
