diff --git a/packages/injective-math/src/fp_decimal/exp.rs b/packages/injective-math/src/fp_decimal/exp.rs index a43ad5e5..f1e42afd 100644 --- a/packages/injective-math/src/fp_decimal/exp.rs +++ b/packages/injective-math/src/fp_decimal/exp.rs @@ -7,23 +7,51 @@ use crate::fp_decimal::{FPDecimal, U256}; impl FPDecimal { #[allow(clippy::many_single_char_names)] - pub fn _exp_taylor_expansion(a: FPDecimal, b: FPDecimal, n: u128) -> FPDecimal { + pub fn _exp_taylor_expansion(a: FPDecimal, b: FPDecimal) -> FPDecimal { //a^b n+1 terms taylor expansion - assert!(n <= 13u128); - if n == 0 { - FPDecimal::ONE - } else { - let base = a.ln() * b; - let mut x = FPDecimal::ONE + base; - let mut numerator = base; - let mut denominator = FPDecimal::ONE; - for i in 2..n + 1 { + let base = a.ln() * b; + let mut numerator = base; + let mut denominator = FPDecimal::ONE; + + let denominator_parts = vec![ + FPDecimal::TWO, + FPDecimal::THREE, + FPDecimal::FOUR, + FPDecimal::FIVE, + FPDecimal::SIX, + FPDecimal::SEVEN, + FPDecimal::EIGHT, + FPDecimal::NINE, + FPDecimal::TEN, + FPDecimal::must_from_str("11"), + FPDecimal::must_from_str("12"), + FPDecimal::must_from_str("13"), + FPDecimal::must_from_str("14"), + FPDecimal::must_from_str("15"), + FPDecimal::must_from_str("16"), + FPDecimal::must_from_str("17"), + FPDecimal::must_from_str("18"), + FPDecimal::must_from_str("19"), + FPDecimal::must_from_str("20"), + FPDecimal::must_from_str("21"), + FPDecimal::must_from_str("22"), + FPDecimal::must_from_str("23"), + FPDecimal::must_from_str("24"), + FPDecimal::must_from_str("25"), + ]; + + denominator_parts + .iter() + .map(|part| { numerator *= base; - denominator *= FPDecimal::from(i); - x += numerator / denominator; - } - x - } + denominator *= *part; + numerator / denominator + }) + .collect::>() + .iter() + .sum::() + + FPDecimal::ONE + + a.ln() * b } // e^(a) pub fn _exp(a: FPDecimal) -> FPDecimal { @@ -416,7 +444,7 @@ impl FPDecimal { // base^exponent // NOTE: only accurate for 1,3,5,7,11, and combinations of these numbers // 14 terms taylor expansion provides a good enough approximation - const N_TERMS: u128 = 13; + //const N_TERMS: u128 = 13; match exponent.cmp(&FPDecimal::ZERO) { Ordering::Equal => Ok(FPDecimal::ONE), // Negative exponent @@ -424,25 +452,25 @@ impl FPDecimal { exponent = -exponent; match exponent.cmp(&(FPDecimal::ONE)) { Ordering::Equal => Ok(FPDecimal::ONE / base), - Ordering::Less => compute_negative_exponent_less_one(base, exponent, N_TERMS), - Ordering::Greater => compute_negative_exponent_greater_one(base, exponent, N_TERMS), + Ordering::Less => compute_negative_exponent_less_one(base, exponent), + Ordering::Greater => compute_negative_exponent_greater_one(base, exponent), } } // Positive exponent Ordering::Greater => match exponent.cmp(&FPDecimal::ONE) { Ordering::Equal => Ok(base), - Ordering::Less => compute_positive_exponent_less_one(base, exponent, N_TERMS), - Ordering::Greater => compute_positive_exponent_greater_one(base, exponent, N_TERMS), + Ordering::Less => compute_positive_exponent_less_one(base, exponent), + Ordering::Greater => compute_positive_exponent_greater_one(base, exponent), }, } } - fn compute_negative_exponent_greater_one(mut base: FPDecimal, exponent: FPDecimal, n_terms: u128) -> Result { + fn compute_negative_exponent_greater_one(mut base: FPDecimal, exponent: FPDecimal) -> Result { let mut int_b = exponent.int(); let rem_b = exponent - int_b; let mut float_exp = FPDecimal::ONE; if rem_b != FPDecimal::ZERO { - float_exp = FPDecimal::_exp_taylor_expansion(FPDecimal::ONE / base, rem_b, n_terms); + float_exp = FPDecimal::_exp_taylor_expansion(FPDecimal::ONE / base, rem_b); } let mut tmp_a = FPDecimal::ONE; while int_b > FPDecimal::ONE { @@ -487,7 +515,7 @@ impl FPDecimal { type BaseCheckFunction<'a> = (&'a dyn Fn(FPDecimal) -> Option, FPDecimal); - fn compute_negative_exponent_less_one(base: FPDecimal, exponent: FPDecimal, n_terms: u128) -> Result { + fn compute_negative_exponent_less_one(base: FPDecimal, exponent: FPDecimal) -> Result { // NOTE: only accurate for 1,3,5,7,11, and combinations of these numbers let abs_error = FPDecimal::must_from_str("0.00000000000000001"); let reciprocal = if (FPDecimal::reciprocal(exponent) - FPDecimal::reciprocal(exponent).int()).abs() < abs_error { @@ -519,7 +547,7 @@ impl FPDecimal { } } - Ok(FPDecimal::_exp_taylor_expansion(FPDecimal::ONE / base, exponent, n_terms)) + Ok(FPDecimal::_exp_taylor_expansion(FPDecimal::ONE / base, exponent)) } fn positive_exponent_check_basic_log( @@ -551,7 +579,7 @@ impl FPDecimal { None } - fn compute_positive_exponent_less_one(base: FPDecimal, exponent: FPDecimal, n_terms: u128) -> Result { + fn compute_positive_exponent_less_one(base: FPDecimal, exponent: FPDecimal) -> Result { // taylor expansion approximation of exponentiation computation with float number exponent // NOTE: only accurate for 1,3,5,7,11, and combinations of these numbers let abs_error = FPDecimal::must_from_str("0.00000000000000001"); @@ -578,7 +606,7 @@ impl FPDecimal { } } - Ok(FPDecimal::_exp_taylor_expansion(base, exponent, n_terms)) + Ok(FPDecimal::_exp_taylor_expansion(base, exponent)) } fn compute_integer_exponentiation(mut base: FPDecimal, mut exponent: FPDecimal) -> FPDecimal { @@ -597,12 +625,12 @@ impl FPDecimal { base * temp_base } - fn compute_positive_exponent_greater_one(mut base: FPDecimal, exponent: FPDecimal, n_terms: u128) -> Result { + fn compute_positive_exponent_greater_one(mut base: FPDecimal, exponent: FPDecimal) -> Result { let integer_part_of_exponent = exponent.int(); let fractional_part_of_exponent = exponent - integer_part_of_exponent; let fractional_exponentiation = if fractional_part_of_exponent != FPDecimal::ZERO { - FPDecimal::_exp_taylor_expansion(base, fractional_part_of_exponent, n_terms) + FPDecimal::_exp_taylor_expansion(base, fractional_part_of_exponent) } else { FPDecimal::ONE }; @@ -936,20 +964,19 @@ impl FPDecimal { fn compute_exponentiation(base: FPDecimal, mut exponent: FPDecimal) -> Result { // a^b // 14 terms taylor expansion provides a good enough approximation - const N_TERMS: u128 = 13; match exponent.cmp(&FPDecimal::ZERO) { Ordering::Equal => Ok(FPDecimal::ONE), Ordering::Less => { exponent = -exponent; match exponent.cmp(&(FPDecimal::ONE)) { Ordering::Equal => Ok(FPDecimal::ONE / base), - Ordering::Less => compute_negative_exponent_less_one(FPDecimal::ONE / base, exponent, N_TERMS), - Ordering::Greater => compute_negative_exponent_greater_one(FPDecimal::ONE / base, exponent, N_TERMS), + Ordering::Less => compute_negative_exponent_less_one(FPDecimal::ONE / base, exponent), + Ordering::Greater => compute_negative_exponent_greater_one(FPDecimal::ONE / base, exponent), } } Ordering::Greater => match exponent.cmp(&FPDecimal::ONE) { Ordering::Equal => Ok(base), - Ordering::Less => compute_positive_exponent_less_one(base, exponent, N_TERMS), + Ordering::Less => compute_positive_exponent_less_one(base, exponent), Ordering::Greater => compute_positive_exponent_greater_one(base, exponent), }, } @@ -957,7 +984,7 @@ impl FPDecimal { type BaseCheckFunction<'a> = (&'a dyn Fn(FPDecimal) -> Option, FPDecimal); - fn compute_negative_exponent_less_one(mut base: FPDecimal, exponent: FPDecimal, n_terms: u128) -> Result { + fn compute_negative_exponent_less_one(mut base: FPDecimal, exponent: FPDecimal) -> Result { // NOTE: only accurate for 1,3,5,7,11, and combinations of these numbers base = -base; @@ -981,7 +1008,7 @@ impl FPDecimal { return Ok(-FPDecimal::ONE / base); } - Ok(FPDecimal::_exp_taylor_expansion(FPDecimal::ONE / base, exponent, n_terms)) + Ok(FPDecimal::_exp_taylor_expansion(FPDecimal::ONE / base, exponent)) } fn check_conditions_and_return_negative(mut base: FPDecimal, divisor: FPDecimal, exponent: &FPDecimal) -> bool { @@ -1008,12 +1035,12 @@ impl FPDecimal { } } - fn compute_negative_exponent_greater_one(mut base: FPDecimal, exponent: FPDecimal, n_terms: u128) -> Result { + fn compute_negative_exponent_greater_one(mut base: FPDecimal, exponent: FPDecimal) -> Result { let integer_part_of_exponent = exponent.int(); let fractional_part_of_exponent = exponent - integer_part_of_exponent; let fractional_exponentiation = if fractional_part_of_exponent != FPDecimal::ZERO { - FPDecimal::_exp_taylor_expansion(FPDecimal::ONE / base, fractional_part_of_exponent, n_terms) + FPDecimal::_exp_taylor_expansion(FPDecimal::ONE / base, fractional_part_of_exponent) } else { FPDecimal::ONE }; @@ -1021,7 +1048,7 @@ impl FPDecimal { Ok(FPDecimal::ONE / base * fractional_exponentiation) } - fn compute_positive_exponent_less_one(mut base: FPDecimal, exponent: FPDecimal, n_terms: u128) -> Result { + fn compute_positive_exponent_less_one(mut base: FPDecimal, exponent: FPDecimal) -> Result { // taylor expansion approximation of exponentiation computation with float number exponent // NOTE: only accurate for 1,3,5,7,11, and combinations of these numbers base = -base; @@ -1043,7 +1070,7 @@ impl FPDecimal { } } - Ok(FPDecimal::_exp_taylor_expansion(base, exponent, n_terms)) + Ok(FPDecimal::_exp_taylor_expansion(base, exponent)) } fn check_conditions_and_return(base: &mut FPDecimal, divisor: &FPDecimal, exponent: &FPDecimal) -> bool { @@ -1111,10 +1138,20 @@ impl FPDecimal { #[cfg(test)] mod tests { + use crate::fp_decimal::U256; use crate::FPDecimal; - use bigint::U256; use std::str::FromStr; + #[test] + fn test_3_pow_2_point_3() { + // a^x = e^(xln(a)) + // 3^2.3 = e(2.3ln(3)) + assert_eq!( + FPDecimal::_exp_taylor_expansion(FPDecimal::THREE, FPDecimal::must_from_str("2.3")), + FPDecimal::must_from_str("12.513502532843181622") + ); + } + #[test] fn test_exp() { assert_eq!(FPDecimal::_exp(FPDecimal::ONE), FPDecimal::E); diff --git a/packages/injective-math/src/fp_decimal/from_str.rs b/packages/injective-math/src/fp_decimal/from_str.rs index 29c25525..e2c4aa7c 100644 --- a/packages/injective-math/src/fp_decimal/from_str.rs +++ b/packages/injective-math/src/fp_decimal/from_str.rs @@ -1,9 +1,7 @@ -use bigint::U256; +use crate::fp_decimal::{FPDecimal, U256}; use cosmwasm_std::StdError; use std::str::FromStr; -use crate::fp_decimal::FPDecimal; - impl FromStr for FPDecimal { type Err = StdError; @@ -14,11 +12,14 @@ impl FromStr for FPDecimal { /// This never performs any kind of rounding. /// More than 18 fractional digits, even zeros, result in an error. fn from_str(input: &str) -> Result { - let sign = if input.starts_with('-') { 0 } else { 1 }; - let parts: Vec<&str> = input.trim_start_matches('-').split('.').collect(); + let mut sign = if input.starts_with('-') { 0 } else { 1 }; + let parts = input.trim_start_matches('-').split('.').collect::>(); match parts.len() { 1 => { let integer = U256::from_dec_str(parts[0]).map_err(|_| StdError::generic_err("Error parsing integer"))?; + if integer == U256([0, 0, 0, 0]) { + sign = 1; + } Ok(FPDecimal { num: integer * FPDecimal::ONE.num, sign, @@ -30,7 +31,9 @@ impl FromStr for FPDecimal { let exp = FPDecimal::DIGITS .checked_sub(parts[1].len()) .ok_or_else(|| StdError::generic_err(format!("Cannot parse more than {} fractional digits", FPDecimal::DIGITS)))?; - + if integer == U256([0, 0, 0, 0]) { + sign = 1; + } Ok(FPDecimal { num: integer * FPDecimal::ONE.num + fraction * U256::exp10(exp), sign, @@ -47,20 +50,31 @@ impl FPDecimal { pub fn must_from_str(input: &str) -> Self { let i = Self::from_str(input).unwrap(); // to handle must_from_str("-0") - if i.num == U256([0, 0, 0, 0]) { - return FPDecimal::ZERO; - } + //if i.num == U256([0, 0, 0, 0]) { + // return FPDecimal::ZERO; + //} i } } #[cfg(test)] mod tests { - use crate::FPDecimal; - use bigint::U256; + use primitive_types::U256; use std::str::FromStr; + #[test] + fn test_from_str_zero() { + let zero = FPDecimal::from_str("0").unwrap(); + let neg_zero = FPDecimal::from_str("-0").unwrap(); + let zero_zero = FPDecimal::from_str("00").unwrap(); + let neg_zero_zero = FPDecimal::from_str("-00").unwrap(); + assert_eq!(zero, FPDecimal::ZERO); + assert_eq!(zero_zero, FPDecimal::ZERO); + assert_eq!(neg_zero, FPDecimal::ZERO); + assert_eq!(neg_zero_zero, FPDecimal::ZERO); + } + #[test] fn test_from_str() { let val = FPDecimal::from_str("-1.23"); diff --git a/packages/injective-math/src/fp_decimal/log.rs b/packages/injective-math/src/fp_decimal/log.rs index 3d389b3b..45581560 100644 --- a/packages/injective-math/src/fp_decimal/log.rs +++ b/packages/injective-math/src/fp_decimal/log.rs @@ -518,8 +518,47 @@ impl FPDecimal { } /// natural logarithm + + #[allow(clippy::many_single_char_names)] + pub fn _ln_greater_than_one(mut a: FPDecimal) -> FPDecimal { + // ln(123.456) => ln(1.23456 * 10^2) =>ln(1.23456)+2ln(10) + let mut exponent = 0u128; + while a > FPDecimal::ONE { + a /= FPDecimal::TEN; + exponent += 1; + } + a.ln() + FPDecimal::from(exponent) * FPDecimal::TEN.ln() + } + + #[allow(clippy::many_single_char_names)] + fn two_agm(mut a0: FPDecimal, mut b0: FPDecimal, tol: FPDecimal) -> FPDecimal { + loop { + if (a0 - b0).abs() > tol { + break; + } + let a1 = (a0 + b0) / FPDecimal::TWO; + let b1 = (a0 * b0).sqrt(); + a0 = a1; + b0 = b1; + } + a0 + b0 + } + + #[allow(clippy::many_single_char_names)] + fn _ln(&self) -> FPDecimal { + // m =8, 2**8=256; + let two_pow_m = FPDecimal::from(256u128); + let s = *self * two_pow_m; + let tol = FPDecimal::must_from_str("0.00000001"); + let a0 = FPDecimal::ONE; + let b0 = FPDecimal::FOUR / s; + let two_agm = FPDecimal::two_agm(a0, b0, tol); + + FPDecimal::PI / two_agm - FPDecimal::EIGHT * FPDecimal::LN2 + } + #[allow(clippy::many_single_char_names)] - pub fn _ln(a: FPDecimal) -> FPDecimal { + pub fn _ln_old(a: FPDecimal) -> FPDecimal { assert!(a.sign != 0); assert!(a != FPDecimal::ZERO); let mut v = a.num; @@ -589,10 +628,13 @@ impl FPDecimal { } pub fn ln(&self) -> FPDecimal { + if *self == FPDecimal::TWO { + return FPDecimal::LN2; + } if let Some(value) = self._log_e() { return value; } - FPDecimal::_ln(*self) + self._ln() } pub fn log(&self, base: FPDecimal) -> FPDecimal { @@ -626,14 +668,20 @@ impl FPDecimal { mod tests { use crate::FPDecimal; - use bigint::U256; + //use bigint::U256; + + use primitive_types::U256; + #[test] + fn test_ln3() { + assert_eq!(FPDecimal::THREE.ln(), FPDecimal::must_from_str("1.09861228866810969")); + } #[test] fn test_ln_sanity() { - let half = FPDecimal::ONE.div(2i128); - println!("{}", FPDecimal::_ln(half)); // works if you comment this out - let num = FPDecimal::ONE.mul(5).div(4); - println!("{}", FPDecimal::checked_positive_pow(num, half).unwrap()); + //let half = FPDecimal::ONE.div(2i128); + assert_eq!((FPDecimal::ONE / FPDecimal::TWO).ln(), FPDecimal::must_from_str("-0.693147180559945307")); + assert_eq!((FPDecimal::FIVE / FPDecimal::FOUR).ln(), FPDecimal::must_from_str("0.223143551314209761")); + assert_eq!((FPDecimal::TEN.pow(FPDecimal::must_from_str("20"))).ln(), FPDecimal::must_from_str("20")); } #[test] diff --git a/packages/injective-math/src/fp_decimal/mod.rs b/packages/injective-math/src/fp_decimal/mod.rs index bb2ffe79..022d1923 100644 --- a/packages/injective-math/src/fp_decimal/mod.rs +++ b/packages/injective-math/src/fp_decimal/mod.rs @@ -1,8 +1,8 @@ use std::str::FromStr; use std::{convert::TryFrom, ops::Neg}; -use bigint::U256; use cosmwasm_std::{Decimal256, StdError, Uint128, Uint256}; +use primitive_types::U256; use schemars::JsonSchema; #[allow(clippy::upper_case_acronyms)] @@ -122,6 +122,7 @@ impl Neg for FPDecimal { } } +// constants impl FPDecimal { pub const MAX: FPDecimal = FPDecimal { num: U256::MAX, sign: 1 }; pub const MIN: FPDecimal = FPDecimal { num: U256::MAX, sign: 0 }; @@ -134,6 +135,12 @@ impl FPDecimal { num: U256([0, 0, 0, 0]), sign: 1, }; + + pub const LN2: FPDecimal = FPDecimal { + num: U256([693_147_180_559_945_309, 0, 0, 0]), + sign: 1, + }; + pub const ONE: FPDecimal = FPDecimal { num: U256([1_000_000_000_000_000_000, 0, 0, 0]), sign: 1, @@ -196,22 +203,22 @@ impl FPDecimal { }; pub const E_10: FPDecimal = FPDecimal { - num: U256([1053370797511854089u64, 1194u64, 0, 0]), + num: U256([1_053_370_797_511_854_089u64, 1194u64, 0, 0]), sign: 1, }; // e^10 pub const LN_1_5: FPDecimal = FPDecimal { - num: U256([405465108108164382, 0, 0, 0]), + num: U256([405_465_108_108_164_382, 0, 0, 0]), sign: 1, }; // ln(1.5) pub const LN_10: FPDecimal = FPDecimal { - num: U256([2302585092994045684, 0, 0, 0]), + num: U256([2_302_585_092_994_045_684, 0, 0, 0]), sign: 1, }; // ln(10) pub const PI: FPDecimal = FPDecimal { - num: U256([3_141_592_653_589_793_238, 0, 0, 0]), + num: U256([3_141_592_653_589_793_115, 0, 0, 0]), sign: 1, }; @@ -286,9 +293,28 @@ mod trigonometry; mod tests { use std::{convert::TryFrom, str::FromStr}; + use crate::fp_decimal::U256; use crate::FPDecimal; - use bigint::U256; use cosmwasm_std::{Decimal256, Uint256}; + + #[test] + fn test_const_pi() { + let pi = FPDecimal::PI; + let three_point_two = FPDecimal::must_from_str("3.2"); + let three_point_one = FPDecimal::must_from_str("3.1"); + let pi_precise = FPDecimal::must_from_str("3.141592653589793115"); + assert!(three_point_one < pi); + assert!(pi < three_point_two); + assert_eq!(pi, pi_precise); + } + + #[test] + fn test_const_ln2() { + let ln2 = FPDecimal::LN2; + let ln2_precise = FPDecimal::must_from_str("0.693147180559945309"); + assert_eq!(ln2, ln2_precise); + } + #[test] fn test_neg_sign() { let lhs = FPDecimal::ZERO - FPDecimal::ONE; diff --git a/packages/injective-math/src/utils.rs b/packages/injective-math/src/utils.rs index c33582a2..cba94ba4 100644 --- a/packages/injective-math/src/utils.rs +++ b/packages/injective-math/src/utils.rs @@ -1,5 +1,6 @@ use crate::FPDecimal; -use bigint::U256; +use primitive_types::U256; + use cosmwasm_std::StdError; use std::cmp::Ordering; use std::{fmt::Display, str::FromStr};