computable_real/

real.rs

use crate::{Approx, Primitive};

use std::cmp::Ordering;
use std::fmt::{self, Debug, Formatter};

use num::bigint::Sign;
use num::{BigInt, Integer, One, Signed, ToPrimitive, Zero};

/// The main structure representing a real number.
///
/// It contains a private field which determines the function and the arguments,
/// and the current best approximation, if any.
///
/// The core idea here is that we represent a real number as a function
/// `f(n: i32) -> BigInt`, where `n` is the desired precision (typically
/// negative). The function is constructed such that `|f(n)*2^n - x| < 2^n`.
/// By using an increasingly negative `n`, we can approximate `x` to an
/// arbitrary precision.
///
/// # Warning on Undefined Behavior
///
/// <div class="warning">
/// Comparing two computable real numbers is inefficient, and indecidable in
/// general. Because of this, we do not validate whether the arguments to the
/// functions are valid. It is the user's responsibility to ensure that the
/// arguments are in the function domain. Anything else will cause undefined
/// behavior.
/// </div>
///
/// The following things will cause undefined behavior, which may return an
/// incorrect result, create an infinite loop, or panic:
///
/// - Division by zero
/// - Calling [`Real::sqrt`] on a negative number
/// - Calling [`Real::ln`] on a non-positive number
/// - Calling [`Real::tan`] on a number where `cos` is zero
/// - Calling [`Real::asin`] or [`Real::acos`] on a number outside [-1, 1]
#[derive(Clone)]
pub struct Real {
    /// The primitive that represents the real number.
    prim: Box<Primitive>,
    /// The current best approximation of the real number.
    pub approx: Option<Approx>,
}

impl Debug for Real {
    fn fmt(&self, f: &mut Formatter) -> fmt::Result {
        write!(f, "{:?}", self.prim)
    }
}

impl Real {
    /// Creates a new real number from a primitive.
    ///
    /// Normally you would not use this directly, but rather use one of the
    /// other associated functions like [`Real::int`] or [`Real::pi`].
    pub fn new(prim: Primitive) -> Self {
        Real {
            prim: Box::new(prim),
            approx: None,
        }
    }

    /// Returns the most significant digit based on the current approximation,
    /// or None if no approximation has been made.
    ///
    /// If msd is `n`, then `2^(n-1) < |x| < 2^(n+1)`.
    pub(crate) fn known_msd(&self) -> Option<i32> {
        let res = self.approx.as_ref()?;
        let size = res.value.abs().bits() as i32;
        Some(res.precision + size - 1)
    }

    /// Evaluates to the desired precision and returns the msd, or None if the
    /// correct msd < precision.
    pub(crate) fn msd(&mut self, precision: i32) -> Option<i32> {
        if let Some(res) = &self.approx {
            if res.value.abs() > BigInt::one() {
                return self.known_msd();
            }
        }
        let res = self.appr(precision - 1);
        if res.value.abs() <= BigInt::one() {
            // msd could still be arbitrarily far to the right
            None
        } else {
            self.known_msd()
        }
    }

    /// Iteratively lowers the precision until the msd is found, or the
    /// provided minimum precision is reached.
    pub(crate) fn iter_msd_n(&mut self, min_precision: i32) -> Option<i32> {
        let mut curr = 0;
        loop {
            let msd = self.msd(curr);
            if msd.is_some() {
                return msd;
            }
            curr = (curr * 3) / 2 - 16;
            if curr <= min_precision + 30 {
                break;
            }
        }
        self.msd(min_precision)
    }

    /// Iteratively lowers the precision until the msd is found. This will
    /// cause an infinite loop if the value is zero.
    pub(crate) fn iter_msd(&mut self) -> Option<i32> {
        self.iter_msd_n(i32::MIN)
    }

