symbolica 2.2.0

use std::{
    ops::{Add, DerefMut, Div, Mul, Sub},
    sync::Arc,
};

use numerica::printer::PrintOptions;

use crate::{
    atom::{Atom, AtomCore, AtomView, FunctionBuilder, Indeterminate, InlineVar, Symbol},
    coefficient::{Coefficient, CoefficientView},
    domains::{Ring, atom::AtomField, integer::Integer, rational::Rational},
    error,
    poly::{
        PolyVariable,
        series::{Series, SeriesDepth, SeriesError},
    },
    state::Workspace,
};

impl AtomView<'_> {
    /// Take a derivative of the expression with respect to `x`.
    pub(crate) fn derivative(&self, x: &Indeterminate) -> Atom {
        Workspace::get_local().with(|ws| {
            let mut out = ws.new_atom();
            self.derivative_with_ws_into(x, ws, &mut out);
            out.into_inner()
        })
    }

    /// Take a derivative of the expression with respect to `x` and
    /// write the result in `out`.
    /// Returns `true` if the derivative is non-zero.
    pub(crate) fn derivative_into(&self, x: &Indeterminate, out: &mut Atom) -> bool {
        Workspace::get_local().with(|ws| self.derivative_with_ws_into(x, ws, out))
    }

    /// Take a derivative of the expression with respect to `x` and
    /// write the result in `out`.
    /// Returns `true` if the derivative is non-zero.
    pub(crate) fn derivative_with_ws_into(
        &self,
        x: &Indeterminate,
        workspace: &Workspace,
        out: &mut Atom,
    ) -> bool {
        if x == self {
            out.to_num(1);
            return true;
        }

        match self {
            AtomView::Num(_) | AtomView::Var(_) => {
                out.to_num(Coefficient::zero());
                false
            }
            AtomView::Fun(f_orig) => {
                // detect if the function to derive is the derivative function itself.
                // der stores: der(depth_1, ..., depth_n, f, arg_1, ..., arg_n).
                let mut der_function_call = workspace.new_atom();
                let (f, is_der, arg_count) = if f_orig.get_symbol() == Symbol::DERIVATIVE {
                    if f_orig.get_nargs() < 3 || f_orig.get_nargs() % 2 == 0 {
                        error!(
                            "Derivative function {} must contain n depths, a function symbol, and n arguments",
                            self.printer(PrintOptions::file())
                        );
                        out.to_num(Coefficient::Indeterminate);
                        return false;
                    }

                    let der_depth_count = (f_orig.get_nargs() - 1) / 2;
                    let function_symbol = match f_orig.iter().nth(der_depth_count).unwrap() {
                        AtomView::Var(v) => v.get_symbol(),
                        _ => {
                            error!(
                                "Derivative function {} argument must be a function symbol",
                                self.printer(PrintOptions::file())
                            );
                            out.to_num(Coefficient::Indeterminate);
                            return false;
                        }
                    };

                    let function_call = der_function_call.to_fun(function_symbol);
                    for arg in f_orig.iter().skip(der_depth_count + 1) {
                        function_call.add_arg(arg);
                    }

                    (
                        match der_function_call.as_view() {
                            AtomView::Fun(f) => f,
                            _ => unreachable!(),
                        },
                        true,
                        der_depth_count,
                    )
                } else {
                    (*f_orig, false, f_orig.get_nargs())
                };

                // take derivative of all the arguments and store it in a list
                let mut args_der = Vec::with_capacity(f.get_nargs());
                for (i, arg) in f.iter().enumerate() {
                    let mut arg_der = workspace.new_atom();
                    if arg.derivative_with_ws_into(x, workspace, &mut arg_der) {
                        args_der.push((i, arg_der));
                    }
                }

                if args_der.is_empty() {
                    out.to_num(Coefficient::zero());
                    return false;
                }

                // derive special functions
                if f.get_nargs() == 1
                    && [Symbol::LOG, Symbol::SIN, Symbol::COS].contains(&f.get_symbol())
                {
                    let mut fn_der = workspace.new_atom();
                    match f.get_symbol_id() {
                        Symbol::LOG_ID => {
                            let mut n = workspace.new_atom();
                            n.to_num(-1);

                            fn_der.to_pow(f.iter().next().unwrap(), n.as_view());
                        }
                        Symbol::SIN_ID => {
                            let p = fn_der.to_fun(Symbol::COS);
                            p.add_arg(f.iter().next().unwrap());
                        }
                        Symbol::COS_ID => {
                            let mut n = workspace.new_atom();
                            n.to_num(-1);

                            let mut sin = workspace.new_atom();
                            let sin_fun = sin.to_fun(Symbol::SIN);
                            sin_fun.add_arg(f.iter().next().unwrap());

                            let m = fn_der.to_mul();
                            m.extend(sin.as_view());
                            m.extend(n.as_view());
                        }
                        _ => unreachable!(),
                    }

                    let (_, mut arg_der) = args_der.pop().unwrap_or((0, Atom::Zero.into()));
                    if let Atom::Mul(m) = arg_der.deref_mut() {
                        m.extend(fn_der.as_view());
                        arg_der.as_view().normalize(workspace, out);
                    } else {
                        let mut mul = workspace.new_atom();
                        let m = mul.to_mul();
                        m.extend(fn_der.as_view());
                        m.extend(arg_der.as_view());
                        mul.as_view().normalize(workspace, out);
                    }

                    return true;
                }

                if f.get_nargs() == 3 && f.get_symbol_id() == Symbol::IF_ID {
                    let mut fn_der = workspace.new_atom();
                    let der = fn_der.to_fun(Symbol::IF);
                    der.add_arg(f.iter().next().unwrap());

                    if let Some((_, d)) = args_der.iter().find(|(i, _)| *i == 1) {
                        der.add_arg(d.as_view());
                    } else {
                        der.add_arg(Atom::Zero.as_view());
                    };

                    if let Some((_, d)) = args_der.iter().find(|(i, _)| *i == 2) {
                        der.add_arg(d.as_view());
                    } else {
                        der.add_arg(Atom::Zero.as_view());
                    };

                    fn_der.as_view().normalize(workspace, out);
                    return true;
                }

                // create a derivative function that tags which index was derived
                let mut add = workspace.new_atom();
                let a = add.to_add();
                let mut fn_der = workspace.new_atom();
                let mut n = workspace.new_atom();
                let mut mul = workspace.new_atom();
                for (index, arg_der) in args_der {
                    if let Some(custom_der) = &f.get_symbol().get_global_data().custom_derivative {
                        let mut setter = fn_der.deref_mut().into();
                        custom_der(*self, index, &mut setter);
                        if setter.is_set() {
                            let m = mul.to_mul();
                            m.extend(fn_der.as_view());
                            m.extend(arg_der.as_view());
                            a.extend(m.as_view());
                            continue;
                        }
                    }

                    let p = fn_der.to_fun(Symbol::DERIVATIVE);

                    if is_der {
                        for (i, x_orig) in f_orig.iter().take(arg_count).enumerate() {
                            n.set_from_view(&x_orig);
                            *n += if i == index { 1 } else { 0 };
                            p.add_arg(n.as_view());
                        }
                    } else {
                        for i in 0..f.get_nargs() {
                            n.to_num((if i == index { 1 } else { 0 }, 1));
                            p.add_arg(n.as_view());
                        }
                    }

                    let function_symbol = InlineVar::new(f.get_symbol());
                    p.add_arg(function_symbol.as_view());
                    for arg in f.iter() {
                        p.add_arg(arg);
                    }

                    let m = mul.to_mul();
                    m.extend(fn_der.as_view());
                    m.extend(arg_der.as_view());
                    mul.as_view().normalize(workspace, out);

                    a.extend(mul.as_view());
                }

                add.as_view().normalize(workspace, out);
                true
            }
            AtomView::Pow(p) => {
                let (base, exp) = p.get_base_exp();

                let mut exp_der = workspace.new_atom();
                let exp_der_non_zero = exp.derivative_with_ws_into(x, workspace, &mut exp_der);

                let mut base_der = workspace.new_atom();
                let base_der_non_zero = base.derivative_with_ws_into(x, workspace, &mut base_der);

                if !exp_der_non_zero && !base_der_non_zero {
                    out.to_num(0);
                    return false;
                }

                let mut exp_der_contrib = workspace.new_atom();

                if exp_der_non_zero {
                    // create log(base)
                    let mut log_base = workspace.new_atom();
                    let lb = log_base.to_fun(Symbol::LOG);
                    lb.add_arg(base);

                    if let Atom::Mul(m) = exp_der.deref_mut() {
                        m.extend(*self);
                        m.extend(log_base.as_view());
                        exp_der.as_view().normalize(workspace, &mut exp_der_contrib);
                    } else {
                        let mut mul = workspace.new_atom();
                        let m = mul.to_mul();
                        m.extend(*self);
                        m.extend(exp_der.as_view());
                        m.extend(log_base.as_view());
                        mul.as_view().normalize(workspace, &mut exp_der_contrib);
                    }

                    if !base_der_non_zero {
                        out.set_from_view(&exp_der_contrib.as_view());
                        return true;
                    }
                }

                let mut mul_h = workspace.new_atom();
                let mul = mul_h.to_mul();
                mul.extend(base_der.as_view());

                let mut new_exp = workspace.new_atom();
                if let AtomView::Num(n) = exp {
                    mul.extend(exp);

                    let res = n.get_coeff_view() + -1;
                    new_exp.to_num(res);
                } else {
                    mul.extend(exp);

                    let ao = new_exp.to_add();
                    ao.extend(exp);

                    let mut min_one = workspace.new_atom();
                    min_one.to_num(-1);

                    ao.extend(min_one.as_view());
                }

                let mut pow_h = workspace.new_atom();
                pow_h.to_pow(base, new_exp.as_view());

                mul.extend(pow_h.as_view());

                if exp_der_non_zero {
                    let mut add = workspace.new_atom();
                    let a = add.to_add();

                    a.extend(mul_h.as_view());
                    a.extend(exp_der_contrib.as_view());

                    add.as_view().normalize(workspace, out);
                } else {
                    mul_h.as_view().normalize(workspace, out);
                }

                true
            }
            AtomView::Mul(args) => {
                let mut add_h = workspace.new_atom();
                let add = add_h.to_add();
                let mut mul_h = workspace.new_atom();
                let mut non_zero = false;
                for arg in args.iter() {
                    let mut arg_der = workspace.new_atom();
                    if arg.derivative_with_ws_into(x, workspace, &mut arg_der) {
                        if let Atom::Mul(mm) = arg_der.deref_mut() {
                            for other_arg in args.iter() {
                                if other_arg != arg {
                                    mm.extend(other_arg);
                                }
                            }

                            add.extend(arg_der.as_view());
                        } else {
                            let mm = mul_h.to_mul();
                            mm.extend(arg_der.as_view());
                            for other_arg in args.iter() {
                                if other_arg != arg {
                                    mm.extend(other_arg);
                                }
                            }
                            add.extend(mul_h.as_view());
                        }

                        non_zero = true;
                    }
                }

                if non_zero {
                    add_h.as_view().normalize(workspace, out);
                    true
                } else {
                    out.to_num(0);
                    false
                }
            }
            AtomView::Add(args) => {
                let mut add_h = workspace.new_atom();
                let add = add_h.to_add();
                let mut arg_der = workspace.new_atom();
                let mut non_zero = false;
                for arg in args.iter() {
                    if arg.derivative_with_ws_into(x, workspace, &mut arg_der) {
                        add.extend(arg_der.as_view());
                        non_zero = true;
                    }
                }

                if non_zero {
                    add_h.as_view().normalize(workspace, out);
                    true
                } else {
                    out.to_num(0);
                    false
                }
            }
        }
    }

    pub(crate) fn series(
        &self,
        x: &Indeterminate,
        expansion_point: AtomView,
        depth: SeriesDepth,
    ) -> Result<Series<AtomField>, SeriesError> {
        let order = depth.order();
        if !depth.is_absolute() && (order.is_negative() || order.is_zero()) {
            return Err(SeriesError::NonPositiveRelativeDepth { depth });
        }

        // heuristic current depth
        let mut current_depth = if order.is_negative() || order.is_zero() {
            Rational::one()
        } else {
            order.clone()
        };

        // do not do an expensive statistical zero check at all stages during the series expansion
        // note that a statistical check will be used in the case of negative exponents in the expansion
        let field = AtomField {
            statistical_zero_test: false,
            ..Default::default()
        };

        loop {
            let info = Series::new(
                &field,
                None,
                Arc::new(x.clone().into()),
                expansion_point.to_owned(),
                &current_depth + &(1.into(), current_depth.denominator()).into(),
            );

            let mut series = self.series_impl(x, expansion_point, &info)?;
            if !depth.is_absolute() && series.relative_order() >= *order {
                series.truncate_relative_order(order.clone());
                break Ok(series);
            } else if depth.is_absolute() && series.absolute_order() > *order {
                series.truncate_absolute_order(order + &(1.into(), order.denominator()).into());
                break Ok(series);
            } else {
                // increase the expansion depth
                // TODO: find better heuristic
                current_depth = &current_depth * &2.into();
            }
        }
    }

    /// Series expand in `x` around `expansion_point` to depth `depth`.
    pub(crate) fn series_impl(
        &self,
        x: &Indeterminate,
        expansion_point: AtomView,
        info: &Series<AtomField>,
    ) -> Result<Series<AtomField>, SeriesError> {
        if !self.contains_indeterminate(x) {
            return Ok(info.constant(self.to_owned()));
        }

        if *x == *self {
            return Ok(info.shifted_variable(expansion_point.to_owned()));
        }

        // TODO: optimize, appending a monomial using addition is slow
        match self {
            AtomView::Num(_) | AtomView::Var(_) => Ok(info.constant(self.to_owned())),
            AtomView::Fun(f) => {
                let mut args_series = Vec::with_capacity(f.get_nargs());
                for arg in f {
                    args_series.push(arg.series_impl(x, expansion_point, info)?);
                }

                if args_series.is_empty() {
                    return Ok(info.constant(f.to_owned().into()));
                }

                let symbol = f.get_symbol();

                if symbol == Symbol::IF && f.iter().skip(1).any(|arg| arg.contains_indeterminate(x))
                {
                    return Err(SeriesError::NonConstantIf {
                        expression: self.to_owned(),
                    });
                }

                if !f.get_symbol().is_fixed_builtin()
                    && args_series
                        .iter()
                        .any(|x| x.get_trailing_exponent().is_negative())
                {
                    // fill in the expanded arguments, perhaps the leading negative exponent will be popped out,
                    // in which case we can proceed
                    let mut f_eval = FunctionBuilder::new(f.get_symbol());
                    for c in &args_series {
                        f_eval = f_eval.add_arg(c.to_atom());
                    }
                    let a = f_eval.finish();

                    if !matches!(a, Atom::Fun(_)) {
                        return a.as_view().series_impl(x, expansion_point, info);
                    }
                }

                match f.get_symbol_id() {
                    Symbol::COS_ID => args_series[0].cos(),
                    Symbol::SIN_ID => args_series[0].sin(),
                    Symbol::LOG_ID => args_series[0].log(),
                    _ => {
                        if let Some(custom_series) = &f.get_symbol().get_series_function()
                            && let Some((singular, regularized)) = custom_series(&args_series)
                        {
                            let singular_series =
                                singular.as_view().series_impl(x, expansion_point, info)?;
                            // TODO: expand deeper?
                            let regularized_series =
                                regularized
                                    .as_view()
                                    .series_impl(x, expansion_point, info)?;
                            return Ok(&singular_series * &regularized_series);
                        }

                        // TODO: also check for log(x)
                        if args_series
                            .iter()
                            .any(|x| x.get_trailing_exponent().is_negative())
                        {
                            return Err(SeriesError::FunctionArgumentPole {
                                expression: self.to_owned(),
                            });
                        }

                        let mut f_eval = FunctionBuilder::new(f.get_symbol());
                        for c in &args_series {
                            f_eval = f_eval.add_arg(c.to_atom());
                        }
                        let a = f_eval.finish();

                        let constant = a.replace(x.clone()).with(expansion_point.to_owned());

                        // TODO: depth is an overestimate
                        let order = info.absolute_order();
                        let depth = order.numerator().to_i64().unwrap() as u32
                            * order.denominator().to_i64().unwrap() as u32;

                        let mut result = info.constant(constant.clone());

                        // TODO: do not allow series expansion in IF if it is not a constant!

                        let mut d = a.clone();
                        for i in 1..=depth {
                            d = d.as_view().derivative(x);

                            if d.is_zero() {
                                break;
                            }

                            let rep = d
                                .replace(x.clone())
                                .with(expansion_point.to_owned())
                                .expand();

                            result = &result
                                + &info
                                    .monomial(info.get_field().one(), i.into())
                                    .mul_coeff(&rep)
                                    .div_coeff(&Atom::num(Integer::factorial(i)));
                        }

                        Ok(result)
                    }
                }
            }
            AtomView::Pow(p) => {
                let (base, exp) = p.get_base_exp();

                let mut base_series = base.series_impl(x, expansion_point, info)?;

                if let AtomView::Num(n) = exp {
                    if let CoefficientView::Natural(n, d, ni, _) = n.get_coeff_view() {
                        if ni != 0 {
                            return Err(SeriesError::UnsupportedComplexExponent {
                                expression: self.to_owned(),
                            });
                        }

                        if n < 0 {
                            // in case of 1/0, grow the expansion depth of the base series
                            // it could be that the base series is exactly zero,
                            // to prevent an infinite loop, we stop the loop at ep^-1000
                            let mut current_depth = info.relative_order();
                            let mut first_attempt = true;

                            // make sure the trailing coefficient is not zero with a statistical test
                            while first_attempt
                                && !base_series
                                    .get_trailing_coefficient()
                                    .zero_test(10, 1e-5)
                                    .is_false()
                                || !first_attempt && base_series.is_zero() && current_depth < 1000
                            {
                                let info = Series::new(
                                    &AtomField {
                                        statistical_zero_test: true,
                                        ..info.get_field().clone()
                                    },
                                    None,
                                    info.get_variable().clone(),
                                    info.get_expansion_point().clone(),
                                    &current_depth
                                        + &(1.into(), current_depth.denominator()).into(),
                                );

                                base_series = base.series_impl(x, expansion_point, &info)?;
                                current_depth = &current_depth * &2.into();
                                first_attempt = false;
                            }
                        }

                        base_series.rpow((n, d).into())
                    } else {
                        Err(SeriesError::UnsupportedLargeExponent {
                            expression: self.to_owned(),
                        })
                    }
                } else {
                    let e = exp.series_impl(x, expansion_point, info)?;
                    base_series.pow(&e)
                }
            }
            AtomView::Mul(args) => {
                let mut iter = args.iter();
                let mut series = iter.next().unwrap().series_impl(x, expansion_point, info)?;
                for arg in iter {
                    series = &series * &arg.series_impl(x, expansion_point, info)?;
                }

                Ok(series)
            }
            AtomView::Add(args) => {
                let mut iter = args.iter();
                let mut series = iter.next().unwrap().series_impl(x, expansion_point, info)?;
                for arg in iter {
                    series = &series + &arg.series_impl(x, expansion_point, info)?;
                }

                Ok(series)
            }
        }
    }
}

