feat(math): add support for some math functions (#483)

pallas-math/Cargo.toml

@@ -10,9 +10,19 @@ license = "Apache-2.0"
 readme = "README.md"
 authors = ["Andrew Westberg <andrewwestberg@gmail.com>"]
 
+[features]
+default = ["gmp"]
+gmp = ["dep:gmp-mpfr-sys"]
+num = ["dep:num-bigint", "dep:num-integer", "dep:num-traits"]
 
 [dependencies]
-# rug = "1.24.1"
+gmp-mpfr-sys = { version = "1.6.4", features = ["mpc"], default-features = false, optional = true }
+once_cell = "1.19.0"
+num-bigint = { version = "0.4.6", optional = true }
+num-integer = { version = "0.1.46", optional = true }
+num-traits = { version = "0.2.19", optional = true }
+regex = "1.10.5"
+thiserror = "1.0.61"
 
 [dev-dependencies]
 quickcheck = "1.0"

pallas-math/src/lib.rs

@@ -1 +1,14 @@
 pub mod math;
+
+// Ensure only one of `gmp` or `num` is enabled, not both.
+#[cfg(all(feature = "gmp", feature = "num"))]
+compile_error!("Features `gmp` and `num` are mutually exclusive.");
+
+#[cfg(all(not(feature = "gmp"), not(feature = "num")))]
+compile_error!("One of the features `gmp` or `num` must be enabled.");
+
+#[cfg(feature = "gmp")]
+pub mod math_gmp;
+
+#[cfg(feature = "num")]
+pub mod math_num;

pallas-math/src/math.rs

@@ -1,17 +1,256 @@
 /*!
 # Cardano Math functions
-*/
+ */
 
-pub fn add(a: i32, b: i32) -> i32 {
-    a + b
+use std::fmt::{Debug, Display};
+use std::ops::{Div, Mul, Neg, Sub};
+
+use thiserror::Error;
+
+#[cfg(feature = "gmp")]
+use crate::math_gmp::Decimal;
+#[cfg(feature = "num")]
+use crate::math_num::Decimal;
+
+#[derive(Debug, Error)]
+pub enum Error {
+    #[error("error in regex")]
+    RegexFailure(#[from] regex::Error),
+
+    #[error("string contained a nul byte")]
+    NulFailure(#[from] std::ffi::NulError),
+}
+
+pub const DEFAULT_PRECISION: u64 = 34;
+
+pub trait FixedPrecision:
+    Neg + Mul + Div + Sub + Display + Clone + PartialEq + PartialOrd + Debug + From<u64> + From<i64>
+{
+    /// Creates a new fixed point number with the given precision
+    fn new(precision: u64) -> Self;
+
+    /// Creates a new fixed point number from an integer string. Precision tells us how many decimals
+    fn from_str(s: &str, precision: u64) -> Result<Self, Error>;
+
+    /// Returns the precision of the fixed point number
+    fn precision(&self) -> u64;
+
+    /// Performs the 'exp' approximation. First does the scaling of 'x' to [0,1]
+    /// and then calls the continued fraction approximation function.
+    fn exp(&self) -> Self;
+
+    /// Entry point for 'ln' approximation. First does the necessary scaling, and
+    /// then calls the continued fraction calculation. For any value outside the
+    /// domain, i.e., 'x in (-inf,0]', the function returns '-INFINITY'.
+    fn ln(&self) -> Self;
+
+    /// Entry point for 'pow' function. x^y = exp(y * ln x)
+    fn pow(&self, y: &Self) -> Self;
+
+    /// Entry point for bounded iterations for comparing two exp values.
+    fn exp_cmp(&self, max_n: u64, bound_self: i64, compare: &Self) -> ExpCmpOrdering;
+}
+
+#[derive(Debug, Clone, PartialEq)]
+pub enum ExpOrdering {
+    GT,
+    LT,
+    UNKNOWN,
+}
+
+impl From<&str> for ExpOrdering {
+    fn from(s: &str) -> Self {
+        match s {
+            "GT" => ExpOrdering::GT,
+            "LT" => ExpOrdering::LT,
+            _ => ExpOrdering::UNKNOWN,
+        }
+    }
+}
+
+#[derive(Debug, Clone, PartialEq)]
+pub struct ExpCmpOrdering {
+    pub iterations: u64,
+    pub estimation: ExpOrdering,
+    pub approx: Decimal,
 }
 
 #[cfg(test)]
 mod tests {
+    use std::fs::File;
+    use std::io::BufRead;
+    use std::path::PathBuf;
+
+    #[cfg(feature = "gmp")]
+    use crate::math_gmp::Decimal;
+    #[cfg(feature = "num")]
+    use crate::math_num::Decimal;
+
     use super::*;
 
     #[test]
-    fn test_add() {
-        assert_eq!(add(1, 2), 3);
+    fn test_fixed_precision() {
+        let fp: Decimal = Decimal::new(34);
+        assert_eq!(fp.precision(), 34);
+        assert_eq!(fp.to_string(), "0.0000000000000000000000000000000000");
+    }
+
+    #[test]
+    fn test_fixed_precision_eq() {
+        let fp1: Decimal = Decimal::new(34);
+        let fp2: Decimal = Decimal::new(34);
+        assert_eq!(fp1, fp2);
+    }
+
+    #[test]
+    fn test_fixed_precision_from_str() {
+        let fp: Decimal = Decimal::from_str("1234567890123456789012345678901234", 34).unwrap();
+        assert_eq!(fp.precision(), 34);
+        assert_eq!(fp.to_string(), "0.1234567890123456789012345678901234");
+
+        let fp: Decimal = Decimal::from_str("-1234567890123456789012345678901234", 30).unwrap();
+        assert_eq!(fp.precision(), 30);
+        assert_eq!(fp.to_string(), "-1234.567890123456789012345678901234");
+
+        let fp: Decimal = Decimal::from_str("-1234567890123456789012345678901234", 34).unwrap();
+        assert_eq!(fp.precision(), 34);
+        assert_eq!(fp.to_string(), "-0.1234567890123456789012345678901234");
+    }
+
+    #[test]
+    fn test_fixed_precision_exp() {
+        let fp: Decimal = Decimal::from(1u64);
+        assert_eq!(fp.to_string(), "1.0000000000000000000000000000000000");
+        let exp_fp = fp.exp();
+        assert_eq!(exp_fp.to_string(), "2.7182818284590452353602874043083282");
+    }
+
+    #[test]
+    fn test_fixed_precision_mul() {
+        let fp1: Decimal = Decimal::from_str("52500000000000000000000000000000000", 34).unwrap();
+        let fp2: Decimal = Decimal::from_str("43000000000000000000000000000000000", 34).unwrap();
+        let fp3 = &fp1 * &fp2;
+        assert_eq!(fp3.to_string(), "22.5750000000000000000000000000000000");
+        let fp4 = fp1 * fp2;
+        assert_eq!(fp4.to_string(), "22.5750000000000000000000000000000000");
+    }
+
+    #[test]
+    fn test_fixed_precision_div() {
+        let fp1: Decimal = Decimal::from_str("1", 34).unwrap();
+        let fp2: Decimal = Decimal::from_str("10", 34).unwrap();
+        let fp3 = &fp1 / &fp2;
+        assert_eq!(fp3.to_string(), "0.1000000000000000000000000000000000");
+        let fp4 = fp1 / fp2;
+        assert_eq!(fp4.to_string(), "0.1000000000000000000000000000000000");
+    }
+
+    #[test]
+    fn test_fixed_precision_sub() {
+        let fp1: Decimal = Decimal::from_str("1", 34).unwrap();
+        assert_eq!(fp1.to_string(), "0.0000000000000000000000000000000001");
+        let fp2: Decimal = Decimal::from_str("10", 34).unwrap();
+        assert_eq!(fp2.to_string(), "0.0000000000000000000000000000000010");
+        let fp3 = &fp1 - &fp2;
+        assert_eq!(fp3.to_string(), "-0.0000000000000000000000000000000009");
+        let fp4 = fp1 - fp2;
+        assert_eq!(fp4.to_string(), "-0.0000000000000000000000000000000009");
+    }
+
+    #[test]
+    fn golden_tests() {
+        let mut data_path = PathBuf::from(env!("CARGO_MANIFEST_DIR"));
+        data_path.push("tests/data/golden_tests.txt");
+
+        // read each line of golden_tests.txt
+        let file = File::open(data_path).expect("golden_tests.txt: file not found");
+        let reader = std::io::BufReader::new(file);
+
+        // read each line of golden_tests_result.txt
+        let mut data_path = PathBuf::from(env!("CARGO_MANIFEST_DIR"));
+        data_path.push("tests/data/golden_tests_result.txt");
+        let file = File::open(data_path).expect("golden_tests_result.txt: file not found");
+        let result_reader = std::io::BufReader::new(file);
+
+        let one: Decimal = Decimal::from(1u64);
+        let ten: Decimal = Decimal::from(10u64);
+        let f: Decimal = &one / &ten;
+        assert_eq!(f.to_string(), "0.1000000000000000000000000000000000");
+
+        for (test_line, result_line) in reader.lines().zip(result_reader.lines()) {
+            let test_line = test_line.expect("failed to read line");
+            // println!("test_line: {}", test_line);
+            let mut parts = test_line.split_whitespace();
+            let x = Decimal::from_str(parts.next().unwrap(), DEFAULT_PRECISION)
+                .expect("failed to parse x");
+            let a = Decimal::from_str(parts.next().unwrap(), DEFAULT_PRECISION)
+                .expect("failed to parse a");
+            let b = Decimal::from_str(parts.next().unwrap(), DEFAULT_PRECISION)
+                .expect("failed to parse b");
+            let result_line = result_line.expect("failed to read line");
+            // println!("result_line: {}", result_line);
+            let mut result_parts = result_line.split_whitespace();
+            let expected_exp_x = result_parts.next().expect("expected_exp_x not found");
+            let expected_ln_a = result_parts.next().expect("expected_ln_a not found");
+            let expected_threshold_b = result_parts.next().expect("expected_threshold_b not found");
+            let expected_approx_exp = result_parts.next().expect("expected_approx_exp not found");
+            let expected_estimation =
+                ExpOrdering::from(result_parts.next().expect("expected_estimation not found"));
+            let expected_iterations = result_parts.next().expect("expected_iterations not found");
+
+            // calculate exp' x
+            let exp_x = x.exp();
+            assert_eq!(exp_x.to_string(), expected_exp_x);
+
+            // calculate ln' a, print -ln' a
+            let ln_a = a.ln();
+            assert_eq!((-ln_a).to_string(), expected_ln_a);
+
+            // calculate (1 - f) *** b
+            let c = &one - &f;
+            assert_eq!(c.to_string(), "0.9000000000000000000000000000000000");
+            let threshold_b = c.pow(&b);
+            assert_eq!((&one - &threshold_b).to_string(), expected_threshold_b);
+
+            // do Taylor approximation for
+            //  a < 1 - (1 - f) *** b <=> 1/(1-a) < exp(-b * ln' (1 - f))
+            // using Lagrange error term calculation
+            let c = &one - &f;
+            let temp = c.ln();
+            let alpha = &b * &temp;
+            let alpha = -alpha;
+            let q_ = &one - &a;
+            let q = &one / &q_;
+            let res = alpha.exp_cmp(1000, 3, &q);
+
+            // println!("alpha: {}", alpha);
+            // println!("q: {}", q);
+            // println!("res.approx: {}", res.approx);
+            // println!("res.estimation: {:?}", res.estimation);
+            // println!("res.iterations: {}", res.iterations);
+
+            // we compare 1/(1-p) < e^-(1-(1-f)^sigma)
+            if a < (&one - &threshold_b) && res.estimation != ExpOrdering::LT {
+                println!(
+                    "wrong result should be leader {} should be more like {}",
+                    &temp,
+                    &one - &threshold_b
+                );
+                assert!(false);
+            }
+
+            if !(a < (&one - &threshold_b)) && res.estimation != ExpOrdering::GT {
+                println!(
+                    "wrong result should not be leader {} should be more like {}",
+                    &temp,
+                    &one - &threshold_b
+                );
+                assert!(false);
+            }
+
+            assert_eq!(res.approx.to_string(), expected_approx_exp);
+            assert_eq!(res.estimation, expected_estimation);
+            assert_eq!(res.iterations.to_string(), expected_iterations);
+        }
     }
 }

pallas-math/src/math_num.rs

@@ -0,0 +1,597 @@
+/*!
+# Cardano Math functions using the num-bigint crate
+ */
+
+use std::cmp::Ordering;
+use std::fmt::{Display, Formatter};
+use std::ops::{Div, Mul, Neg, Sub};
+use std::str::FromStr;
+
+use num_bigint::BigInt;
+use num_integer::Integer;
+use num_traits::{Signed, ToPrimitive};
+use once_cell::sync::Lazy;
+use regex::Regex;
+
+use crate::math::{Error, ExpCmpOrdering, ExpOrdering, FixedPrecision, DEFAULT_PRECISION};
+
+#[derive(Debug, Clone)]
+pub struct Decimal {
+    precision: u64,
+    precision_multiplier: BigInt,
+    data: BigInt,
+}
+
+impl PartialEq for Decimal {
+    fn eq(&self, other: &Self) -> bool {
+        self.precision == other.precision
+            && self.precision_multiplier == other.precision_multiplier
+            && self.data == other.data
+    }
+}
+
+impl PartialOrd for Decimal {
+    fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
+        if self.precision != other.precision
+            || self.precision_multiplier != other.precision_multiplier
+        {
+            return None;
+        }
+        Some(self.data.cmp(&other.data))
+    }
+}
+
+impl Display for Decimal {
+    fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result {
+        write!(
+            f,
+            "{}",
+            print_fixedp(
+                &self.data,
+                &self.precision_multiplier,
+                self.precision as usize,
+            )
+        )
+    }
+}
+
+impl From<u64> for Decimal {
+    fn from(n: u64) -> Self {
+        let mut result = Decimal::new(DEFAULT_PRECISION);
+        result.data = BigInt::from(n) * &result.precision_multiplier;
+        result
+    }
+}
+
+impl From<i64> for Decimal {
+    fn from(n: i64) -> Self {
+        let mut result = Decimal::new(DEFAULT_PRECISION);
+        result.data = BigInt::from(n) * &result.precision_multiplier;
+        result
+    }
+}
+
+impl From<&BigInt> for Decimal {
+    fn from(n: &BigInt) -> Self {
+        let mut result = Decimal::new(DEFAULT_PRECISION);
+        result.data.clone_from(n);
+        result
+    }
+}
+
+impl Neg for Decimal {
+    type Output = Self;
+
+    fn neg(self) -> Self::Output {
+        let mut result = Decimal::new(self.precision);
+        result.data = -self.data;
+        result
+    }
+}
+
+impl Mul for Decimal {
+    type Output = Self;
+
+    fn mul(self, rhs: Self) -> Self::Output {
+        let mut result = Decimal::new(self.precision);
+        result.data = &self.data * &rhs.data;
+        scale(&mut result.data);
+        result
+    }
+}
+
+// Implement Mul for a reference to Decimal
+impl<'a, 'b> Mul<&'b Decimal> for &'a Decimal {
+    type Output = Decimal;
+
+    fn mul(self, rhs: &'b Decimal) -> Self::Output {
+        let mut result = Decimal::new(self.precision);
+        result.data = &self.data * &rhs.data;
+        scale(&mut result.data);
+        result
+    }
+}
+
+impl Div for Decimal {
+    type Output = Self;
+
+    fn div(self, rhs: Self) -> Self::Output {
+        let mut result = Decimal::new(self.precision);
+        div(&mut result.data, &self.data, &rhs.data);
+        result
+    }
+}
+
+// Implement Div for a reference to Decimal
+impl<'a, 'b> Div<&'b Decimal> for &'a Decimal {
+    type Output = Decimal;
+
+    fn div(self, rhs: &'b Decimal) -> Self::Output {
+        let mut result = Decimal::new(self.precision);
+        div(&mut result.data, &self.data, &rhs.data);
+        result
+    }
+}
+
+impl Sub for Decimal {
+    type Output = Self;
+
+    fn sub(self, rhs: Self) -> Self::Output {
+        let mut result = Decimal::new(self.precision);
+        result.data = &self.data - &rhs.data;
+        result
+    }
+}
+
+// Implement Sub for a reference to Decimal
+impl<'a, 'b> Sub<&'b Decimal> for &'a Decimal {
+    type Output = Decimal;
+
+    fn sub(self, rhs: &'b Decimal) -> Self::Output {
+        let mut result = Decimal::new(self.precision);
+        result.data = &self.data - &rhs.data;
+        result
+    }
+}
+
+impl FixedPrecision for Decimal {
+    fn new(precision: u64) -> Self {
+        let ten = BigInt::from(10);
+        let precision_multiplier = ten.pow(precision as u32);
+        let data = BigInt::from(0);
+        Decimal {
+            precision,
+            precision_multiplier,
+            data,
+        }
+    }
+
+    fn from_str(s: &str, precision: u64) -> Result<Self, Error> {
+        // assert that s contains only digits using a regex
+        if !DIGITS_REGEX.is_match(s) {
+            return Err(Error::RegexFailure(regex::Error::Syntax(
+                "string contained non-digits".to_string(),
+            )));
+        }
+
+        let mut decimal = Decimal::new(precision);
+        decimal.data = BigInt::from_str(s).unwrap();
+        Ok(decimal)
+    }
+
+    fn precision(&self) -> u64 {
+        self.precision
+    }
+
+    fn exp(&self) -> Self {
+        let mut exp_x = Decimal::new(self.precision);
+        ref_exp(&mut exp_x.data, &self.data);
+        exp_x
+    }
+
+    fn ln(&self) -> Self {
+        let mut ln_x = Decimal::new(self.precision);
+        ref_ln(&mut ln_x.data, &self.data);
+        ln_x
+    }
+
+    fn pow(&self, rhs: &Self) -> Self {
+        let mut pow_x = Decimal::new(self.precision);
+        ref_pow(&mut pow_x.data, &self.data, &rhs.data);
+        pow_x
+    }
+
+    fn exp_cmp(&self, max_n: u64, bound_self: i64, compare: &Self) -> ExpCmpOrdering {
+        let mut output = Decimal::new(self.precision);
+        ref_exp_cmp(
+            &mut output.data,
+            max_n,
+            &self.data,
+            bound_self,
+            &compare.data,
+        )
+    }
+}
+
+fn print_fixedp(n: &BigInt, precision: &BigInt, width: usize) -> String {
+    let (mut temp_q, mut temp_r) = n.div_rem(precision);
+
+    let is_negative_q = temp_q < ZERO.value;
+    let is_negative_r = temp_r < ZERO.value;
+
+    if is_negative_q {
+        temp_q = temp_q.abs();
+    }
+    if is_negative_r {
+        temp_r = temp_r.abs();
+    }
+
+    let mut s = String::new();
+    if is_negative_q || is_negative_r {
+        s.push('-');
+    }
+    s.push_str(&temp_q.to_string());
+    s.push('.');
+    let r = temp_r.to_string();
+    let r_len = r.len();
+    // fill with zeroes up to width for the fractional part
+    if r_len < width {
+        s.push_str(&"0".repeat(width - r_len));
+    }
+    s.push_str(&r);
+    s
+}
+
+struct Constant {
+    value: BigInt,
+}
+
+impl Constant {
+    pub fn new(init: fn() -> BigInt) -> Constant {
+        Constant { value: init() }
+    }
+}
+
+unsafe impl Sync for Constant {}
+unsafe impl Send for Constant {}
+
+static DIGITS_REGEX: Lazy<Regex> = Lazy::new(|| Regex::new(r"^-?\d+$").unwrap());
+static TEN: Lazy<Constant> = Lazy::new(|| Constant::new(|| BigInt::from(10)));
+static PRECISION: Lazy<Constant> = Lazy::new(|| Constant::new(|| TEN.value.pow(34)));
+static EPS: Lazy<Constant> = Lazy::new(|| Constant::new(|| TEN.value.pow(34 - 24)));
+static ONE: Lazy<Constant> = Lazy::new(|| Constant::new(|| BigInt::from(1) * &PRECISION.value));
+static ZERO: Lazy<Constant> = Lazy::new(|| Constant::new(|| BigInt::from(0)));
+static E: Lazy<Constant> = Lazy::new(|| {
+    Constant::new(|| {
+        let mut e = BigInt::from(0);
+        ref_exp(&mut e, &ONE.value);
+        e
+    })
+});
+
+/// Entry point for 'exp' approximation. First does the scaling of 'x' to [0,1]
+/// and then calls the continued fraction approximation function.
+fn ref_exp(rop: &mut BigInt, x: &BigInt) -> i32 {
+    let mut iterations = 0;
+    match x.cmp(&ZERO.value) {
+        std::cmp::Ordering::Equal => {
+            // rop = 1
+            rop.clone_from(&ONE.value);
+        }
+        std::cmp::Ordering::Less => {
+            let x_ = -x;
+            let mut temp = BigInt::from(0);
+            iterations = ref_exp(&mut temp, &x_);
+            // rop = 1 / temp
+            div(rop, &ONE.value, &temp);
+        }
+        std::cmp::Ordering::Greater => {
+            let mut n_exponent = x.div_ceil(&PRECISION.value);
+            let n = n_exponent.to_u32().expect("n_exponent to_u32 failed");
+            n_exponent *= &PRECISION.value; /* ceil(x) */
+            let x_ = x / n;
+            iterations = mp_exp_taylor(rop, 1000, &x_, &EPS.value);
+
+            // rop = rop.pow(n)
+            ipow(rop, &rop.clone(), n as i64);
+        }
+    }
+
+    iterations
+}
+
+/// Division with quotent and remainder
+#[inline]
+fn div_qr(q: &mut BigInt, r: &mut BigInt, x: &BigInt, y: &BigInt) {
+    (*q, *r) = x.div_rem(y);
+}
+
+/// Division
+pub fn div(rop: &mut BigInt, x: &BigInt, y: &BigInt) {
+    let mut temp_q = BigInt::from(0);
+    let mut temp_r = BigInt::from(0);
+    let mut temp: BigInt;
+    div_qr(&mut temp_q, &mut temp_r, x, y);
+
+    temp = &temp_q * &PRECISION.value;
+    temp_r = &temp_r * &PRECISION.value;
+    let temp_r2 = temp_r.clone();
+    div_qr(&mut temp_q, &mut temp_r, &temp_r2, y);
+
+    temp += &temp_q;
+    *rop = temp;
+}
+/// Taylor / MacLaurin series approximation
+fn mp_exp_taylor(rop: &mut BigInt, max_n: i32, x: &BigInt, epsilon: &BigInt) -> i32 {
+    let mut divisor = ONE.value.clone();
+    let mut last_x = ONE.value.clone();
+    rop.clone_from(&ONE.value);
+    let mut n = 0;
+    while n < max_n {
+        let mut next_x = x * &last_x;
+        scale(&mut next_x);
+        let next_x2 = next_x.clone();
+        div(&mut next_x, &next_x2, &divisor);
+
+        if next_x.abs() < epsilon.abs() {
+            break;
+        }
+
+        divisor += &ONE.value;
+        *rop += &next_x;
+        last_x.clone_from(&next_x);
+        n += 1;
+    }
+
+    n
+}
+
+fn scale(rop: &mut BigInt) {
+    let mut temp = BigInt::from(0);
+    let mut a = BigInt::from(0);
+    div_qr(&mut a, &mut temp, rop, &PRECISION.value);
+    if *rop < ZERO.value && temp != ZERO.value {
+        a -= 1;
+    }
+    *rop = a;
+}
+
+/// Integer power internal function
+fn ipow_(rop: &mut BigInt, x: &BigInt, n: i64) {
+    if n == 0 {
+        rop.clone_from(&ONE.value);
+    } else if n % 2 == 0 {
+        let mut res = BigInt::from(0);
+        ipow_(&mut res, x, n / 2);
+        *rop = &res * &res;
+        scale(rop);
+    } else {
+        let mut res = BigInt::from(0);
+        ipow_(&mut res, x, n - 1);
+        *rop = res * x;
+        scale(rop);
+    }
+}
+
+/// Integer power
+fn ipow(rop: &mut BigInt, x: &BigInt, n: i64) {
+    if n < 0 {
+        let mut temp = BigInt::from(0);
+        ipow_(&mut temp, x, -n);
+        div(rop, &ONE.value, &temp);
+    } else {
+        ipow_(rop, x, n);
+    }
+}
+
+/// Compute an approximation of 'ln(1 + x)' via continued fractions. Either for a
+///    maximum of 'maxN' iterations or until the absolute difference between two
+///    succeeding convergents is smaller than 'eps'. Assumes 'x' to be within
+///    [1,e).
+fn mp_ln_n(rop: &mut BigInt, max_n: i32, x: &BigInt, epsilon: &BigInt) {
+    let mut ba: BigInt;
+    let mut aa: BigInt;
+    let mut ab: BigInt;
+    let mut bb: BigInt;
+    let mut a_: BigInt;
+    let mut b_: BigInt;
+    let mut diff: BigInt;
+    let mut convergent: BigInt = BigInt::from(0);
+    let mut last: BigInt = BigInt::from(0);
+    let mut first = true;
+    let mut n = 1;
+
+    let mut a: BigInt;
+    let mut b = ONE.value.clone();
+
+    let mut an_m2 = ONE.value.clone();
+    let mut bn_m2 = BigInt::from(0);
+    let mut an_m1 = BigInt::from(0);
+    let mut bn_m1 = ONE.value.clone();
+
+    let mut curr_a = 1;
+
+    while n <= max_n + 2 {
+        let curr_a_2 = curr_a * curr_a;
+        a = x * curr_a_2;
+        if n > 1 && n % 2 == 1 {
+            curr_a += 1;
+        }
+
+        ba = &b * &an_m1;
+        scale(&mut ba);
+        aa = &a * &an_m2;
+        scale(&mut aa);
+        a_ = &ba + &aa;
+
+        bb = &b * &bn_m1;
+        scale(&mut bb);
+        ab = &a * &bn_m2;
+        scale(&mut ab);
+        b_ = &bb + &ab;
+
+        div(&mut convergent, &a_, &b_);
+
+        if first {
+            first = false;
+        } else {
+            diff = &convergent - &last;
+            if diff.abs() < epsilon.abs() {
+                break;
+            }
+        }
+
+        last.clone_from(&convergent);
+
+        n += 1;
+        an_m2.clone_from(&an_m1);
+        bn_m2.clone_from(&bn_m1);
+        an_m1.clone_from(&a_);
+        bn_m1.clone_from(&b_);
+
+        b += &ONE.value;
+    }
+
+    *rop = convergent;
+}
+
+fn find_e(x: &BigInt) -> i64 {
+    let mut x_: BigInt = BigInt::from(0);
+    let mut x__: BigInt;
+
+    div(&mut x_, &ONE.value, &E.value);
+    x__ = E.value.clone();
+
+    let mut l = -1;
+    let mut u = 1;
+    while &x_ > x || &x__ < x {
+        x_ = &x_ * &x_;
+        scale(&mut x_);
+
+        x__ = &x__ * &x__;
+        scale(&mut x__);
+
+        l *= 2;
+        u *= 2;
+    }
+
+    while l + 1 != u {
+        let mid = l + ((u - l) / 2);
+
+        ipow(&mut x_, &E.value, mid);
+        if x < &x_ {
+            u = mid;
+        } else {
+            l = mid;
+        }
+    }
+    l
+}
+
+/// Entry point for 'ln' approximation. First does the necessary scaling, and
+/// then calls the continued fraction calculation. For any value outside the
+/// domain, i.e., 'x in (-inf,0]', the function returns '-INFINITY'.
+fn ref_ln(rop: &mut BigInt, x: &BigInt) -> bool {
+    let mut factor = BigInt::from(0);
+    let mut x_ = BigInt::from(0);
+    if x <= &ZERO.value {
+        return false;
+    }
+
+    let n = find_e(x);
+
+    *rop = BigInt::from(n);
+    *rop = rop.clone() * &PRECISION.value;
+    ref_exp(&mut factor, rop);
+
+    div(&mut x_, x, &factor);
+
+    x_ = &x_ - &ONE.value;
+
+    let x_2 = x_.clone();
+    mp_ln_n(&mut x_, 1000, &x_2, &EPS.value);
+    *rop = rop.clone() + &x_;
+
+    true
+}
+
+fn ref_pow(rop: &mut BigInt, base: &BigInt, exponent: &BigInt) {
+    /* x^y = exp(y * ln x) */
+    let mut tmp: BigInt = BigInt::from(0);
+    ref_ln(&mut tmp, base);
+    tmp *= exponent;
+    scale(&mut tmp);
+    ref_exp(rop, &tmp);
+}
+
+/// `bound_x` is the bound for exp in the interval x is chosen from
+/// `compare` the value to compare to
+///
+/// if the result is GT, then the computed value is guaranteed to be greater, if
+/// the result is LT, the computed value is guaranteed to be less than
+/// `compare`. In the case of `UNKNOWN` no conclusion was possible for the
+/// selected precision.
+///
+/// Lagrange remainder require knowledge of the maximum value to compute the
+/// maximal error of the remainder.
+fn ref_exp_cmp(
+    rop: &mut BigInt,
+    max_n: u64,
+    x: &BigInt,
+    bound_x: i64,
+    compare: &BigInt,
+) -> ExpCmpOrdering {
+    rop.clone_from(&ONE.value);
+    let mut n = 0u64;
+    let mut divisor: BigInt;
+    let mut next_x: BigInt;
+    let mut error: BigInt;
+    let mut upper: BigInt;
+    let mut lower: BigInt;
+    let mut error_term: BigInt;
+
+    divisor = ONE.value.clone();
+    error = x.clone();
+
+    let mut estimate = ExpOrdering::UNKNOWN;
+    while n < max_n {
+        next_x = error.clone();
+        if next_x.abs() < EPS.value.abs() {
+            break;
+        }
+        divisor += &ONE.value;
+
+        // update error estimation, this is initially bound_x * x and in general
+        // bound_x * x^(n+1)/(n + 1)!  we use `error` to store the x^n part and a
+        // single integral multiplication with the bound
+        error *= x;
+        scale(&mut error);
+        let e2 = error.clone();
+        div(&mut error, &e2, &divisor);
+        error_term = &error * bound_x;
+        *rop += &next_x;
+
+        /* compare is guaranteed to be above overall result */
+        upper = &*rop + &error_term;
+        if compare > &upper {
+            estimate = ExpOrdering::GT;
+            n += 1;
+            break;
+        }
+
+        /* compare is guaranteed to be below overall result */
+        lower = &*rop - &error_term;
+        if compare < &lower {
+            estimate = ExpOrdering::LT;
+            n += 1;
+            break;
+        }
+        n += 1;
+    }
+
+    ExpCmpOrdering {
+        iterations: n,
+        estimation: estimate,
+        approx: Decimal::from(&*rop),
+    }
+}