    /// Generates an approximation of the real number.
    ///
    /// If we have already approximated the number with a sufficient precision,
    /// we simply scale the approximation. Otherwise we call the underlying
    /// approximate method.
    ///
    /// # Examples
    /// ```
    /// use computable_real::Real;
    ///
    /// let mut n = Real::from(42);
    /// assert_eq!(n.appr(0).value, 42.into()); // 42*2^0 = 42
    /// assert_eq!(n.appr(-1).value, 84.into()); // 84*2^-1 = 42
    /// ```
    pub fn appr(&mut self, precision: i32) -> Approx {
        if let Some(res) = &self.approx {
            if precision >= res.precision {
                return res.scale_to(precision);
            }
        }
        let res = self.prim.approximate(precision, self.clone());
        self.approx = Some(res.clone());
        res
    }

    /// Renders the number in `base` with `n` digits to the right of
    /// decimal point.
    ///
    /// # Examples
    /// ```
    /// use computable_real::Real;
    ///
    /// let n = Real::from(42);
    /// assert_eq!(n.render_base(0, 10), "42");
    /// assert_eq!(n.render_base(2, 10), "42.00");
    /// assert_eq!(n.render_base(4, 2), "101010.0000");
    /// assert_eq!(n.render_base(0, 16), "2a");
    ///
    /// let pi = Real::pi();
    /// assert_eq!(pi.render_base(10, 10), "3.1415926536");
    /// assert_eq!(pi.render_base(10, 16), "3.243f6a8886");
    /// assert_eq!(pi.render_base(10, 2), "11.0010010001");
    /// ```
    pub fn render_base(&self, n: u32, base: u32) -> String {
        let scaled = match base {
            16 => self.clone().shifted(n as i32 * 4),
            _ => self.clone() * Real::int(BigInt::from(base).pow(n)),
        }
        .appr(0);

        let value_str = scaled.value.abs().to_str_radix(base);

        let padded_str = match value_str.len() as u32 {
            len if len <= n => format!("{}{}", "0".repeat((n - len + 1) as usize), value_str),
            _ => value_str,
        };

        let len = padded_str.len() as u32;
        let (whole, frac) = padded_str.split_at((len - n) as usize);

        let sign = if scaled.value.sign() == Sign::Minus {
            "-"
        } else {
            ""
        };
        let dot = if n == 0 { "" } else { "." };
        format!("{}{}{}{}", sign, whole, dot, frac)
    }

    /// Renders the number with `n` digits to the right of the decimal point.
    ///
    /// Equivalent to `render_base(n, 10)`.
    ///
    /// # Examples
    /// ```
    /// use computable_real::Real;
    ///
    /// let n = Real::from(42);
    /// assert_eq!(n.render(0), "42");
    /// assert_eq!(n.render(2), "42.00");
    ///
    /// let zero = Real::from(0);
    /// assert_eq!(zero.render(0), "0");
    /// assert_eq!(zero.render(2), "0.00");
    ///
    /// let minus_one = Real::from(-1);
    /// assert_eq!(minus_one.render(0), "-1");
    /// assert_eq!(minus_one.render(2), "-1.00");
    /// ```
    pub fn render(&self, n: u32) -> String {
        self.render_base(n, 10)
    }
}

/// Convenient methods and associated functions to create computable reals.
impl Real {
    /// Creates a computable real representing an integer.
    ///
    /// # Examples
    /// ```
    /// use computable_real::Real;
    ///
    /// let n = Real::int(42.into());
    /// assert_eq!(n.render(0), "42");
    /// ```
    pub fn int(n: BigInt) -> Self {
        Self::new(Primitive::Int(n))
    }

    /// Creates a computable real representing zero.
    ///
    /// Equivalent to `Real::int(BigInt::ZERO)`.
    pub fn zero() -> Self {
        Self::int(BigInt::ZERO)
    }

    /// Creates a computable real representing one.
    ///
    /// Equivalent to `Real::int(BigInt::one())`.
    pub fn one() -> Self {
        Self::int(BigInt::one())
    }