impl Mul<&Atom> for Series<AtomField> {
    type Output = Result<Series<AtomField>, SeriesError>;

    fn mul(self, rhs: &Atom) -> Result<Series<AtomField>, SeriesError> {
        (&self) * rhs
    }
}

impl Mul<&Series<AtomField>> for &Atom {
    type Output = Result<Series<AtomField>, SeriesError>;

    fn mul(self, rhs: &Series<AtomField>) -> Result<Series<AtomField>, SeriesError> {
        rhs * self
    }
}

impl Mul<&Series<AtomField>> for Atom {
    type Output = Result<Series<AtomField>, SeriesError>;

    fn mul(self, rhs: &Series<AtomField>) -> Result<Series<AtomField>, SeriesError> {
        rhs * &self
    }
}

impl Mul<&Atom> for &Series<AtomField> {
    type Output = Result<Series<AtomField>, SeriesError>;

    fn mul(self, rhs: &Atom) -> Result<Series<AtomField>, SeriesError> {
        let x = match self.get_variable().as_ref() {
            PolyVariable::Symbol(x) => Indeterminate::from(*x),
            PolyVariable::Function(_, x) => Indeterminate::try_from(x.clone()).map_err(|_| {
                SeriesError::NonIndeterminateSeriesVariable {
                    variable: self.get_variable().as_ref().clone(),
                }
            })?,
            PolyVariable::Power(_) | PolyVariable::Temporary(_) => {
                return Err(SeriesError::NonIndeterminateSeriesVariable {
                    variable: self.get_variable().as_ref().clone(),
                });
            }
        };

        let expansion_point = self.get_expansion_point();
        let mut current_depth = self.relative_order();

        if current_depth.is_zero() {
            current_depth = (2, 1).into();
        }

        loop {
            let info = Series::new(
                self.get_field(),
                None,
                self.get_variable().clone(),
                expansion_point.to_owned(),
                current_depth.clone(),
            );

            let series = rhs
                .as_view()
                .series_impl(&x, expansion_point.as_view(), &info)?
                * self;
            if series.relative_order() >= self.relative_order() {
                return Ok(series);
            } else {
                // increase the expansion depth
                current_depth = &current_depth * &2.into();
            }
        }
    }
}

