computable_real/

real.rs

use super::{Approx, Primitive};
use num::bigint::Sign;
use num::{BigInt, One, Signed, ToPrimitive};
use std::cmp::Ordering;
use std::fmt::{self, Debug, Formatter};

/// 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 arbitrary precision.
#[derive(Clone)]
pub struct Real {
    pub prim: Box<Primitive>,
    pub approx: Option<Approx>,
}

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

const LN_MIN: i32 = 8; // 0.5
const LN_MAX: i32 = 24; // 1.5

impl Real {
    pub fn new(prim: Primitive) -> Self {
        Real {
            prim: Box::new(prim),
            approx: None,
        }
    }

    /// Returns the most significant digit based on the current approximation.
    /// If msd is n, then 2^(n-1) < |x| < 2^(n+1).
    /// Returns None if no approximation has been made.
    pub 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)
    }

    // TODO: understand the comparision with 1
    /// Evaluates to the desired precision and returns the msd.
    /// Returns None if the correct msd < precision
    pub 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 precision is reached.
    pub 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 create an infinite loop if the value is zero.
    pub 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.
    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
    }

    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)
    }

    pub fn render(&self, n: u32) -> String {
        self.render_base(n, 10)
    }

    pub fn cmp_ra(&mut self, other: &mut Real, rel_tol: i32, abs_tol: i32) -> Ordering {
        let msd1 = self.iter_msd_n(abs_tol).unwrap_or(i32::MIN);
        let msd2 = other
            .iter_msd_n(if msd1 > abs_tol { msd1 } else { abs_tol })
            .unwrap_or(i32::MIN);
        let msd = msd1.max(msd2);
        if msd == i32::MIN {
            return Ordering::Equal;
        }
        let rel = msd + rel_tol;
        self.cmp_a(other, rel.max(abs_tol))
    }

    pub fn cmp_a(&mut self, other: &mut Real, abs_tol: i32) -> Ordering {
        let prec = abs_tol - 1;
        let a = self.appr(prec).value;
        let b = other.appr(prec).value;
        a.cmp(&b)
    }

    pub fn cmp_until(&self, other: &Real, min_precision: i32) -> Ordering {
        let mut a = -16;
        let mut clone = self.clone();
        let mut other = other.clone();
        loop {
            let cmp = clone.cmp_a(&mut other, a);
            if cmp != Ordering::Equal || a < min_precision {
                return cmp;
            }
            a *= 2;
        }
    }

    pub fn signum_a(&self, a: i32) -> i32 {
        if let Some(res) = &self.approx {
            if res.value.signum() != BigInt::ZERO {
                return res.value.signum().to_i32().unwrap();
            }
        }
        let prec = a - 1;
        self.clone().appr(prec).value.signum().to_i32().unwrap()
    }
}

impl Real {
    pub fn int(value: BigInt) -> Self {
        Self::new(Primitive::Int(value))
    }

    pub fn squared(self) -> Self {
        Self::new(Primitive::Square(self))
    }

    pub fn atan(value: BigInt) -> Self {
        Self::new(Primitive::Atan(value))
    }

    pub fn exp(mut self) -> Self {
        let rough_appr = self.appr(-10).value;
        let two = BigInt::from(2);
        if rough_appr.abs() > two {
            self.shifted(-1).exp().squared()
        } else {
            Self::new(Primitive::Exp(self))
        }
    }

    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 {
            let pi_multiples: BigInt = (half_pi_multiples + 1) / 2;
            let adj = Self::pi() * Real::int(pi_multiples.clone());
            if (pi_multiples & BigInt::one()).signum() != BigInt::ZERO {
                -(self - adj).cos()
            } else {
                (self - adj).cos()
            }
        } else if self.appr(-1).value.abs() >= two {
            self.shifted(-1).cos().squared().shifted(1) - Self::one()
        } else {
            Self::new(Primitive::Cos(self))
        }
    }

    pub fn sin(self) -> Self {
        (Self::pi().shifted(-1) - self).cos()
    }

    pub fn asin(mut self) -> Self {
        let rough_appr = self.appr(-10).value;
        if rough_appr > 750.into() {
            Self::one() - self.squared().sqrt().acos()
        } else if rough_appr < (-750).into() {
            -(-self.asin())
        } else {
            Self::new(Primitive::Asin(self))
        }
    }

    pub fn acos(self) -> Self {
        Self::pi().shifted(-1) - self.asin()
    }

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

    pub fn ln(mut self) -> Self {
        let rough_appr = self.appr(-4).value;
        if rough_appr <= LN_MIN.into() {
            -((Self::one() / self).ln())
        } else if rough_appr >= LN_MAX.into() {
            if rough_appr <= 64.into() {
                self.sqrt().sqrt().ln().shifted(2)
            } else {
                let ln2_1 =
                    Real::int(7.into()) * (Real::int(10.into()) / Real::int(9.into())).ln2();
                let ln2_2 =
                    Real::int(2.into()) * (Real::int(25.into()) / Real::int(24.into())).ln2();
                let ln2_3 =
                    Real::int(3.into()) * (Real::int(81.into()) / Real::int(80.into())).ln2();
                let ln2 = ln2_1 - ln2_2 + ln2_3;
                let extra_bits = rough_appr.bits() as i32 - 3;
                self.shifted(-extra_bits).ln() + Real::int(extra_bits.into()) * ln2
            }
        } else {
            self.ln2()
        }
    }

    pub fn sqrt(self) -> Self {
        Self::new(Primitive::Sqrt(self))
    }

    pub fn shifted(self, bits: i32) -> Self {
        Self::new(Primitive::Shift(self, bits))
    }

    /// John Machin's formula from 1706:
    /// pi/4 = 4 * atan(1/5) - atan(1/239)
    pub fn pi() -> Self {
        let atan5 = Real::atan(BigInt::from(5));
        let atan239 = Real::atan(BigInt::from(239));
        let four = Real::int(BigInt::from(4));
        let four_atan5 = four.clone() * atan5;
        let diff = four_atan5 - atan239;
        diff * four
    }

    pub fn inv(self) -> Self {
        Self::new(Primitive::Inv(self))
    }
}