    /// Creates a computable real representing pi.
    ///
    /// This uses John Machin's formula from 1706:
    /// `pi/4 = 4 * atan(1/5) - atan(1/239)`
    ///
    /// # Examples
    /// ```
    /// use computable_real::Real;
    ///
    /// let pi = Real::pi();
    /// assert_eq!(pi.render(10), "3.1415926536");
    /// ```
    pub fn pi() -> Self {
        let atan5 = Self::new(Primitive::Atan(BigInt::from(5)));
        let atan239 = Self::new(Primitive::Atan(BigInt::from(239)));
        let four_atan5 = Real::from(4) * atan5;
        let diff = four_atan5 - atan239;
        diff * Real::from(4)
    }

    /// Creates a computable real representing `self^2`.
    ///
    /// # Examples
    /// ```
    /// use computable_real::Real;
    ///
    /// assert_eq!(Real::from(42).squared().render(0), "1764");
    /// assert_eq!(Real::from(-42).squared().render(0), "1764");
    /// assert_eq!(Real::from(1).squared().render(0), "1");
    /// ```
    pub fn squared(self) -> Self {
        Self::new(Primitive::Square(self))
    }

    /// Creates a computable real representing `a` if `self < 0`, and `b`
    /// otherwise.
    ///
    /// # Examples
    /// ```
    /// use computable_real::Real;
    ///
    /// let n = Real::from(42);
    /// let m = Real::from(-42);
    /// assert_eq!(n.clone().select(n.clone(), m.clone()).render(0), "-42");
    /// assert_eq!(m.clone().select(n.clone(), m.clone()).render(0), "42");
    /// ```
    pub fn select(self, a: Self, b: Self) -> Self {
        Self::new(Primitive::Select(self, a, b))
    }

    /// Creates a computable real representing `exp(self)`.
    ///
    /// # Examples
    /// ```
    /// use computable_real::Real;
    ///
    /// let n = Real::one();
    /// let e = n.exp();
    /// assert_eq!(e.render(10), "2.7182818285");
    /// ```
    pub fn exp(mut self) -> Self {
        let rough_appr = self.appr(-10).value;
        let two = BigInt::from(2);
        if rough_appr.abs() > two {
            // If x >= 2*2^-10, we scale the argument into range
            // exp(x) = exp(x/2)^2
            self.shifted(-1).exp().squared()
        } else {
            Self::new(Primitive::Exp(self))
        }
    }

    /// Creates a computable real representing `cos(self)`.
    ///
    /// # Examples
    /// ```
    /// use computable_real::Real;
    ///
    /// let cos_pi = Real::pi().cos();
    /// assert_eq!(cos_pi.render(0), "-1");
    ///
    /// let cos_pi_half = Real::pi().shifted(-1).cos();
    /// assert_eq!(cos_pi_half.render(0), "0");
    ///
    /// let cos_n3_pi = (Real::pi() * Real::from(-3)).cos();
    /// assert_eq!(cos_n3_pi.render(0), "-1");
    ///
    /// let cos_pi_over_3 = (Real::pi() / Real::from(3)).cos();
    /// assert_eq!(cos_pi_over_3.render(2), "0.50");
    /// ```
    pub fn cos(mut self) -> Self {
        let half_pi_multiples = (self.clone() / Self::pi()).appr(-1).value;
        let two = BigInt::from(2);

        if half_pi_multiples.abs() >= two {
            // If |arg| >= pi, we move the argument into range
            let pi_multiples: BigInt = super::scale(half_pi_multiples.clone(), -1);
            let adj = Self::pi() * Real::from(pi_multiples.clone());
            if pi_multiples.is_even() {
                (self - adj).cos()
            } else {
                -(self - adj).cos()
            }
        } else if self.appr(-1).value.abs() >= two {
            // If 1 <= |arg| <= pi, use cos double angle formula:
            // cos(x) = 2 * cos(x/2)^2  - 1
            self.shifted(-1).cos().squared().shifted(1) - Self::one()
        } else {
            // If |arg| < 1
            Self::new(Primitive::Cos(self))
        }
    }