impl Add<&Atom> for Series<AtomField> {
    type Output = Result<Series<AtomField>, SeriesError>;

    fn add(self, rhs: &Atom) -> Result<Series<AtomField>, SeriesError> {
        (&self) + rhs
    }
}

impl Add<&Series<AtomField>> for &Atom {
    type Output = Result<Series<AtomField>, SeriesError>;

    fn add(self, rhs: &Series<AtomField>) -> Result<Series<AtomField>, SeriesError> {
        rhs + self
    }
}

impl Add<&Series<AtomField>> for Atom {
    type Output = Result<Series<AtomField>, SeriesError>;

    fn add(self, rhs: &Series<AtomField>) -> Result<Series<AtomField>, SeriesError> {
        rhs + &self
    }
}

impl Add<&Atom> for &Series<AtomField> {
    type Output = Result<Series<AtomField>, SeriesError>;

    fn add(self, rhs: &Atom) -> Result<Series<AtomField>, SeriesError> {
        let x = match self.get_variable().as_ref() {
            PolyVariable::Symbol(x) => Indeterminate::from(*x),
            PolyVariable::Function(_, x) => Indeterminate::try_from(x.clone()).map_err(|_| {
                SeriesError::NonIndeterminateSeriesVariable {
                    variable: self.get_variable().as_ref().clone(),
                }
            })?,
            PolyVariable::Power(_) | PolyVariable::Temporary(_) => {
                return Err(SeriesError::NonIndeterminateSeriesVariable {
                    variable: self.get_variable().as_ref().clone(),
                });
            }
        };

        let expansion_point = self.get_expansion_point();
        let mut current_depth = self.relative_order();

        if current_depth.is_zero() {
            current_depth = (2, 1).into();
        }

        loop {
            let info = Series::new(
                self.get_field(),
                None,
                self.get_variable().clone(),
                expansion_point.to_owned(),
                current_depth.clone(),
            );

            let series = rhs
                .as_view()
                .series_impl(&x, expansion_point.as_view(), &info)?
                + self.clone();
            if series.absolute_order() >= self.absolute_order() {
                return Ok(series);
            } else {
                // increase the expansion depth
                current_depth = &current_depth * &2.into();
            }
        }
    }
}