    /// Creates a computable real representing `sin(self)`.
    ///
    /// # Examples
    /// ```
    /// use computable_real::{Real, Comparator};
    /// use std::cmp::Ordering;
    ///
    /// let sin_pi = Real::pi().sin();
    /// assert_eq!(sin_pi.render(0), "0");
    ///
    /// let sin_pi_half = Real::pi().shifted(-1).sin();
    /// assert_eq!(sin_pi_half.render(0), "1");
    ///
    /// let sin_n3_pi = (Real::pi() * Real::from(-3)).sin();
    /// assert_eq!(sin_n3_pi.render(0), "0");
    ///
    /// // sin(pi/3) = sqrt(3)/2
    /// let cmp = Comparator::new().min_prec(-10);
    /// let mut sin_pi_over_3 = (Real::pi() / Real::from(3)).sin();
    /// let mut sqrt_3_over_2 = Real::from(3).sqrt() / Real::from(2);
    /// assert_eq!(cmp.cmp(&mut sin_pi_over_3, &mut sqrt_3_over_2), Ordering::Equal);
    /// ```
    pub fn sin(self) -> Self {
        // sin(x) = cos(pi/2 - x)
        (Self::pi().shifted(-1) - self).cos()
    }

    /// Creates a computable real representing `tan(self)`.
    ///
    /// # Examples
    /// ```
    /// use computable_real::Real;
    ///
    /// let tan_pi = Real::pi().tan();
    /// assert_eq!(tan_pi.render(0), "0");
    ///
    /// let tan_n3_pi = (Real::pi() * Real::from(-3)).tan();
    /// assert_eq!(tan_n3_pi.render(0), "0");
    ///
    /// let tan_pi_over_3 = (Real::pi() / Real::from(3)).tan();
    /// assert_eq!(tan_pi_over_3.render(0), "2");
    /// ```
    pub fn tan(self) -> Self {
        // tan(x) = sin(x) / cos(x)
        self.clone().sin() / self.cos()
    }

    /// Creates a computable real representing `asin(self)`.
    ///
    /// # Examples
    /// ```
    /// use computable_real::{Real, Comparator};
    /// use std::cmp::Ordering;
    ///
    /// let asin_zero = Real::from(0).asin();
    /// assert_eq!(asin_zero.render(0), "0");
    ///
    /// let cmp = Comparator::new().min_prec(-10);
    ///
    /// let mut asin_one = Real::from(1).asin();
    /// let mut half_pi = Real::pi().shifted(-1);
    /// assert_eq!(cmp.cmp(&mut asin_one, &mut half_pi), Ordering::Equal);
    ///
    /// let n_one_sqrt_2 = Real::from(-1) / Real::from(2).sqrt();
    /// let mut val = n_one_sqrt_2.asin();
    /// let mut n_pi_over_4 = -Real::pi().shifted(-2);
    /// assert_eq!(cmp.cmp(&mut val, &mut n_pi_over_4), Ordering::Equal);
    /// ```
    pub fn asin(mut self) -> Self {
        let rough_appr = self.appr(-10).value;
        // 750*2^-10 is slightly larger than sin(pi/4) = 1/sqrt(2)
        if rough_appr > 750.into() {
            // Use asin(x) = acos(sqrt(1 - x^2)), follows from sin^2 + cos^2 = 1
            (Self::one() - self.squared()).sqrt().acos()
        } else if rough_appr < (-750).into() {
            -(-self.asin())
        } else {
            Self::new(Primitive::Asin(self))
        }
    }

    /// Creates a computable real representing `acos(self)`.
    ///
    /// # Examples
    /// ```
    /// use computable_real::{Real, Comparator};
    /// use std::cmp::Ordering;
    ///
    /// let cmp = Comparator::new().min_prec(-10);
    ///
    /// let mut acos_zero = Real::from(0).acos();
    /// let mut half_pi = Real::pi().shifted(-1);
    /// assert_eq!(cmp.cmp(&mut acos_zero, &mut half_pi), Ordering::Equal);
    ///
    /// let mut acos_one = Real::from(1).acos();
    /// assert_eq!(cmp.cmp(&mut acos_one, &mut Real::from(0)), Ordering::Equal);
    ///
    /// let mut acos_half = Real::from(1).shifted(-1).acos();
    /// let mut pi_over_3 = Real::pi() / Real::from(3);
    /// assert_eq!(cmp.cmp(&mut acos_half, &mut pi_over_3), Ordering::Equal);
    /// ```
    pub fn acos(self) -> Self {
        // acos(x) = pi/2 - asin(x)
        Self::pi().shifted(-1) - self.asin()
    }

    /// Creates a computable real representing `atan(self)`.
    ///
    /// # Examples
    /// ```
    /// use computable_real::{Real, Comparator};
    /// use std::cmp::Ordering;
    ///
    /// let atan_zero = Real::from(0).atan();
    /// assert_eq!(atan_zero.render(0), "0");
    ///
    /// let cmp = Comparator::new().min_prec(-10);
    ///
    /// let mut atan_one = Real::from(1).atan();
    /// let mut pi_over_4 = Real::pi().shifted(-2);
    /// assert_eq!(cmp.cmp(&mut atan_one, &mut pi_over_4), Ordering::Equal);
    ///
    /// let n_one_sqrt_3 = Real::from(-1) / Real::from(3).sqrt();
    /// let mut val = n_one_sqrt_3.atan();
    /// let mut n_pi_over_6 = -Real::pi() / Real::from(6);
    /// assert_eq!(cmp.cmp(&mut val, &mut n_pi_over_6), Ordering::Equal);
    /// ```
    pub fn atan(self) -> Self {
        // sin(x)^2 = (tan(x)^2) / (1 + tan(x)^2)
        let square = self.clone().squared();
        let abs_sin_atan = (square.clone() / (Self::one() + square)).sqrt();
        let sin_atan = self.select(-abs_sin_atan.clone(), abs_sin_atan);
        sin_atan.asin()
    }

    fn simple_ln(self) -> Self {
        Self::new(Primitive::Ln(self - Self::one()))
    }

    /// Creates a computable real representing `ln(self)`, the natural logarithm.
    ///
    /// # Examples
    /// ```
    /// use computable_real::Real;
    ///
    /// let ln1 = Real::from(1).ln();
    /// assert_eq!(ln1.render(0), "0");
    ///
    /// let e = Real::from(1).exp();
    /// let ln_e = e.ln();
    /// assert_eq!(ln_e.render(0), "1");
    /// ```
    pub fn ln(mut self) -> Self {
        let rough_appr = self.appr(-4).value;
        if rough_appr <= 8.into() {
            // If x <= 8*2^-4 = 0.5, we use ln(x) = -ln(1/x)
            -((Self::one() / self).ln())
        } else if rough_appr >= 24.into() {
            // If x >= 24*2^-4 = 1.5:
            if rough_appr <= 64.into() {
                // If 1.5 <= x <= 64*2^-4 = 4, we use ln(x) = 4*ln(sqrt(sqrt(x)))
                self.sqrt().sqrt().ln().shifted(2)
            } else {
                // If x > 4:
                let ln2_1 = Real::from(7) * (Real::from(10) / Real::from(9)).simple_ln();
                let ln2_2 = Real::from(2) * (Real::from(25) / Real::from(24)).simple_ln();
                let ln2_3 = Real::from(3) * (Real::from(81) / Real::from(80)).simple_ln();
                // 7*ln(10/9) - 2*ln(25/24) + 3*ln(81/80) = ln(2)
                let ln2 = ln2_1 - ln2_2 + ln2_3;
                let extra_bits = rough_appr.bits() as i32 - 3;
                // ln(x) = ln(x/2^extra_bits) + extra_bits*ln(2)
                self.shifted(-extra_bits).ln() + Real::from(extra_bits) * ln2
            }
        } else {
            // If 0.5 < x < 1.5, we can directly use the primitive
            self.simple_ln()
        }
    }