impl Div<&Atom> for Series<AtomField> {
    type Output = Result<Series<AtomField>, SeriesError>;

    fn div(self, rhs: &Atom) -> Result<Series<AtomField>, SeriesError> {
        (&self) / rhs
    }
}

impl Div<&Series<AtomField>> for &Atom {
    type Output = Result<Series<AtomField>, SeriesError>;

    fn div(self, rhs: &Series<AtomField>) -> Result<Series<AtomField>, SeriesError> {
        rhs.rpow((-1, 1).into())? * self
    }
}

impl Div<&Series<AtomField>> for Atom {
    type Output = Result<Series<AtomField>, SeriesError>;

    fn div(self, rhs: &Series<AtomField>) -> Result<Series<AtomField>, SeriesError> {
        rhs.rpow((-1, 1).into())? * &self
    }
}

impl Div<&Atom> for &Series<AtomField> {
    type Output = Result<Series<AtomField>, SeriesError>;

    fn div(self, rhs: &Atom) -> Result<Series<AtomField>, SeriesError> {
        let x = match self.get_variable().as_ref() {
            PolyVariable::Symbol(x) => Indeterminate::from(*x),
            PolyVariable::Function(_, x) => Indeterminate::try_from(x.clone()).map_err(|_| {
                SeriesError::NonIndeterminateSeriesVariable {
                    variable: self.get_variable().as_ref().clone(),
                }
            })?,
            PolyVariable::Power(_) | PolyVariable::Temporary(_) => {
                return Err(SeriesError::NonIndeterminateSeriesVariable {
                    variable: self.get_variable().as_ref().clone(),
                });
            }
        };

        let expansion_point = self.get_expansion_point();
        let mut current_depth = self.relative_order();

        if current_depth.is_zero() {
            current_depth = (2, 1).into();
        }

        loop {
            let info = Series::new(
                self.get_field(),
                None,
                self.get_variable().clone(),
                expansion_point.to_owned(),
                current_depth.clone(),
            );

            let series = self
                / &rhs
                    .as_view()
                    .series_impl(&x, expansion_point.as_view(), &info)?;
            if series.relative_order() >= self.relative_order() {
                return Ok(series);
            } else {
                // increase the expansion depth
                current_depth = &current_depth * &2.into();
            }
        }
    }
}

impl Sub<&Atom> for Series<AtomField> {
    type Output = Result<Series<AtomField>, SeriesError>;

    fn sub(self, rhs: &Atom) -> Result<Series<AtomField>, SeriesError> {
        (&self) + &(-rhs)
    }
}

impl Sub<&Series<AtomField>> for &Atom {
    type Output = Result<Series<AtomField>, SeriesError>;

    fn sub(self, rhs: &Series<AtomField>) -> Result<Series<AtomField>, SeriesError> {
        -rhs.clone() + self
    }
}

impl Sub<&Series<AtomField>> for Atom {
    type Output = Result<Series<AtomField>, SeriesError>;

    fn sub(self, rhs: &Series<AtomField>) -> Result<Series<AtomField>, SeriesError> {
        -rhs.clone() + &self
    }
}

impl Sub<&Atom> for &Series<AtomField> {
    type Output = Result<Series<AtomField>, SeriesError>;

    fn sub(self, rhs: &Atom) -> Result<Series<AtomField>, SeriesError> {
        self + &(-rhs)
    }
}

#[cfg(test)]
mod test {
    use crate::{
        atom::{Atom, AtomCore, AtomView},
        parse,
        poly::series::SeriesDepth,
        symbol,
    };

    #[test]
    fn derivative() {
        let v1 = symbol!("v1");
        let inputs = [
            "(1+2*v1)^(5+v1)",
            "log(2*v1) + exp(3*v1) + sin(4*v1) + cos(y*v1)",
            "f(v1^2,v1)",
            "der(0,1,f,v1,v1^3)",
        ];
        let r = inputs.map(|input| parse!(input).derivative(v1));

        let res = [
            "(2*v1+1)^(v1+5)*log(2*v1+1)+2*(v1+5)*(2*v1+1)^(v1+4)",
            "2*(2*v1)^-1+3*exp(3*v1)+4*cos(4*v1)-y*sin(v1*y)",
            "der(0,1,f,v1^2,v1)+2*v1*der(1,0,f,v1^2,v1)",
            "der(1,1,f,v1,v1^3)+3*v1^2*der(0,2,f,v1,v1^3)",
        ];
        let res = res.map(|input| parse!(input));

        assert_eq!(r, res);
    }