    /// Creates a computable real representing `sqrt(self)`.
    ///
    /// # Examples
    /// ```
    /// use computable_real::Real;
    ///
    /// let sqrt_2 = Real::from(2).sqrt();
    /// assert_eq!(sqrt_2.render(10), "1.4142135624");
    ///
    /// let sqrt_1 = Real::from(1).sqrt();
    /// assert_eq!(sqrt_1.render(0), "1");
    ///
    /// let sqrt_0 = Real::from(0).sqrt();
    /// assert_eq!(sqrt_0.render(0), "0");
    ///
    /// let sqrt_quarter = Real::from(1).shifted(-2).sqrt();
    /// assert_eq!(sqrt_quarter.render(3), "0.500");
    /// ```
    pub fn sqrt(self) -> Self {
        Self::new(Primitive::Sqrt(self))
    }

    /// Creates a computable real representing `self * 2^n`.
    ///
    /// # Examples
    /// ```
    /// use computable_real::Real;
    ///
    /// assert_eq!(Real::from(42).shifted(0).render(0), "42");
    /// assert_eq!(Real::from(42).shifted(2).render(0), "168");
    /// assert_eq!(Real::from(42).shifted(-2).render(3), "10.500");
    /// assert_eq!(Real::from(42).shifted(-2).shifted(2).render(0), "42");
    /// ```
    pub fn shifted(self, n: i32) -> Self {
        Self::new(Primitive::Shift(self, n))
    }

    /// Creates a computable real representing `1/self`.
    ///
    /// # Examples
    /// ```
    /// use computable_real::Real;
    ///
    /// assert_eq!(Real::from(42).inv().render(10), "0.0238095238");
    /// assert_eq!(Real::from(1).inv().render(0), "1");
    /// assert_eq!(Real::from(-1).inv().render(0), "-1");
    /// assert_eq!(Real::from(42).inv().inv().render(0), "42");
    /// ```
    pub fn inv(self) -> Self {
        Self::new(Primitive::Inv(self))
    }
}

/// Utility struct to compare real numbers.
///
/// One limitiation of computable real numbers is that, in general whether two
/// numbers are equal is undecidable. We could evaluate them to incraesingly
/// high precision until they diverge, but if they are in fact equal, this
/// would never terminate.
///
/// Given this, we provide the [Comparator] struct to compare real numbers and
/// get the sign of a real number. You will have to provide some kind of
/// tolerance level for the comparison.
///
/// Three different precision levels can be configured, and they are all
/// typically negative:
/// - [`Comparator::rel`] sets the relative precision level
/// - [`Comparator::abs`] sets the absolute precision level
/// - [`Comparator::min_prec`] sets the minimum precision level
///
/// The relative level is added on top of `max(msd(a), msd(b))` to determine the
/// precision we evaluate the numbers to to perform the comparison. Here `msd`
/// is defined as the binary position of the most significant digit, such that
/// `2^(msd-1) < |x| < 2^(msd+1)`. Note that the relative level is only used
/// if the absolute level is set.
///
/// If no relative level is set, the numbers will be evaluated to the absolute
/// level for comparison. If the relative level is set, the precision will be
/// taken as the maximum between the relative and absolute levels.
///
/// The minimum precision level sets a bound of the evaluation precision. We
/// will never evaluate the numbers to a precision beyond this. When no
/// absolute level is set, we will iteratively evaluate the numbers until they
/// diverge, or until the minimum precision level is reached.
///
/// # Examples
///
/// ```
/// use computable_real::{Comparator, Real};
/// use std::cmp::Ordering;
///
/// let mut one = Real::from(1);
/// let mut pi = Real::pi();
///
/// let cmp = Comparator::new();
/// assert_eq!(cmp.cmp(&mut one, &mut pi), Ordering::Less);
///
/// // A comparator that will only evaluate to 2^-10 precision
/// let lazy_cmp = Comparator::new().abs(-10);
/// let mut smol = Real::from(1).shifted(-100);
/// // It cannot tell smol from 0
/// assert_eq!(lazy_cmp.signum(&mut smol), 0);
/// let mut pi_smol = pi.clone() + smol.clone();
/// // It cannot tell pi + smol from pi
/// assert_eq!(lazy_cmp.cmp(&mut pi_smol, &mut pi), Ordering::Equal);
///
/// // The default comparator can correctly compare
/// assert_eq!(cmp.signum(&mut smol), 1);
/// assert_eq!(cmp.cmp(&mut pi_smol, &mut pi), Ordering::Greater);
/// ```
pub struct Comparator {
    rel: Option<i32>,
    abs: Option<i32>,
    min: i32,
}