    #[test]
    fn series() {
        let v1 = symbol!("v1");

        let input = parse!("exp(v1^2+1)*log(v1+3)/v1/(v1+1)");
        let t = input
            .series(v1, Atom::num(0).as_view(), SeriesDepth::absolute(2))
            .unwrap()
            .to_atom();

        let res = parse!(
            "1/3*exp(1)+v1*(-7/18*exp(1)+2*exp(1)*log(3))+v1^2*(119/162*exp(1)-2*exp(1)*log(3))-exp(1)*log(3)+v1^-1*exp(1)*log(3)"
        );
        assert_eq!(t, res);
    }

    #[test]
    fn series_shift() {
        let v1 = symbol!("v1");
        let input = parse!("1/(v1+1)");
        let t = input
            .series(v1, Atom::num(-1).as_view(), SeriesDepth::absolute(5))
            .unwrap()
            .to_atom();

        let res = parse!("1/(v1+1)");
        assert_eq!(t, res);
    }

    #[test]
    fn series_spurious_pole() {
        let v1 = symbol!("v1");
        let input = parse!("(1-cos(v1))/sin(v1)");
        let t = input
            .series(v1, Atom::num(0).as_view(), SeriesDepth::absolute(5))
            .unwrap()
            .to_atom();

        let res = parse!("1/2*v1+1/24*v1^3+1/240*v1^5");
        assert_eq!(t, res);
    }

    #[test]
    fn series_logx() {
        let v1 = symbol!("v1");
        let input = parse!("log(v1)*(1+v1)");
        let t = input
            .series(v1, Atom::num(0).as_view(), SeriesDepth::absolute(4))
            .unwrap()
            .to_atom();

        let res = parse!("log(v1)+v1*log(v1)");
        assert_eq!(t, res);
    }

    #[test]
    fn series_sqrt() {
        let v1 = symbol!("v1");
        let input = parse!("(v1^3+v1+1)^(1/2)");
        let t = input
            .series(v1, Atom::num(0).as_view(), SeriesDepth::absolute(4))
            .unwrap()
            .to_atom();

        let res = parse!("1+1/2*v1-1/8*v1^2+9/16*v1^3-37/128*v1^4");
        assert_eq!(t, res);
    }

    #[test]
    fn series_fractions() {
        let v1 = symbol!("v1");
        let input = parse!("1/v1^5");

        let t = input
            .series(v1, Atom::num(0).as_view(), SeriesDepth::absolute(3))
            .unwrap();

        let t2 = t.rpow((1, 3).into()).unwrap();

        assert_eq!(t2.absolute_order(), (22, 3));
    }

    #[test]
    fn series_zero() {
        let v1 = symbol!("v1");

        let input = parse!("1/v1^2+1/v1+v1");
        let t = input
            .series(v1, Atom::num(0).as_view(), SeriesDepth::absolute(0))
            .unwrap();

        assert_eq!(t.to_atom().expand(), parse!("v1^-2+v1^-1"));
    }

    #[test]
    fn series_poles() {
        let v1 = symbol!("v1");
        let input = parse!("1/(v1^10+v1^20)");

        let t = input
            .series(v1, Atom::num(0).as_view(), SeriesDepth::absolute(-1))
            .unwrap()
            .to_atom();

        assert_eq!(t, parse!("v1^-10"))
    }

    #[test]
    fn series_user_function() {
        let v1 = symbol!("v1");

        let input = parse!("f(exp(v1),sin(v1))");
        let t = input
            .series(v1, Atom::num(0).as_view(), SeriesDepth::absolute(2))
            .unwrap()
            .to_atom();

        let res = parse!(
            "f(1,0)+v1*(der(0,1,f,1,0)+der(1,0,f,1,0))+1/2*v1^2*(der(0,2,f,1,0)+der(1,0,f,1,0)+2*der(1,1,f,1,0)+der(2,0,f,1,0))"
        );
        assert_eq!(t, res);
    }