impl Default for Comparator {
    /// Equivalent to [`Comparator::new`].
    fn default() -> Self {
        Self::new()
    }
}

impl Comparator {
    /// Returns a new comparator with no relative and absolute precision levels,
    /// and a minimum precision level of `i32::MIN`.
    ///
    /// <div class="warning">
    /// You can use the returned comparator as is, but if you compare two equal
    /// numbers or call [`Comparator::signum`] on zero, it can cause an
    /// infinite loop.
    ///
    /// When two numbers may be equal, you should set at least the absolute
    /// precision or the minimum precision level.
    /// </div>
    pub fn new() -> Self {
        Self {
            rel: None,
            abs: None,
            min: i32::MIN,
        }
    }

    /// Sets the relative precision level.
    pub fn rel(mut self, value: i32) -> Self {
        self.rel = Some(value);
        self
    }

    /// Sets the absolute precision level.
    pub fn abs(mut self, value: i32) -> Self {
        self.abs = Some(value);
        self
    }

    /// Sets the minimum precision level.
    pub fn min_prec(mut self, value: i32) -> Self {
        self.min = value;
        self
    }

    /// Compares two real numbers.
    ///
    /// See examples in [`Comparator#examples`].
    pub fn cmp(&self, a: &mut Real, b: &mut Real) -> Ordering {
        if let Some(abs) = self.abs {
            let tol = if let Some(rel) = self.rel {
                let msd = match (a.iter_msd_n(abs), b.iter_msd_n(abs)) {
                    (Some(msd1), Some(msd2)) => msd1.max(msd2),
                    (Some(msd1), None) => msd1,
                    (None, Some(msd2)) => msd2,
                    (None, None) => return Ordering::Equal,
                };
                abs.max(msd + rel)
            } else {
                abs
            };

            let prec = self.min.max(tol - 1);
            let a = a.appr(prec).value;
            let b = b.appr(prec).value;
            return a.cmp(&b);
        }

        let mut prec = -16;
        let mut cmp = Comparator::new().min_prec(self.min);
        loop {
            cmp.abs = Some(prec);
            let res = cmp.cmp(a, b);
            prec *= 2;
            if res != Ordering::Equal || prec <= self.min {
                return res;
            }
        }
    }

    /// Returns the sign of a real number.
    ///
    /// See examples in [`Comparator#examples`].
    pub fn signum(&self, n: &mut Real) -> i32 {
        if let Some(res) = &n.approx {
            if !res.value.is_zero() {
                return res.value.signum().to_i32().unwrap();
            }
        }

        match self.cmp(n, &mut Real::from(0)) {
            Ordering::Less => -1,
            Ordering::Equal => 0,
            Ordering::Greater => 1,
        }
    }
}