    #[test]
    fn series_derivative_keeps_function_unevaluated() {
        let v1 = symbol!("v1");
        let _ = symbol!(
            "series_derivative_keeps_function_unevaluated::f",
            norm = |input, out| {
                if let AtomView::Fun(f) = input
                    && f.get_nargs() == 1
                    && f.iter().next().unwrap() == Atom::num(2).as_view()
                {
                    out.to_num(1);
                }
            }
        );

        let input = parse!("series_derivative_keeps_function_unevaluated::f(v1)");
        let t = input
            .series(v1, Atom::num(2).as_view(), SeriesDepth::absolute(2))
            .unwrap()
            .to_atom();

        let res = parse!(
            "1+(-2+v1)*der(1,series_derivative_keeps_function_unevaluated::f,2)+1/2*(-2+v1)^2*der(2,series_derivative_keeps_function_unevaluated::f,2)"
        );
        assert_eq!(t, res);
    }

    #[test]
    fn series_custom_function() {
        let v1 = symbol!("v1");
        let _ = symbol!(
            "series_custom_function::f",
            series = |args| {
                let [arg] = args else { unreachable!() };
                let point = arg.coefficient(0.into());
                let delta = arg.to_atom() - &point;
                Some((Atom::num(1), Atom::num(1) / (Atom::num(1) - &delta)))
            }
        );

        let t = parse!("series_custom_function::f(v1)")
            .series(v1, Atom::num(0).as_view(), SeriesDepth::absolute(3))
            .unwrap()
            .to_atom();

        assert_eq!(t, parse!("1+v1+v1^2+v1^3"));
    }

    #[test]
    fn series_exp_log() {
        let v1 = symbol!("v1");

        let input = parse!("1+2*log(v1^4)");
        let t = input
            .series(v1, Atom::num(0).as_view(), SeriesDepth::absolute(4))
            .unwrap()
            .exp()
            .unwrap();

        assert_eq!(t.to_atom().expand(), parse!("v1^8*exp(1)"));
    }

    #[test]
    fn series_sub_atom() {
        let v1 = symbol!("v1");

        let input = parse!("1/(1-v1)");
        let t = input
            .series(v1, Atom::num(0).as_view(), SeriesDepth::absolute(4))
            .unwrap();

        let r = (t - &parse!("1/v1+1")).unwrap();

        assert_eq!(r.absolute_order(), (5, 1));
        assert_eq!(r.to_atom(), parse!("-1*v1^-1+v1+v1^2+v1^3+v1^4"));
    }

    #[test]
    fn series_div_atom() {
        let v1 = symbol!("v1");

        let input = parse!("v1");
        let t = input
            .series(v1, Atom::num(0).as_view(), SeriesDepth::absolute(4))
            .unwrap();

        let r = ((t / &parse!("exp(v1)-1")).unwrap() * &parse!("v1")).unwrap();

        assert_eq!(r.relative_order(), (4, 1));
        assert_eq!(r.to_atom(), parse!("v1+-1/2*v1^2+1/12*v1^3"));
    }

    #[test]
    fn series_relative_order() {
        let v1 = symbol!("v1");

        let input = parse!("exp(v1)/v1-1/6*v1^2");
        let t = input
            .series(v1, Atom::num(0).as_view(), SeriesDepth::relative(4))
            .unwrap();

        assert_eq!(t.relative_order(), (4, 1));
        assert_eq!(t.to_atom(), parse!("v1^-1+1+1/2*v1"));
    }

    #[test]
    fn series_truncate() {
        let v1 = symbol!("v1");

        let input = parse!("v1^10");
        let t = input
            .series(v1, Atom::num(0).as_view(), SeriesDepth::absolute(4))
            .unwrap();
        assert_eq!(t.absolute_order(), (10, 1));
        assert_eq!(t.relative_order(), (0, 1));
        assert_eq!(t.to_atom(), parse!("0"));

        let r = (&t * &input).unwrap();
        assert_eq!(r.absolute_order(), (20, 1));
    }

    #[test]
    fn series_empty() {
        let v1 = symbol!("v1");

        let input = parse!("v1");
        let t = input
            .series(v1, Atom::num(0).as_view(), SeriesDepth::absolute(4))
            .unwrap();

        let r = &t - &t;

        let t2 = parse!("v1^6")
            .series(v1, Atom::num(0).as_view(), SeriesDepth::relative(4))
            .unwrap();

        let x = (&r + &parse!("v1^6")).unwrap();
        assert_eq!(r.absolute_order(), (5, 1));

        let c = x.cos().unwrap();
        assert_eq!(c.absolute_order(), (10, 1));
        assert_eq!(c.relative_order(), (10, 1));

        let s = x.sin().unwrap();
        assert_eq!(s.absolute_order(), (5, 1));
        assert_eq!(s.relative_order(), (5, 1));

        let e = x.exp().unwrap();
        assert_eq!(e.absolute_order(), (5, 1));
        assert_eq!(e.relative_order(), (5, 1));

        let add = &r + &t2;
        assert_eq!(add.absolute_order(), (5, 1));
        let mul = &r * &t2;
        assert_eq!(mul.absolute_order(), (11, 1));
    }
}