From f79e0f5efff681eca6e23d5313f23c3a5d92ffdb Mon Sep 17 00:00:00 2001 From: raphaelDkhn Date: Mon, 23 Oct 2023 10:35:23 +0300 Subject: [PATCH] implement fp8x23wide --- src/numbers/fixed_point/implementations.cairo | 3 +- .../implementations/fp8x23wide.cairo | 4 + .../implementations/fp8x23wide/core.cairo | 378 +++++ .../implementations/fp8x23wide/helpers.cairo | 39 + .../implementations/fp8x23wide/math.cairo | 5 + .../fp8x23wide/math/comp.cairo | 76 + .../fp8x23wide/math/core.cairo | 660 +++++++++ .../implementations/fp8x23wide/math/hyp.cairo | 159 +++ .../implementations/fp8x23wide/math/lut.cairo | 1229 +++++++++++++++++ .../fp8x23wide/math/trig.cairo | 448 ++++++ 10 files changed, 3000 insertions(+), 1 deletion(-) create mode 100644 src/numbers/fixed_point/implementations/fp8x23wide.cairo create mode 100644 src/numbers/fixed_point/implementations/fp8x23wide/core.cairo create mode 100644 src/numbers/fixed_point/implementations/fp8x23wide/helpers.cairo create mode 100644 src/numbers/fixed_point/implementations/fp8x23wide/math.cairo create mode 100644 src/numbers/fixed_point/implementations/fp8x23wide/math/comp.cairo create mode 100644 src/numbers/fixed_point/implementations/fp8x23wide/math/core.cairo create mode 100644 src/numbers/fixed_point/implementations/fp8x23wide/math/hyp.cairo create mode 100644 src/numbers/fixed_point/implementations/fp8x23wide/math/lut.cairo create mode 100644 src/numbers/fixed_point/implementations/fp8x23wide/math/trig.cairo diff --git a/src/numbers/fixed_point/implementations.cairo b/src/numbers/fixed_point/implementations.cairo index e6152e25a..d7617f9c8 100644 --- a/src/numbers/fixed_point/implementations.cairo +++ b/src/numbers/fixed_point/implementations.cairo @@ -2,4 +2,5 @@ mod fp8x23; mod fp16x16; mod fp64x64; mod fp32x32; -mod fp16x16wide; \ No newline at end of file +mod fp16x16wide; +mod fp8x23wide; \ No newline at end of file diff --git a/src/numbers/fixed_point/implementations/fp8x23wide.cairo b/src/numbers/fixed_point/implementations/fp8x23wide.cairo new file mode 100644 index 000000000..2cc1d5085 --- /dev/null +++ b/src/numbers/fixed_point/implementations/fp8x23wide.cairo @@ -0,0 +1,4 @@ +mod core; +mod math; +mod helpers; + diff --git a/src/numbers/fixed_point/implementations/fp8x23wide/core.cairo b/src/numbers/fixed_point/implementations/fp8x23wide/core.cairo new file mode 100644 index 000000000..36b64ce5e --- /dev/null +++ b/src/numbers/fixed_point/implementations/fp8x23wide/core.cairo @@ -0,0 +1,378 @@ +use debug::PrintTrait; + +use option::OptionTrait; +use result::{ResultTrait, ResultTraitImpl}; +use traits::{TryInto, Into}; + +use orion::numbers::signed_integer::{i32::i32, i8::i8}; +use orion::numbers::fixed_point::core::{FixedTrait}; +use orion::numbers::fixed_point::implementations::fp8x23wide::math::{core, trig, hyp}; +use orion::numbers::fixed_point::utils; + +/// A struct representing a fixed point number. +#[derive(Serde, Copy, Drop)] +struct FP8x23W { + mag: u64, + sign: bool +} + +// CONSTANTS + +const TWO: u64 = 16777216; // 2 ** 24 +const ONE: u64 = 8388608; // 2 ** 23 +const HALF: u64 = 4194304; // 2 ** 22 +const MAX: u64 = 2147483648; // 2 ** 31 + + +impl FP8x23WImpl of FixedTrait { + fn ZERO() -> FP8x23W { + return FP8x23W { mag: 0, sign: false }; + } + + fn ONE() -> FP8x23W { + return FP8x23W { mag: ONE, sign: false }; + } + + fn MAX() -> FP8x23W { + return FP8x23W { mag: MAX, sign: false }; + } + + fn new(mag: u64, sign: bool) -> FP8x23W { + return FP8x23W { mag: mag, sign: sign }; + } + + fn new_unscaled(mag: u64, sign: bool) -> FP8x23W { + return FP8x23W { mag: mag * ONE, sign: sign }; + } + + fn from_felt(val: felt252) -> FP8x23W { + let mag = integer::u64_try_from_felt252(utils::felt_abs(val)).unwrap(); + return FixedTrait::new(mag, utils::felt_sign(val)); + } + + fn abs(self: FP8x23W) -> FP8x23W { + return core::abs(self); + } + + fn acos(self: FP8x23W) -> FP8x23W { + return trig::acos_fast(self); + } + + fn acos_fast(self: FP8x23W) -> FP8x23W { + return trig::acos_fast(self); + } + + fn acosh(self: FP8x23W) -> FP8x23W { + return hyp::acosh(self); + } + + fn asin(self: FP8x23W) -> FP8x23W { + return trig::asin_fast(self); + } + + fn asin_fast(self: FP8x23W) -> FP8x23W { + return trig::asin_fast(self); + } + + fn asinh(self: FP8x23W) -> FP8x23W { + return hyp::asinh(self); + } + + fn atan(self: FP8x23W) -> FP8x23W { + return trig::atan_fast(self); + } + + fn atan_fast(self: FP8x23W) -> FP8x23W { + return trig::atan_fast(self); + } + + fn atanh(self: FP8x23W) -> FP8x23W { + return hyp::atanh(self); + } + + fn ceil(self: FP8x23W) -> FP8x23W { + return core::ceil(self); + } + + fn cos(self: FP8x23W) -> FP8x23W { + return trig::cos_fast(self); + } + + fn cos_fast(self: FP8x23W) -> FP8x23W { + return trig::cos_fast(self); + } + + fn cosh(self: FP8x23W) -> FP8x23W { + return hyp::cosh(self); + } + + fn floor(self: FP8x23W) -> FP8x23W { + return core::floor(self); + } + + // Calculates the natural exponent of x: e^x + fn exp(self: FP8x23W) -> FP8x23W { + return core::exp(self); + } + + // Calculates the binary exponent of x: 2^x + fn exp2(self: FP8x23W) -> FP8x23W { + return core::exp2(self); + } + + // Calculates the natural logarithm of x: ln(x) + // self must be greater than zero + fn ln(self: FP8x23W) -> FP8x23W { + return core::ln(self); + } + + // Calculates the binary logarithm of x: log2(x) + // self must be greather than zero + fn log2(self: FP8x23W) -> FP8x23W { + return core::log2(self); + } + + // Calculates the base 10 log of x: log10(x) + // self must be greater than zero + fn log10(self: FP8x23W) -> FP8x23W { + return core::log10(self); + } + + // Calclates the value of x^y and checks for overflow before returning + // self is a fixed point value + // b is a fixed point value + fn pow(self: FP8x23W, b: FP8x23W) -> FP8x23W { + return core::pow(self, b); + } + + fn round(self: FP8x23W) -> FP8x23W { + return core::round(self); + } + + fn sin(self: FP8x23W) -> FP8x23W { + return trig::sin_fast(self); + } + + fn sin_fast(self: FP8x23W) -> FP8x23W { + return trig::sin_fast(self); + } + + fn sinh(self: FP8x23W) -> FP8x23W { + return hyp::sinh(self); + } + + // Calculates the square root of a fixed point value + // x must be positive + fn sqrt(self: FP8x23W) -> FP8x23W { + return core::sqrt(self); + } + + fn tan(self: FP8x23W) -> FP8x23W { + return trig::tan_fast(self); + } + + fn tan_fast(self: FP8x23W) -> FP8x23W { + return trig::tan_fast(self); + } + + fn tanh(self: FP8x23W) -> FP8x23W { + return hyp::tanh(self); + } + + fn sign(self: FP8x23W) -> FP8x23W { + return core::sign(self); + } +} + + +impl FP8x23WPrint of PrintTrait { + fn print(self: FP8x23W) { + self.sign.print(); + self.mag.print(); + } +} + +// Into a raw felt without unscaling +impl FP8x23WIntoFelt252 of Into { + fn into(self: FP8x23W) -> felt252 { + let mag_felt = self.mag.into(); + + if self.sign { + return mag_felt * -1; + } else { + return mag_felt * 1; + } + } +} + +impl FP8x23WTryIntoU128 of TryInto { + fn try_into(self: FP8x23W) -> Option { + if self.sign { + return Option::None(()); + } else { + // Unscale the magnitude and round down + return Option::Some((self.mag / ONE).into()); + } + } +} + +impl FP8x23WTryIntoU64 of TryInto { + fn try_into(self: FP8x23W) -> Option { + if self.sign { + return Option::None(()); + } else { + // Unscale the magnitude and round down + return Option::Some((self.mag / ONE).into()); + } + } +} + + +impl FP8x23WTryIntoU16 of TryInto { + fn try_into(self: FP8x23W) -> Option { + if self.sign { + Option::None(()) + } else { + // Unscale the magnitude and round down + return (self.mag / ONE).try_into(); + } + } +} + +impl FP8x23WTryIntoU8 of TryInto { + fn try_into(self: FP8x23W) -> Option { + if self.sign { + Option::None(()) + } else { + // Unscale the magnitude and round down + return (self.mag / ONE).try_into(); + } + } +} + +impl FP8x23WIntoI32 of Into { + fn into(self: FP8x23W) -> i32 { + _i32_into_fp(self) + } +} + +impl FP8x23WTryIntoI8 of TryInto { + fn try_into(self: FP8x23W) -> Option { + _i8_try_from_fp(self) + } +} + +impl FP8x23WPartialEq of PartialEq { + #[inline(always)] + fn eq(lhs: @FP8x23W, rhs: @FP8x23W) -> bool { + return core::eq(lhs, rhs); + } + + #[inline(always)] + fn ne(lhs: @FP8x23W, rhs: @FP8x23W) -> bool { + return core::ne(lhs, rhs); + } +} + +impl FP8x23WAdd of Add { + fn add(lhs: FP8x23W, rhs: FP8x23W) -> FP8x23W { + return core::add(lhs, rhs); + } +} + +impl FP8x23WAddEq of AddEq { + #[inline(always)] + fn add_eq(ref self: FP8x23W, other: FP8x23W) { + self = Add::add(self, other); + } +} + +impl FP8x23WSub of Sub { + fn sub(lhs: FP8x23W, rhs: FP8x23W) -> FP8x23W { + return core::sub(lhs, rhs); + } +} + +impl FP8x23WSubEq of SubEq { + #[inline(always)] + fn sub_eq(ref self: FP8x23W, other: FP8x23W) { + self = Sub::sub(self, other); + } +} + +impl FP8x23WMul of Mul { + fn mul(lhs: FP8x23W, rhs: FP8x23W) -> FP8x23W { + return core::mul(lhs, rhs); + } +} + +impl FP8x23WMulEq of MulEq { + #[inline(always)] + fn mul_eq(ref self: FP8x23W, other: FP8x23W) { + self = Mul::mul(self, other); + } +} + +impl FP8x23WDiv of Div { + fn div(lhs: FP8x23W, rhs: FP8x23W) -> FP8x23W { + return core::div(lhs, rhs); + } +} + +impl FP8x23WDivEq of DivEq { + #[inline(always)] + fn div_eq(ref self: FP8x23W, other: FP8x23W) { + self = Div::div(self, other); + } +} + +impl FP8x23WPartialOrd of PartialOrd { + #[inline(always)] + fn ge(lhs: FP8x23W, rhs: FP8x23W) -> bool { + return core::ge(lhs, rhs); + } + + #[inline(always)] + fn gt(lhs: FP8x23W, rhs: FP8x23W) -> bool { + return core::gt(lhs, rhs); + } + + #[inline(always)] + fn le(lhs: FP8x23W, rhs: FP8x23W) -> bool { + return core::le(lhs, rhs); + } + + #[inline(always)] + fn lt(lhs: FP8x23W, rhs: FP8x23W) -> bool { + return core::lt(lhs, rhs); + } +} + +impl FP8x23WNeg of Neg { + #[inline(always)] + fn neg(a: FP8x23W) -> FP8x23W { + return core::neg(a); + } +} + +impl FP8x23WRem of Rem { + #[inline(always)] + fn rem(lhs: FP8x23W, rhs: FP8x23W) -> FP8x23W { + return core::rem(lhs, rhs); + } +} + +/// INTERNAL + +fn _i32_into_fp(x: FP8x23W) -> i32 { + i32 { mag: (x.mag / ONE).try_into().unwrap(), sign: x.sign } +} + +fn _i8_try_from_fp(x: FP8x23W) -> Option { + let unscaled_mag: Option = (x.mag / ONE).try_into(); + + match unscaled_mag { + Option::Some(val) => Option::Some(i8 { mag: unscaled_mag.unwrap(), sign: x.sign }), + Option::None(_) => Option::None(()) + } +} diff --git a/src/numbers/fixed_point/implementations/fp8x23wide/helpers.cairo b/src/numbers/fixed_point/implementations/fp8x23wide/helpers.cairo new file mode 100644 index 000000000..a627803be --- /dev/null +++ b/src/numbers/fixed_point/implementations/fp8x23wide/helpers.cairo @@ -0,0 +1,39 @@ +use debug::PrintTrait; +use traits::Into; + +use orion::numbers::fixed_point::implementations::fp8x23wide::core::{ + HALF, ONE, TWO, FP8x23W, FP8x23WSub, FP8x23WDiv, FixedTrait, FP8x23WPrint +}; + +const DEFAULT_PRECISION: u64 = 8; // 1e-6 + +// To use `DEFAULT_PRECISION`, final arg is: `Option::None(())`. +// To use `custom_precision` of 430_u64: `Option::Some(430_u64)`. +fn assert_precise(result: FP8x23W, expected: felt252, msg: felt252, custom_precision: Option) { + let precision = match custom_precision { + Option::Some(val) => val, + Option::None(_) => DEFAULT_PRECISION, + }; + + let diff = (result - FixedTrait::from_felt(expected)).mag; + + if (diff > precision) { + result.print(); + assert(diff <= precision, msg); + } +} + +fn assert_relative(result: FP8x23W, expected: felt252, msg: felt252, custom_precision: Option) { + let precision = match custom_precision { + Option::Some(val) => val, + Option::None(_) => DEFAULT_PRECISION, + }; + + let diff = result - FixedTrait::from_felt(expected); + let rel_diff = (diff / result).mag; + + if (rel_diff > precision) { + result.print(); + assert(rel_diff <= precision, msg); + } +} diff --git a/src/numbers/fixed_point/implementations/fp8x23wide/math.cairo b/src/numbers/fixed_point/implementations/fp8x23wide/math.cairo new file mode 100644 index 000000000..970c65f30 --- /dev/null +++ b/src/numbers/fixed_point/implementations/fp8x23wide/math.cairo @@ -0,0 +1,5 @@ +mod core; +mod comp; +mod lut; +mod trig; +mod hyp; diff --git a/src/numbers/fixed_point/implementations/fp8x23wide/math/comp.cairo b/src/numbers/fixed_point/implementations/fp8x23wide/math/comp.cairo new file mode 100644 index 000000000..95b329109 --- /dev/null +++ b/src/numbers/fixed_point/implementations/fp8x23wide/math/comp.cairo @@ -0,0 +1,76 @@ +use orion::numbers::fixed_point::implementations::fp8x23wide::core::{ + FP8x23W, FixedTrait, FP8x23WPartialOrd, FP8x23WPartialEq +}; + +fn max(a: FP8x23W, b: FP8x23W) -> FP8x23W { + if (a >= b) { + return a; + } else { + return b; + } +} + +fn min(a: FP8x23W, b: FP8x23W) -> FP8x23W { + if (a <= b) { + return a; + } else { + return b; + } +} + +fn xor(a: FP8x23W, b: FP8x23W) -> bool { + if (a == FixedTrait::new(0, false) || b == FixedTrait::new(0, false)) && (a != b) { + return true; + } else { + return false; + } +} + +fn or(a: FP8x23W, b: FP8x23W) -> bool { + let zero = FixedTrait::new(0, false); + if a == zero && b == zero { + return false; + } else { + return true; + } +} + +// Tests -------------------------------------------------------------------------------------------------------------- + +#[test] +fn test_max() { + let a = FixedTrait::new_unscaled(1, false); + let b = FixedTrait::new_unscaled(0, false); + let c = FixedTrait::new_unscaled(1, true); + + assert(max(a, a) == a, 'max(a, a)'); + assert(max(a, b) == a, 'max(a, b)'); + assert(max(a, c) == a, 'max(a, c)'); + + assert(max(b, a) == a, 'max(b, a)'); + assert(max(b, b) == b, 'max(b, b)'); + assert(max(b, c) == b, 'max(b, c)'); + + assert(max(c, a) == a, 'max(c, a)'); + assert(max(c, b) == b, 'max(c, b)'); + assert(max(c, c) == c, 'max(c, c)'); +} + +#[test] +fn test_min() { + let a = FixedTrait::new_unscaled(1, false); + let b = FixedTrait::new_unscaled(0, false); + let c = FixedTrait::new_unscaled(1, true); + + assert(min(a, a) == a, 'min(a, a)'); + assert(min(a, b) == b, 'min(a, b)'); + assert(min(a, c) == c, 'min(a, c)'); + + assert(min(b, a) == b, 'min(b, a)'); + assert(min(b, b) == b, 'min(b, b)'); + assert(min(b, c) == c, 'min(b, c)'); + + assert(min(c, a) == c, 'min(c, a)'); + assert(min(c, b) == c, 'min(c, b)'); + assert(min(c, c) == c, 'min(c, c)'); +} diff --git a/src/numbers/fixed_point/implementations/fp8x23wide/math/core.cairo b/src/numbers/fixed_point/implementations/fp8x23wide/math/core.cairo new file mode 100644 index 000000000..129ff02c8 --- /dev/null +++ b/src/numbers/fixed_point/implementations/fp8x23wide/math/core.cairo @@ -0,0 +1,660 @@ +use core::debug::PrintTrait; +use option::OptionTrait; +use result::{ResultTrait, ResultTraitImpl}; +use traits::{Into, TryInto}; +use integer::{u64_safe_divmod, u64_as_non_zero, u64_wide_mul}; + +use orion::numbers::fixed_point::implementations::fp8x23wide::core::{ + HALF, ONE, MAX, FP8x23W, FP8x23WAdd, FP8x23WImpl, FP8x23WAddEq, FP8x23WSub, FP8x23WMul, + FP8x23WMulEq, FP8x23WTryIntoU128, FP8x23WPartialEq, FP8x23WPartialOrd, FP8x23WSubEq, FP8x23WNeg, + FP8x23WDiv, FP8x23WIntoFelt252, FixedTrait +}; +use orion::numbers::fixed_point::implementations::fp8x23wide::math::lut; + +// PUBLIC + +fn abs(a: FP8x23W) -> FP8x23W { + return FixedTrait::new(a.mag, false); +} + +fn add(a: FP8x23W, b: FP8x23W) -> FP8x23W { + if a.sign == b.sign { + return FixedTrait::new(a.mag + b.mag, a.sign); + } + + if a.mag == b.mag { + return FixedTrait::ZERO(); + } + + if (a.mag > b.mag) { + return FixedTrait::new(a.mag - b.mag, a.sign); + } else { + return FixedTrait::new(b.mag - a.mag, b.sign); + } +} + +fn ceil(a: FP8x23W) -> FP8x23W { + let (div, rem) = u64_safe_divmod(a.mag, u64_as_non_zero(ONE)); + + if rem == 0 { + return a; + } else if !a.sign { + return FixedTrait::new_unscaled(div + 1, false); + } else if div == 0 { + return FixedTrait::new_unscaled(0, false); + } else { + return FixedTrait::new_unscaled(div, true); + } +} + +fn div(a: FP8x23W, b: FP8x23W) -> FP8x23W { + let a_u64 = integer::u64_wide_mul(a.mag, ONE); + let res_u64 = a_u64 / b.mag.into(); + + // Re-apply sign + return FixedTrait::new(res_u64.try_into().unwrap(), a.sign ^ b.sign); +} + +fn eq(a: @FP8x23W, b: @FP8x23W) -> bool { + return (*a.mag == *b.mag) && (*a.sign == *b.sign); +} + +// Calculates the natural exponent of x: e^x +fn exp(a: FP8x23W) -> FP8x23W { + return exp2(FixedTrait::new(12102203, false) * a); // log2(e) * 2^23 ≈ 12102203 +} + +// Calculates the binary exponent of x: 2^x +fn exp2(a: FP8x23W) -> FP8x23W { + if (a.mag == 0) { + return FixedTrait::ONE(); + } + + let (int_part, frac_part) = integer::u64_safe_divmod(a.mag, u64_as_non_zero(ONE)); + let int_res = FixedTrait::new_unscaled(lut::exp2(int_part), false); + let mut res_u = int_res; + + if frac_part != 0 { + let frac = FixedTrait::new(frac_part, false); + let r8 = FixedTrait::new(19, false) * frac; + let r7 = (r8 + FixedTrait::new(105, false)) * frac; + let r6 = (r7 + FixedTrait::new(1324, false)) * frac; + let r5 = (r6 + FixedTrait::new(11159, false)) * frac; + let r4 = (r5 + FixedTrait::new(80695, false)) * frac; + let r3 = (r4 + FixedTrait::new(465599, false)) * frac; + let r2 = (r3 + FixedTrait::new(2015166, false)) * frac; + let r1 = (r2 + FixedTrait::new(5814540, false)) * frac; + res_u = res_u * (r1 + FixedTrait::ONE()); + } + + if (a.sign == true) { + return FixedTrait::ONE() / res_u; + } else { + return res_u; + } +} + +fn exp2_int(exp: u64) -> FP8x23W { + return FixedTrait::new_unscaled(lut::exp2(exp), false); +} + +fn floor(a: FP8x23W) -> FP8x23W { + let (div, rem) = integer::u64_safe_divmod(a.mag, u64_as_non_zero(ONE)); + + if rem == 0 { + return a; + } else if !a.sign { + return FixedTrait::new_unscaled(div, false); + } else { + return FixedTrait::new_unscaled(div + 1, true); + } +} + +fn ge(a: FP8x23W, b: FP8x23W) -> bool { + if a.sign != b.sign { + return !a.sign; + } else { + return (a.mag == b.mag) || ((a.mag > b.mag) ^ a.sign); + } +} + +fn gt(a: FP8x23W, b: FP8x23W) -> bool { + if a.sign != b.sign { + return !a.sign; + } else { + return (a.mag != b.mag) && ((a.mag > b.mag) ^ a.sign); + } +} + +fn le(a: FP8x23W, b: FP8x23W) -> bool { + if a.sign != b.sign { + return a.sign; + } else { + return (a.mag == b.mag) || ((a.mag < b.mag) ^ a.sign); + } +} + +// Calculates the natural logarithm of x: ln(x) +// self must be greater than zero +fn ln(a: FP8x23W) -> FP8x23W { + return FixedTrait::new(5814540, false) * log2(a); // ln(2) = 0.693... +} + +// Calculates the binary logarithm of x: log2(x) +// self must be greather than zero +fn log2(a: FP8x23W) -> FP8x23W { + assert(a.sign == false, 'must be positive'); + + if (a.mag == ONE) { + return FixedTrait::ZERO(); + } else if (a.mag < ONE) { + // Compute true inverse binary log if 0 < x < 1 + let div = FixedTrait::ONE() / a; + return -log2(div); + } + + let whole = a.mag / ONE; + let (msb, div) = lut::msb(whole); + + if a.mag == div * ONE { + return FixedTrait::new_unscaled(msb, false); + } else { + let norm = a / FixedTrait::new_unscaled(div, false); + let r8 = FixedTrait::new(76243, true) * norm; + let r7 = (r8 + FixedTrait::new(1038893, false)) * norm; + let r6 = (r7 + FixedTrait::new(6277679, true)) * norm; + let r5 = (r6 + FixedTrait::new(22135645, false)) * norm; + let r4 = (r5 + FixedTrait::new(50444339, true)) * norm; + let r3 = (r4 + FixedTrait::new(77896489, false)) * norm; + let r2 = (r3 + FixedTrait::new(83945943, true)) * norm; + let r1 = (r2 + FixedTrait::new(68407458, false)) * norm; + return r1 + FixedTrait::new(28734280, true) + FixedTrait::new_unscaled(msb, false); + } +} + +// Calculates the base 10 log of x: log10(x) +// self must be greater than zero +fn log10(a: FP8x23W) -> FP8x23W { + return FixedTrait::new(2525223, false) * log2(a); // log10(2) = 0.301... +} + +fn lt(a: FP8x23W, b: FP8x23W) -> bool { + if a.sign != b.sign { + return a.sign; + } else { + return (a.mag != b.mag) && ((a.mag < b.mag) ^ a.sign); + } +} + +fn mul(a: FP8x23W, b: FP8x23W) -> FP8x23W { + let prod_u128 = integer::u64_wide_mul(a.mag, b.mag); + + // Re-apply sign + return FixedTrait::new((prod_u128 / ONE.into()).try_into().unwrap(), a.sign ^ b.sign); +} + +fn ne(a: @FP8x23W, b: @FP8x23W) -> bool { + return (*a.mag != *b.mag) || (*a.sign != *b.sign); +} + +fn neg(a: FP8x23W) -> FP8x23W { + if a.mag == 0 { + return a; + } else if !a.sign { + return FixedTrait::new(a.mag, !a.sign); + } else { + return FixedTrait::new(a.mag, false); + } +} + +// Calclates the value of x^y and checks for overflow before returning +// self is a FP8x23W point value +// b is a FP8x23W point value +fn pow(a: FP8x23W, b: FP8x23W) -> FP8x23W { + let (div, rem) = integer::u64_safe_divmod(b.mag, u64_as_non_zero(ONE)); + + // use the more performant integer pow when y is an int + if (rem == 0) { + return pow_int(a, b.mag / ONE, b.sign); + } + + // x^y = exp(y*ln(x)) for x > 0 will error for x < 0 + return exp(b * ln(a)); +} + +// Calclates the value of a^b and checks for overflow before returning +fn pow_int(a: FP8x23W, b: u64, sign: bool) -> FP8x23W { + let mut x = a; + let mut n = b; + + if sign == true { + x = FixedTrait::ONE() / x; + } + + if n == 0 { + return FixedTrait::ONE(); + } + + let mut y = FixedTrait::ONE(); + let two = integer::u64_as_non_zero(2); + + loop { + if n <= 1 { + break; + } + + let (div, rem) = integer::u64_safe_divmod(n, two); + + if rem == 1 { + y = x * y; + } + + x = x * x; + n = div; + }; + + return x * y; +} + +fn rem(a: FP8x23W, b: FP8x23W) -> FP8x23W { + return a - floor(a / b) * b; +} + +fn round(a: FP8x23W) -> FP8x23W { + let (div, rem) = integer::u64_safe_divmod(a.mag, u64_as_non_zero(ONE)); + + if (HALF <= rem) { + return FixedTrait::new_unscaled(div + 1, a.sign); + } else { + return FixedTrait::new_unscaled(div, a.sign); + } +} + +// Calculates the square root of a FP8x23W point value +// x must be positive +fn sqrt(a: FP8x23W) -> FP8x23W { + assert(a.sign == false, 'must be positive'); + + let root = integer::u64_sqrt(a.mag.into() * ONE.into()); + return FixedTrait::new(root.into(), false); +} + +fn sub(a: FP8x23W, b: FP8x23W) -> FP8x23W { + return add(a, -b); +} + +fn sign(a: FP8x23W) -> FP8x23W { + if a.mag == 0 { + FixedTrait::new(0, false) + } else { + FixedTrait::new(ONE, a.sign) + } +} + +// Tests -------------------------------------------------------------------------------------------------------------- + +use orion::numbers::fixed_point::implementations::fp8x23wide::helpers::{ + assert_precise, assert_relative +}; +use orion::numbers::fixed_point::implementations::fp8x23wide::math::trig::{PI, HALF_PI}; + +#[test] +fn test_into() { + let a = FixedTrait::::new_unscaled(5, false); + assert(a.mag == 5 * ONE, 'invalid result'); +} + +#[test] +fn test_try_into_u128() { + // Positive unscaled + let a = FixedTrait::::new_unscaled(5, false); + assert(a.try_into().unwrap() == 5_u128, 'invalid result'); + + // Positive scaled + let b = FixedTrait::::new(5 * ONE, false); + assert(b.try_into().unwrap() == 5_u128, 'invalid result'); + + // Zero + let d = FixedTrait::::new_unscaled(0, false); + assert(d.try_into().unwrap() == 0_u128, 'invalid result'); +} + +#[test] +#[should_panic] +fn test_negative_try_into_u128() { + let a = FixedTrait::::new_unscaled(1, true); + let a: u128 = a.try_into().unwrap(); +} + +#[test] +#[available_gas(1000000)] +fn test_acos() { + let a = FixedTrait::::ONE(); + assert(a.acos().into() == 0, 'invalid one'); +} + +#[test] +#[available_gas(1000000)] +fn test_asin() { + let a = FixedTrait::ONE(); + assert_precise(a.asin(), HALF_PI.into(), 'invalid one', Option::None(())); // PI / 2 +} + +#[test] +#[available_gas(2000000)] +fn test_atan() { + let a = FixedTrait::new(2 * ONE, false); + assert_relative(a.atan(), 9287469, 'invalid two', Option::None(())); +} + +#[test] +fn test_ceil() { + let a = FixedTrait::new(24326963, false); // 2.9 + assert(ceil(a).mag == 3 * ONE, 'invalid pos decimal'); +} + +#[test] +fn test_floor() { + let a = FixedTrait::new(24326963, false); // 2.9 + assert(floor(a).mag == 2 * ONE, 'invalid pos decimal'); +} + +#[test] +fn test_round() { + let a = FixedTrait::new(24326963, false); // 2.9 + assert(round(a).mag == 3 * ONE, 'invalid pos decimal'); +} + +#[test] +#[should_panic] +fn test_sqrt_fail() { + let a = FixedTrait::new_unscaled(25, true); + sqrt(a); +} + +#[test] +fn test_sqrt() { + let mut a = FixedTrait::new_unscaled(0, false); + assert(sqrt(a).mag == 0, 'invalid zero root'); + a = FixedTrait::new_unscaled(25, false); + assert(sqrt(a).mag == 5 * ONE, 'invalid pos root'); +} + + +#[test] +#[available_gas(100000)] +fn test_msb() { + let a = FixedTrait::::new_unscaled(100, false); + let (msb, div) = lut::msb(a.mag / ONE); + assert(msb == 6, 'invalid msb'); + assert(div == 64, 'invalid msb ceil'); +} + +#[test] +#[available_gas(600000)] +fn test_pow() { + let a = FixedTrait::new_unscaled(3, false); + let b = FixedTrait::new_unscaled(4, false); + assert(pow(a, b).mag == 81 * ONE, 'invalid pos base power'); +} + +#[test] +#[available_gas(900000)] +fn test_pow_frac() { + let a = FixedTrait::new_unscaled(3, false); + let b = FixedTrait::new(4194304, false); // 0.5 + assert_relative( + pow(a, b), 14529495, 'invalid pos base power', Option::None(()) + ); // 1.7320508075688772 +} + +#[test] +#[available_gas(1000000)] +fn test_exp() { + let a = FixedTrait::new_unscaled(2, false); + assert_relative(exp(a), 61983895, 'invalid exp of 2', Option::None(())); // 7.389056098793725 +} + +#[test] +#[available_gas(400000)] +fn test_exp2() { + let a = FixedTrait::new_unscaled(5, false); + assert(exp2(a).mag == 268435456, 'invalid exp2 of 2'); +} + +#[test] +#[available_gas(20000)] +fn test_exp2_int() { + assert(exp2_int(5).into() == 268435456, 'invalid exp2 of 2'); +} + +#[test] +#[available_gas(1000000)] +fn test_ln() { + let mut a = FixedTrait::new_unscaled(1, false); + assert(ln(a).mag == 0, 'invalid ln of 1'); + + a = FixedTrait::new(22802601, false); + assert_relative(ln(a), ONE.into(), 'invalid ln of 2.7...', Option::None(())); +} + +#[test] +#[available_gas(1000000)] +fn test_log2() { + let mut a = FixedTrait::new_unscaled(32, false); + assert(log2(a) == FixedTrait::new_unscaled(5, false), 'invalid log2 32'); + + a = FixedTrait::new_unscaled(10, false); + assert_relative(log2(a), 27866353, 'invalid log2 10', Option::None(())); // 3.321928094887362 +} + +#[test] +#[available_gas(1000000)] +fn test_log10() { + let a = FixedTrait::new_unscaled(100, false); + assert_relative(log10(a), 2 * ONE.into(), 'invalid log10', Option::None(())); +} + +#[test] +fn test_eq() { + let a = FixedTrait::new_unscaled(42, false); + let b = FixedTrait::new_unscaled(42, false); + let c = eq(@a, @b); + assert(c == true, 'invalid result'); +} + +#[test] +fn test_ne() { + let a = FixedTrait::new_unscaled(42, false); + let b = FixedTrait::new_unscaled(42, false); + let c = ne(@a, @b); + assert(c == false, 'invalid result'); +} + +#[test] +fn test_add() { + let a = FixedTrait::new_unscaled(1, false); + let b = FixedTrait::new_unscaled(2, false); + assert(add(a, b) == FixedTrait::new_unscaled(3, false), 'invalid result'); +} + +#[test] +fn test_add_eq() { + let mut a = FixedTrait::new_unscaled(1, false); + let b = FixedTrait::new_unscaled(2, false); + a += b; + assert(a == FixedTrait::::new_unscaled(3, false), 'invalid result'); +} + +#[test] +fn test_sub() { + let a = FixedTrait::new_unscaled(5, false); + let b = FixedTrait::new_unscaled(2, false); + let c = a - b; + assert(c == FixedTrait::::new_unscaled(3, false), 'false result invalid'); +} + +#[test] +fn test_sub_eq() { + let mut a = FixedTrait::new_unscaled(5, false); + let b = FixedTrait::new_unscaled(2, false); + a -= b; + assert(a == FixedTrait::::new_unscaled(3, false), 'invalid result'); +} + +#[test] +#[available_gas(100000)] +fn test_mul_pos() { + let a = FP8x23W { mag: 24326963, sign: false }; + let b = FP8x23W { mag: 24326963, sign: false }; + let c = a * b; + assert(c.mag == 70548192, 'invalid result'); +} + +#[test] +fn test_mul_neg() { + let a = FixedTrait::new_unscaled(5, false); + let b = FixedTrait::new_unscaled(2, true); + let c = a * b; + assert(c == FixedTrait::::new_unscaled(10, true), 'invalid result'); +} + +#[test] +fn test_mul_eq() { + let mut a = FixedTrait::new_unscaled(5, false); + let b = FixedTrait::new_unscaled(2, true); + a *= b; + assert(a == FixedTrait::::new_unscaled(10, true), 'invalid result'); +} + +#[test] +fn test_div() { + let a = FixedTrait::new_unscaled(10, false); + let b = FixedTrait::::new(24326963, false); // 2.9 + let c = a / b; + assert(c.mag == 28926234, 'invalid pos decimal'); // 3.4482758620689653 +} + +#[test] +fn test_le() { + let a = FixedTrait::new_unscaled(1, false); + let b = FixedTrait::new_unscaled(0, false); + let c = FixedTrait::::new_unscaled(1, true); + + assert(a <= a, 'a <= a'); + assert(a <= b == false, 'a <= b'); + assert(a <= c == false, 'a <= c'); + + assert(b <= a, 'b <= a'); + assert(b <= b, 'b <= b'); + assert(b <= c == false, 'b <= c'); + + assert(c <= a, 'c <= a'); + assert(c <= b, 'c <= b'); + assert(c <= c, 'c <= c'); +} + +#[test] +fn test_lt() { + let a = FixedTrait::new_unscaled(1, false); + let b = FixedTrait::new_unscaled(0, false); + let c = FixedTrait::::new_unscaled(1, true); + + assert(a < a == false, 'a < a'); + assert(a < b == false, 'a < b'); + assert(a < c == false, 'a < c'); + + assert(b < a, 'b < a'); + assert(b < b == false, 'b < b'); + assert(b < c == false, 'b < c'); + + assert(c < a, 'c < a'); + assert(c < b, 'c < b'); + assert(c < c == false, 'c < c'); +} + +#[test] +fn test_ge() { + let a = FixedTrait::new_unscaled(1, false); + let b = FixedTrait::new_unscaled(0, false); + let c = FixedTrait::::new_unscaled(1, true); + + assert(a >= a, 'a >= a'); + assert(a >= b, 'a >= b'); + assert(a >= c, 'a >= c'); + + assert(b >= a == false, 'b >= a'); + assert(b >= b, 'b >= b'); + assert(b >= c, 'b >= c'); + + assert(c >= a == false, 'c >= a'); + assert(c >= b == false, 'c >= b'); + assert(c >= c, 'c >= c'); +} + +#[test] +fn test_gt() { + let a = FixedTrait::new_unscaled(1, false); + let b = FixedTrait::new_unscaled(0, false); + let c = FixedTrait::::new_unscaled(1, true); + + assert(a > a == false, 'a > a'); + assert(a > b, 'a > b'); + assert(a > c, 'a > c'); + + assert(b > a == false, 'b > a'); + assert(b > b == false, 'b > b'); + assert(b > c, 'b > c'); + + assert(c > a == false, 'c > a'); + assert(c > b == false, 'c > b'); + assert(c > c == false, 'c > c'); +} + +#[test] +#[available_gas(1000000)] +fn test_cos() { + let a = FixedTrait::::new(HALF_PI, false); + assert(a.cos().into() == 0, 'invalid half pi'); +} + +#[test] +#[available_gas(1000000)] +fn test_sin() { + let a = FixedTrait::new(HALF_PI, false); + assert_precise(a.sin(), ONE.into(), 'invalid half pi', Option::None(())); +} + +#[test] +#[available_gas(2000000)] +fn test_tan() { + let a = FixedTrait::::new(HALF_PI / 2, false); + assert(a.tan().mag == 8388608, 'invalid quarter pi'); +} + +#[test] +#[available_gas(2000000)] +fn test_sign() { + let a = FixedTrait::::new(0, false); + assert(a.sign().mag == 0 && !a.sign().sign, 'invalid sign (0, true)'); + + let a = FixedTrait::::new(HALF, true); + assert(a.sign().mag == ONE && a.sign().sign, 'invalid sign (HALF, true)'); + + let a = FixedTrait::::new(HALF, false); + assert(a.sign().mag == ONE && !a.sign().sign, 'invalid sign (HALF, false)'); + + let a = FixedTrait::::new(ONE, true); + assert(a.sign().mag == ONE && a.sign().sign, 'invalid sign (ONE, true)'); + + let a = FixedTrait::::new(ONE, false); + assert(a.sign().mag == ONE && !a.sign().sign, 'invalid sign (ONE, false)'); +} + +#[test] +#[should_panic] +#[available_gas(2000000)] +fn test_sign_fail() { + let a = FixedTrait::::new(HALF, true); + assert(a.sign().mag != ONE && !a.sign().sign, 'invalid sign (HALF, true)'); +} diff --git a/src/numbers/fixed_point/implementations/fp8x23wide/math/hyp.cairo b/src/numbers/fixed_point/implementations/fp8x23wide/math/hyp.cairo new file mode 100644 index 000000000..ed9b66391 --- /dev/null +++ b/src/numbers/fixed_point/implementations/fp8x23wide/math/hyp.cairo @@ -0,0 +1,159 @@ +use core::debug::PrintTrait; +use orion::numbers::fixed_point::implementations::fp8x23wide::core::{ + HALF, ONE, TWO, FP8x23W, FP8x23WImpl, FP8x23WAdd, FP8x23WAddEq, FP8x23WSub, FP8x23WMul, FP8x23WMulEq, + FP8x23WTryIntoU128, FP8x23WPartialEq, FP8x23WPartialOrd, FP8x23WSubEq, FP8x23WNeg, FP8x23WDiv, + FP8x23WIntoFelt252, FixedTrait +}; + +// Calculates hyperbolic cosine of a (fixed point) +fn cosh(a: FP8x23W) -> FP8x23W { + let ea = a.exp(); + return (ea + (FixedTrait::ONE() / ea)) / FixedTrait::new(TWO, false); +} + +// Calculates hyperbolic sine of a (fixed point) +fn sinh(a: FP8x23W) -> FP8x23W { + let ea = a.exp(); + return (ea - (FixedTrait::ONE() / ea)) / FixedTrait::new(TWO, false); +} + +// Calculates hyperbolic tangent of a (fixed point) +fn tanh(a: FP8x23W) -> FP8x23W { + let ea = a.exp(); + let ea_i = FixedTrait::ONE() / ea; + return (ea - ea_i) / (ea + ea_i); +} + +// Calculates inverse hyperbolic cosine of a (fixed point) +fn acosh(a: FP8x23W) -> FP8x23W { + let root = (a * a - FixedTrait::ONE()).sqrt(); + return (a + root).ln(); +} + +// Calculates inverse hyperbolic sine of a (fixed point) +fn asinh(a: FP8x23W) -> FP8x23W { + let root = (a * a + FixedTrait::ONE()).sqrt(); + return (a + root).ln(); +} + +// Calculates inverse hyperbolic tangent of a (fixed point) +fn atanh(a: FP8x23W) -> FP8x23W { + let one = FixedTrait::ONE(); + let ln_arg = (one + a) / (one - a); + return ln_arg.ln() / FixedTrait::new(TWO, false); +} + +// Tests -------------------------------------------------------------------------------------------------------------- + +use option::OptionTrait; +use traits::Into; + +use orion::numbers::fixed_point::implementations::fp8x23wide::helpers::assert_precise; + +#[test] +#[available_gas(10000000)] +fn test_cosh() { + let a = FixedTrait::new(TWO, false); + assert_precise(cosh(a), 31559585, 'invalid two', Option::None(())); // 3.762195691016423 + + let a = FixedTrait::ONE(); + assert_precise(cosh(a), 12944299, 'invalid one', Option::None(())); // 1.5430806347841253 + + let a = FixedTrait::ZERO(); + assert_precise(cosh(a), ONE.into(), 'invalid zero', Option::None(())); + + let a = FixedTrait::ONE(); + assert_precise(cosh(a), 12944299, 'invalid neg one', Option::None(())); // 1.5430806347841253 + + let a = FixedTrait::new(TWO, true); + assert_precise(cosh(a), 31559602, 'invalid neg two', Option::None(())); // 3.762195691016423 +} + +#[test] +#[available_gas(10000000)] +fn test_sinh() { + let a = FixedTrait::new(TWO, false); + assert_precise(sinh(a), 30424310, 'invalid two', Option::None(())); // 3.6268604077773023 + + let a = FixedTrait::ONE(); + assert_precise(sinh(a), 9858302, 'invalid one', Option::None(())); // 1.1752011936029418 + + let a = FixedTrait::ZERO(); + assert(sinh(a).into() == 0, 'invalid zero'); + + let a = FixedTrait::new(ONE, true); + assert_precise(sinh(a), -9858302, 'invalid neg one', Option::None(())); // -1.1752011936029418 + + let a = FixedTrait::new(TWO, true); + assert_precise(sinh(a), -30424328, 'invalid neg two', Option::None(())); // -3.6268604077773023 +} + +#[test] +#[available_gas(10000000)] +fn test_tanh() { + let a = FixedTrait::new(TWO, false); + assert_precise(tanh(a), 8086849, 'invalid two', Option::None(())); // 0.9640275800745076 + + let a = FixedTrait::ONE(); + assert_precise(tanh(a), 6388715, 'invalid one', Option::None(())); // 0.7615941559446443 + + let a = FixedTrait::ZERO(); + assert(tanh(a).into() == 0, 'invalid zero'); + + let a = FixedTrait::new(ONE, true); + assert_precise(tanh(a), -6388715, 'invalid neg one', Option::None(())); // -0.7615941559446443 + + let a = FixedTrait::new(TWO, true); + assert_precise(tanh(a), -8086849, 'invalid neg two', Option::None(())); // 0.9640275800745076 +} + +#[test] +#[available_gas(10000000)] +fn test_acosh() { + let a = FixedTrait::new(31559585, false); // 3.762195691016423 + assert_precise(acosh(a), 16777257, 'invalid two', Option::None(())); + + let a = FixedTrait::new(12944299, false); // 1.5430806347841253 + assert_precise(acosh(a), ONE.into(), 'invalid one', Option::None(())); + + let a = FixedTrait::ONE(); // 1 + assert(acosh(a).into() == 0, 'invalid zero'); +} + +#[test] +#[available_gas(10000000)] +fn test_asinh() { + let a = FixedTrait::new(30424310, false); // 3.6268604077773023 + assert_precise(asinh(a), 16777257, 'invalid two', Option::None(())); + + let a = FixedTrait::new(9858302, false); // 1.1752011936029418 + assert_precise(asinh(a), ONE.into(), 'invalid one', Option::None(())); + + let a = FixedTrait::ZERO(); + assert(asinh(a).into() == 0, 'invalid zero'); + + let a = FixedTrait::new(9858302, true); // -1.1752011936029418 + assert_precise(asinh(a), -ONE.into(), 'invalid neg one', Option::None(())); + + let a = FixedTrait::new(30424310, true); // -3.6268604077773023 + assert_precise(asinh(a), -16777238, 'invalid neg two', Option::None(())); +} + +#[test] +#[available_gas(10000000)] +fn test_atanh() { + let a = FixedTrait::new(7549747, false); // 0.9 + assert_precise(atanh(a), 12349872, 'invalid 0.9', Option::None(())); // 1.4722194895832204 + + let a = FixedTrait::new(HALF, false); // 0.5 + assert_precise(atanh(a), 4607914, 'invalid half', Option::None(())); // 0.5493061443340548 + + let a = FixedTrait::ZERO(); + assert(atanh(a).into() == 0, 'invalid zero'); + + let a = FixedTrait::new(HALF, true); // 0.5 + assert_precise(atanh(a), -4607914, 'invalid neg half', Option::None(())); // 0.5493061443340548 + + let a = FixedTrait::new(7549747, true); // 0.9 + assert_precise(atanh(a), -12349872, 'invalid -0.9', Option::None(())); // 1.4722194895832204 +} diff --git a/src/numbers/fixed_point/implementations/fp8x23wide/math/lut.cairo b/src/numbers/fixed_point/implementations/fp8x23wide/math/lut.cairo new file mode 100644 index 000000000..157499b5b --- /dev/null +++ b/src/numbers/fixed_point/implementations/fp8x23wide/math/lut.cairo @@ -0,0 +1,1229 @@ +// Calculates the most significant bit +fn msb(whole: u64) -> (u64, u64) { + if whole < 256 { + if whole < 2 { + return (0, 1); + } + if whole < 4 { + return (1, 2); + } + if whole < 8 { + return (2, 4); + } + if whole < 16 { + return (3, 8); + } + if whole < 32 { + return (4, 16); + } + if whole < 64 { + return (5, 32); + } + if whole < 128 { + return (6, 64); + } + if whole < 256 { + return (7, 128); + } + } + + return (8, 256); +} + +fn exp2(exp: u64) -> u64 { + if exp <= 16 { + if exp == 0 { + return 1; + } + if exp == 1 { + return 2; + } + if exp == 2 { + return 4; + } + if exp == 3 { + return 8; + } + if exp == 4 { + return 16; + } + if exp == 5 { + return 32; + } + if exp == 6 { + return 64; + } + if exp == 7 { + return 128; + } + if exp == 8 { + return 256; + } + if exp == 9 { + return 512; + } + if exp == 10 { + return 1024; + } + if exp == 11 { + return 2048; + } + if exp == 12 { + return 4096; + } + if exp == 13 { + return 8192; + } + if exp == 14 { + return 16384; + } + if exp == 15 { + return 32768; + } + if exp == 16 { + return 65536; + } + } else if exp <= 32 { + if exp == 17 { + return 131072; + } + if exp == 18 { + return 262144; + } + if exp == 19 { + return 524288; + } + if exp == 20 { + return 1048576; + } + if exp == 21 { + return 2097152; + } + if exp == 22 { + return 4194304; + } + } + + return 8388608; +} + +fn sin(a: u64) -> (u64, u64, u64) { + let slot = a / 51472; + + if slot < 128 { + if slot < 64 { + if slot < 32 { + if slot < 16 { + if slot == 0 { + return (0, 0, 51472); + } + if slot == 1 { + return (51472, 51472, 102941); + } + if slot == 2 { + return (102944, 102941, 154407); + } + if slot == 3 { + return (154416, 154407, 205867); + } + if slot == 4 { + return (205887, 205867, 257319); + } + if slot == 5 { + return (257359, 257319, 308761); + } + if slot == 6 { + return (308831, 308761, 360192); + } + if slot == 7 { + return (360303, 360192, 411609); + } + if slot == 8 { + return (411775, 411609, 463011); + } + if slot == 9 { + return (463247, 463011, 514396); + } + if slot == 10 { + return (514723, 514396, 565761); + } + if slot == 11 { + return (566190, 565761, 617104); + } + if slot == 12 { + return (617662, 617104, 668425); + } + if slot == 13 { + return (669134, 668425, 719720); + } + if slot == 14 { + return (720606, 719720, 770988); + } + if slot == 15 { + return (772078, 770988, 822227); + } + } else { + if slot == 16 { + return (823550, 822227, 873436); + } + if slot == 17 { + return (875022, 873436, 924611); + } + if slot == 18 { + return (926493, 924611, 975751); + } + if slot == 19 { + return (977965, 975751, 1026855); + } + if slot == 20 { + return (1029437, 1026855, 1077920); + } + if slot == 21 { + return (1080909, 1077920, 1128945); + } + if slot == 22 { + return (1132381, 1128945, 1179927); + } + if slot == 23 { + return (1183853, 1179927, 1230864); + } + if slot == 24 { + return (1235324, 1230864, 1281756); + } + if slot == 25 { + return (1286796, 1281756, 1332599); + } + if slot == 26 { + return (1338268, 1332599, 1383392); + } + if slot == 27 { + return (1389740, 1383392, 1434132); + } + if slot == 28 { + return (1441212, 1434132, 1484819); + } + if slot == 29 { + return (1492684, 1484819, 1535450); + } + if slot == 30 { + return (1544156, 1535450, 1586023); + } + if slot == 31 { + return (1595627, 1586023, 1636536); + } + } + } else { + if slot < 48 { + if slot == 32 { + return (1647099, 1636536, 1686988); + } + if slot == 33 { + return (1698571, 1686988, 1737376); + } + if slot == 34 { + return (1750043, 1737376, 1787699); + } + if slot == 35 { + return (1801515, 1787699, 1837954); + } + if slot == 36 { + return (1852987, 1837954, 1888141); + } + if slot == 37 { + return (1904459, 1888141, 1938256); + } + if slot == 38 { + return (1955930, 1938256, 1988298); + } + if slot == 39 { + return (2007402, 1988298, 2038265); + } + if slot == 40 { + return (2058871, 2038265, 2088156); + } + if slot == 41 { + return (2110346, 2088156, 2137968); + } + if slot == 42 { + return (2161818, 2137968, 2187700); + } + if slot == 43 { + return (2213290, 2187700, 2237349); + } + if slot == 44 { + return (2264762, 2237349, 2286914); + } + if slot == 45 { + return (2316233, 2286914, 2336392); + } + if slot == 46 { + return (2367705, 2336392, 2385783); + } + if slot == 47 { + return (2419177, 2385783, 2435084); + } + } else { + if slot == 48 { + return (2470649, 2435084, 2484294); + } + if slot == 49 { + return (2522121, 2484294, 2533410); + } + if slot == 50 { + return (2573593, 2533410, 2582430); + } + if slot == 51 { + return (2625065, 2582430, 2631353); + } + if slot == 52 { + return (2676536, 2631353, 2680177); + } + if slot == 53 { + return (2728008, 2680177, 2728901); + } + if slot == 54 { + return (2779480, 2728901, 2777521); + } + if slot == 55 { + return (2830952, 2777521, 2826037); + } + if slot == 56 { + return (2882424, 2826037, 2874446); + } + if slot == 57 { + return (2933896, 2874446, 2922748); + } + if slot == 58 { + return (2985368, 2922748, 2970939); + } + if slot == 59 { + return (3036839, 2970939, 3019018); + } + if slot == 60 { + return (3088311, 3019018, 3066984); + } + if slot == 61 { + return (3139783, 3066984, 3114834); + } + if slot == 62 { + return (3191255, 3114834, 3162567); + } + if slot == 63 { + return (3242727, 3162567, 3210181); + } + } + } + } else { + if slot < 96 { + if slot < 80 { + if slot == 64 { + return (3294199, 3210181, 3257674); + } + if slot == 65 { + return (3345671, 3257674, 3305045); + } + if slot == 66 { + return (3397142, 3305045, 3352291); + } + if slot == 67 { + return (3448614, 3352291, 3399411); + } + if slot == 68 { + return (3500086, 3399411, 3446402); + } + if slot == 69 { + return (3551558, 3446402, 3493264); + } + if slot == 70 { + return (3603030, 3493264, 3539995); + } + if slot == 71 { + return (3654502, 3539995, 3586592); + } + if slot == 72 { + return (3705973, 3586592, 3633054); + } + if slot == 73 { + return (3757445, 3633054, 3679380); + } + if slot == 74 { + return (3808917, 3679380, 3725567); + } + if slot == 75 { + return (3860389, 3725567, 3771613); + } + if slot == 76 { + return (3911861, 3771613, 3817518); + } + if slot == 77 { + return (3963333, 3817518, 3863279); + } + if slot == 78 { + return (4014805, 3863279, 3908894); + } + if slot == 79 { + return (4066276, 3908894, 3954362); + } + } else { + if slot == 80 { + return (4117751, 3954362, 3999682); + } + if slot == 81 { + return (4169220, 3999682, 4044851); + } + if slot == 82 { + return (4220692, 4044851, 4089867); + } + if slot == 83 { + return (4272164, 4089867, 4134730); + } + if slot == 84 { + return (4323636, 4134730, 4179437); + } + if slot == 85 { + return (4375108, 4179437, 4223986); + } + if slot == 86 { + return (4426579, 4223986, 4268377); + } + if slot == 87 { + return (4478051, 4268377, 4312606); + } + if slot == 88 { + return (4529523, 4312606, 4356674); + } + if slot == 89 { + return (4580995, 4356674, 4400577); + } + if slot == 90 { + return (4632474, 4400577, 4444315); + } + if slot == 91 { + return (4683939, 4444315, 4487885); + } + if slot == 92 { + return (4735411, 4487885, 4531287); + } + if slot == 93 { + return (4786882, 4531287, 4574518); + } + if slot == 94 { + return (4838354, 4574518, 4617576); + } + if slot == 95 { + return (4889826, 4617576, 4660461); + } + } + } else { + if slot < 112 { + if slot == 96 { + return (4941298, 4660461, 4703170); + } + if slot == 97 { + return (4992770, 4703170, 4745702); + } + if slot == 98 { + return (5044242, 4745702, 4788056); + } + if slot == 99 { + return (5095714, 4788056, 4830229); + } + if slot == 100 { + return (5147227, 4830229, 4872221); + } + if slot == 101 { + return (5198657, 4872221, 4914029); + } + if slot == 102 { + return (5250129, 4914029, 4955652); + } + if slot == 103 { + return (5301601, 4955652, 4997088); + } + if slot == 104 { + return (5353073, 4997088, 5038336); + } + if slot == 105 { + return (5404545, 5038336, 5079395); + } + if slot == 106 { + return (5456017, 5079395, 5120262); + } + if slot == 107 { + return (5507488, 5120262, 5160937); + } + if slot == 108 { + return (5558960, 5160937, 5201417); + } + if slot == 109 { + return (5610432, 5201417, 5241701); + } + if slot == 110 { + return (5661904, 5241701, 5281788); + } + if slot == 111 { + return (5713376, 5281788, 5321677); + } + } else { + if slot == 112 { + return (5764848, 5321677, 5361364); + } + if slot == 113 { + return (5816320, 5361364, 5400850); + } + if slot == 114 { + return (5867791, 5400850, 5440133); + } + if slot == 115 { + return (5919263, 5440133, 5479211); + } + if slot == 116 { + return (5970735, 5479211, 5518082); + } + if slot == 117 { + return (6022207, 5518082, 5556746); + } + if slot == 118 { + return (6073679, 5556746, 5595201); + } + if slot == 119 { + return (6125151, 5595201, 5633445); + } + if slot == 120 { + return (6176622, 5633445, 5671477); + } + if slot == 121 { + return (6228094, 5671477, 5709295); + } + if slot == 122 { + return (6279566, 5709295, 5746898); + } + if slot == 123 { + return (6331038, 5746898, 5784285); + } + if slot == 124 { + return (6382510, 5784285, 5821455); + } + if slot == 125 { + return (6433982, 5821455, 5858405); + } + if slot == 126 { + return (6485454, 5858405, 5895134); + } + if slot == 127 { + return (6536925, 5895134, 5931642); + } + } + } + } + } else { + if slot < 192 { + if slot < 160 { + if slot < 144 { + if slot == 128 { + return (6588397, 5931642, 5967926); + } + if slot == 129 { + return (6639869, 5967926, 6003985); + } + if slot == 130 { + return (6691345, 6003985, 6039819); + } + if slot == 131 { + return (6742813, 6039819, 6075425); + } + if slot == 132 { + return (6794285, 6075425, 6110802); + } + if slot == 133 { + return (6845757, 6110802, 6145949); + } + if slot == 134 { + return (6897228, 6145949, 6180865); + } + if slot == 135 { + return (6948700, 6180865, 6215549); + } + if slot == 136 { + return (7000172, 6215549, 6249998); + } + if slot == 137 { + return (7051644, 6249998, 6284212); + } + if slot == 138 { + return (7103116, 6284212, 6318189); + } + if slot == 139 { + return (7154588, 6318189, 6351928); + } + if slot == 140 { + return (7206060, 6351928, 6385428); + } + if slot == 141 { + return (7257531, 6385428, 6418688); + } + if slot == 142 { + return (7309003, 6418688, 6451706); + } + if slot == 143 { + return (7360475, 6451706, 6484482); + } + } else { + if slot == 144 { + return (7411947, 6484482, 6517013); + } + if slot == 145 { + return (7463419, 6517013, 6549299); + } + if slot == 146 { + return (7514891, 6549299, 6581338); + } + if slot == 147 { + return (7566363, 6581338, 6613129); + } + if slot == 148 { + return (7617834, 6613129, 6644672); + } + if slot == 149 { + return (7669306, 6644672, 6675964); + } + if slot == 150 { + return (7720780, 6675964, 6707005); + } + if slot == 151 { + return (7772250, 6707005, 6737793); + } + if slot == 152 { + return (7823722, 6737793, 6768328); + } + if slot == 153 { + return (7875194, 6768328, 6798608); + } + if slot == 154 { + return (7926666, 6798608, 6828632); + } + if slot == 155 { + return (7978137, 6828632, 6858399); + } + if slot == 156 { + return (8029609, 6858399, 6887907); + } + if slot == 157 { + return (8081081, 6887907, 6917156); + } + if slot == 158 { + return (8132553, 6917156, 6946145); + } + if slot == 159 { + return (8184025, 6946145, 6974873); + } + if slot == 160 { + return (8235503, 6974873, 7003337); + } + } + } else { + if slot < 176 { + if slot == 161 { + return (8286968, 7003337, 7031538); + } + if slot == 162 { + return (8338440, 7031538, 7059475); + } + if slot == 163 { + return (8389912, 7059475, 7087145); + } + if slot == 164 { + return (8441384, 7087145, 7114549); + } + if slot == 165 { + return (8492856, 7114549, 7141685); + } + if slot == 166 { + return (8544328, 7141685, 7168552); + } + if slot == 167 { + return (8595800, 7168552, 7195149); + } + if slot == 168 { + return (8647271, 7195149, 7221475); + } + if slot == 169 { + return (8698743, 7221475, 7247530); + } + if slot == 170 { + return (8750215, 7247530, 7273311); + } + if slot == 171 { + return (8801687, 7273311, 7298819); + } + if slot == 172 { + return (8853159, 7298819, 7324052); + } + if slot == 173 { + return (8904631, 7324052, 7349009); + } + if slot == 174 { + return (8956103, 7349009, 7373689); + } + if slot == 175 { + return (9007574, 7373689, 7398092); + } + } else { + if slot == 176 { + return (9059046, 7398092, 7422216); + } + if slot == 177 { + return (9110518, 7422216, 7446061); + } + if slot == 178 { + return (9161990, 7446061, 7469625); + } + if slot == 179 { + return (9213462, 7469625, 7492909); + } + if slot == 180 { + return (9264934, 7492909, 7515910); + } + if slot == 181 { + return (9316406, 7515910, 7538628); + } + if slot == 182 { + return (9367877, 7538628, 7561062); + } + if slot == 183 { + return (9419349, 7561062, 7583212); + } + if slot == 184 { + return (9470821, 7583212, 7605076); + } + if slot == 185 { + return (9522293, 7605076, 7626654); + } + if slot == 186 { + return (9573765, 7626654, 7647945); + } + if slot == 187 { + return (9625237, 7647945, 7668947); + } + if slot == 188 { + return (9676709, 7668947, 7689661); + } + if slot == 189 { + return (9728180, 7689661, 7710086); + } + if slot == 190 { + return (9779651, 7710086, 7730220); + } + if slot == 191 { + return (9831124, 7730220, 7750063); + } + } + } + } else { + if slot < 224 { + if slot < 208 { + if slot == 192 { + return (9882596, 7750063, 7769615); + } + if slot == 193 { + return (9934068, 7769615, 7788874); + } + if slot == 194 { + return (9985540, 7788874, 7807839); + } + if slot == 195 { + return (10037012, 7807839, 7826511); + } + if slot == 196 { + return (10088483, 7826511, 7844888); + } + if slot == 197 { + return (10139955, 7844888, 7862970); + } + if slot == 198 { + return (10191427, 7862970, 7880755); + } + if slot == 199 { + return (10242899, 7880755, 7898244); + } + if slot == 200 { + return (10294373, 7898244, 7915436); + } + if slot == 201 { + return (10345843, 7915436, 7932329); + } + if slot == 202 { + return (10397315, 7932329, 7948924); + } + if slot == 203 { + return (10448786, 7948924, 7965220); + } + if slot == 204 { + return (10500258, 7965220, 7981215); + } + if slot == 205 { + return (10551730, 7981215, 7996911); + } + if slot == 206 { + return (10603202, 7996911, 8012305); + } + if slot == 207 { + return (10654674, 8012305, 8027397); + } + } else { + if slot == 208 { + return (10706146, 8027397, 8042188); + } + if slot == 209 { + return (10757617, 8042188, 8056675); + } + if slot == 210 { + return (10809089, 8056675, 8070859); + } + if slot == 211 { + return (10860561, 8070859, 8084740); + } + if slot == 212 { + return (10912033, 8084740, 8098316); + } + if slot == 213 { + return (10963505, 8098316, 8111587); + } + if slot == 214 { + return (11014977, 8111587, 8124552); + } + if slot == 215 { + return (11066449, 8124552, 8137212); + } + if slot == 216 { + return (11117920, 8137212, 8149565); + } + if slot == 217 { + return (11169392, 8149565, 8161612); + } + if slot == 218 { + return (11220864, 8161612, 8173351); + } + if slot == 219 { + return (11272336, 8173351, 8184783); + } + if slot == 220 { + return (11323808, 8184783, 8195906); + } + if slot == 221 { + return (11375280, 8195906, 8206721); + } + if slot == 222 { + return (11426752, 8206721, 8217227); + } + if slot == 223 { + return (11478223, 8217227, 8227423); + } + } + } else { + if slot < 240 { + if slot == 224 { + return (11529695, 8227423, 8237310); + } + if slot == 225 { + return (11581167, 8237310, 8246887); + } + if slot == 226 { + return (11632639, 8246887, 8256153); + } + if slot == 227 { + return (11684111, 8256153, 8265108); + } + if slot == 228 { + return (11735583, 8265108, 8273752); + } + if slot == 229 { + return (11787055, 8273752, 8282085); + } + if slot == 230 { + return (11838531, 8282085, 8290105); + } + if slot == 231 { + return (11889998, 8290105, 8297814); + } + if slot == 232 { + return (11941470, 8297814, 8305210); + } + if slot == 233 { + return (11992942, 8305210, 8312294); + } + if slot == 234 { + return (12044414, 8312294, 8319064); + } + if slot == 235 { + return (12095886, 8319064, 8325522); + } + if slot == 236 { + return (12147358, 8325522, 8331666); + } + if slot == 237 { + return (12198829, 8331666, 8337496); + } + if slot == 238 { + return (12250301, 8337496, 8343012); + } + if slot == 239 { + return (12301773, 8343012, 8348215); + } + } else { + if slot == 240 { + return (12353244, 8348215, 8353102); + } + if slot == 241 { + return (12404717, 8353102, 8357676); + } + if slot == 242 { + return (12456189, 8357676, 8361935); + } + if slot == 243 { + return (12507661, 8361935, 8365879); + } + if slot == 244 { + return (12559132, 8365879, 8369508); + } + if slot == 245 { + return (12610604, 8369508, 8372822); + } + if slot == 246 { + return (12662076, 8372822, 8375820); + } + if slot == 247 { + return (12713548, 8375820, 8378504); + } + if slot == 248 { + return (12765020, 8378504, 8380871); + } + if slot == 249 { + return (12816492, 8380871, 8382924); + } + if slot == 250 { + return (12867964, 8382924, 8384660); + } + if slot == 251 { + return (12919435, 8384660, 8386082); + } + if slot == 252 { + return (12970907, 8386082, 8387187); + } + if slot == 253 { + return (13022379, 8387187, 8387976); + } + if slot == 254 { + return (13073851, 8387976, 8388450); + } + } + } + } + } + + return (13125323, 8388450, 8388608); +} + +fn atan(a: u64) -> (u64, u64, u64) { + let slot = a / 58720; + + if slot == 0 { + return (0, 0, 58719); + } + if slot == 1 { + return (58720, 58719, 117433); + } + if slot == 2 { + return (117441, 117433, 176135); + } + if slot == 3 { + return (176161, 176135, 234820); + } + if slot == 4 { + return (234881, 234820, 293481); + } + if slot == 5 { + return (293601, 293481, 352115); + } + if slot == 6 { + return (352322, 352115, 410713); + } + if slot == 7 { + return (411042, 410713, 469272); + } + if slot == 8 { + return (469762, 469272, 527785); + } + if slot == 9 { + return (528482, 527785, 586246); + } + if slot == 10 { + return (587201, 586246, 644651); + } + if slot == 11 { + return (645923, 644651, 702993); + } + if slot == 12 { + return (704643, 702993, 761267); + } + if slot == 13 { + return (763363, 761267, 819467); + } + if slot == 14 { + return (822084, 819467, 877588); + } + if slot == 15 { + return (880804, 877588, 935625); + } + if slot == 16 { + return (939524, 935625, 993572); + } + if slot == 17 { + return (998244, 993572, 1051424); + } + if slot == 18 { + return (1056965, 1051424, 1109175); + } + if slot == 19 { + return (1115685, 1109175, 1166821); + } + if slot == 20 { + return (1174411, 1166821, 1224357); + } + if slot == 21 { + return (1233125, 1224357, 1281776); + } + if slot == 22 { + return (1291846, 1281776, 1339075); + } + if slot == 23 { + return (1350566, 1339075, 1396248); + } + if slot == 24 { + return (1409286, 1396248, 1453290); + } + if slot == 25 { + return (1468006, 1453290, 1510197); + } + if slot == 26 { + return (1526727, 1510197, 1566964); + } + if slot == 27 { + return (1585447, 1566964, 1623585); + } + if slot == 28 { + return (1644167, 1623585, 1680058); + } + if slot == 29 { + return (1702887, 1680058, 1736376); + } + if slot == 30 { + return (1761612, 1736376, 1792537); + } + if slot == 31 { + return (1820328, 1792537, 1848534); + } + if slot == 32 { + return (1879048, 1848534, 1904364); + } + if slot == 33 { + return (1937768, 1904364, 1960024); + } + if slot == 34 { + return (1996489, 1960024, 2015508); + } + if slot == 35 { + return (2055209, 2015508, 2070813); + } + if slot == 36 { + return (2113929, 2070813, 2125935); + } + if slot == 37 { + return (2172649, 2125935, 2180869); + } + if slot == 38 { + return (2231370, 2180869, 2235613); + } + if slot == 39 { + return (2290090, 2235613, 2290163); + } + if slot == 40 { + return (2348813, 2290163, 2344515); + } + if slot == 41 { + return (2407530, 2344515, 2398665); + } + if slot == 42 { + return (2466251, 2398665, 2452611); + } + if slot == 43 { + return (2524971, 2452611, 2506348); + } + if slot == 44 { + return (2583691, 2506348, 2559875); + } + if slot == 45 { + return (2642412, 2559875, 2613187); + } + if slot == 46 { + return (2701132, 2613187, 2666281); + } + if slot == 47 { + return (2759852, 2666281, 2719156); + } + if slot == 48 { + return (2818572, 2719156, 2771807); + } + if slot == 49 { + return (2877293, 2771807, 2824233); + } + if slot == 50 { + return (2936014, 2824233, 2876431); + } + if slot == 51 { + return (2994733, 2876431, 2928397); + } + if slot == 52 { + return (3053453, 2928397, 2980130); + } + if slot == 53 { + return (3112174, 2980130, 3031628); + } + if slot == 54 { + return (3170894, 3031628, 3082888); + } + if slot == 55 { + return (3229614, 3082888, 3133907); + } + if slot == 56 { + return (3288334, 3133907, 3184685); + } + if slot == 57 { + return (3347055, 3184685, 3235218); + } + if slot == 58 { + return (3405775, 3235218, 3285506); + } + if slot == 59 { + return (3464495, 3285506, 3335545); + } + if slot == 60 { + return (3523224, 3335545, 3385336); + } + if slot == 61 { + return (3581936, 3385336, 3434875); + } + if slot == 62 { + return (3640656, 3434875, 3484161); + } + if slot == 63 { + return (3699376, 3484161, 3533193); + } + if slot == 64 { + return (3758096, 3533193, 3581970); + } + if slot == 65 { + return (3816817, 3581970, 3630491); + } + if slot == 66 { + return (3875537, 3630491, 3678753); + } + if slot == 67 { + return (3934257, 3678753, 3726756); + } + if slot == 68 { + return (3992977, 3726756, 3774499); + } + if slot == 69 { + return (4051698, 3774499, 3821981); + } + if slot == 70 { + return (4110418, 3821981, 3869201); + } + if slot == 71 { + return (4169138, 3869201, 3916159); + } + if slot == 72 { + return (4227858, 3916159, 3962853); + } + if slot == 73 { + return (4286579, 3962853, 4009282); + } + if slot == 74 { + return (4345299, 4009282, 4055447); + } + if slot == 75 { + return (4404019, 4055447, 4101347); + } + if slot == 76 { + return (4462739, 4101347, 4146981); + } + if slot == 77 { + return (4521460, 4146981, 4192350); + } + if slot == 78 { + return (4580180, 4192350, 4237451); + } + if slot == 79 { + return (4638900, 4237451, 4282286); + } + if slot == 80 { + return (4697620, 4282286, 4326855); + } + if slot == 81 { + return (4756341, 4326855, 4371156); + } + if slot == 82 { + return (4815061, 4371156, 4415191); + } + if slot == 83 { + return (4873781, 4415191, 4458958); + } + if slot == 84 { + return (4932502, 4458958, 4502459); + } + if slot == 85 { + return (4991222, 4502459, 4545693); + } + if slot == 86 { + return (5049942, 4545693, 4588660); + } + if slot == 87 { + return (5108662, 4588660, 4631361); + } + if slot == 88 { + return (5167383, 4631361, 4673795); + } + if slot == 89 { + return (5226103, 4673795, 4715964); + } + if slot == 90 { + return (5284823, 4715964, 4757868); + } + if slot == 91 { + return (5343543, 4757868, 4799506); + } + if slot == 92 { + return (5402264, 4799506, 4840880); + } + if slot == 93 { + return (5460984, 4840880, 4881990); + } + if slot == 94 { + return (5519704, 4881990, 4922837); + } + if slot == 95 { + return (5578424, 4922837, 4963420); + } + if slot == 96 { + return (5637145, 4963420, 5003742); + } + if slot == 97 { + return (5695865, 5003742, 5043802); + } + if slot == 98 { + return (5754585, 5043802, 5083601); + } + + return (5813305, 5083601, 5123141); +} diff --git a/src/numbers/fixed_point/implementations/fp8x23wide/math/trig.cairo b/src/numbers/fixed_point/implementations/fp8x23wide/math/trig.cairo new file mode 100644 index 000000000..025b79bb2 --- /dev/null +++ b/src/numbers/fixed_point/implementations/fp8x23wide/math/trig.cairo @@ -0,0 +1,448 @@ +use debug::PrintTrait; +use integer::{u64_safe_divmod, u64_as_non_zero}; +use option::OptionTrait; + +use orion::numbers::fixed_point::implementations::fp8x23wide::math::lut; +use orion::numbers::fixed_point::implementations::fp8x23wide::core::{ + HALF, ONE, TWO, FP8x23W, FP8x23WImpl, FP8x23WAdd, FP8x23WSub, FP8x23WMul, FP8x23WDiv, + FP8x23WIntoFelt252, FixedTrait +}; + +// CONSTANTS + +const TWO_PI: u64 = 52707178; +const PI: u64 = 26353589; +const HALF_PI: u64 = 13176795; + +// PUBLIC + +// Calculates arccos(a) for -1 <= a <= 1 (fixed point) +// arccos(a) = arcsin(sqrt(1 - a^2)) - arctan identity has discontinuity at zero +fn acos(a: FP8x23W) -> FP8x23W { + let asin_arg = (FixedTrait::ONE() - a * a).sqrt(); // will fail if a > 1 + let asin_res = asin(asin_arg); + + if (a.sign) { + return FixedTrait::new(PI, false) - asin_res; + } else { + return asin_res; + } +} + +fn acos_fast(a: FP8x23W) -> FP8x23W { + let asin_arg = (FixedTrait::ONE() - a * a).sqrt(); // will fail if a > 1 + let asin_res = asin_fast(asin_arg); + + if (a.sign) { + return FixedTrait::new(PI, false) - asin_res; + } else { + return asin_res; + } +} + +// Calculates arcsin(a) for -1 <= a <= 1 (fixed point) +// arcsin(a) = arctan(a / sqrt(1 - a^2)) +fn asin(a: FP8x23W) -> FP8x23W { + if (a.mag == ONE) { + return FixedTrait::new(HALF_PI, a.sign); + } + + let div = (FixedTrait::ONE() - a * a).sqrt(); // will fail if a > 1 + return atan(a / div); +} + +fn asin_fast(a: FP8x23W) -> FP8x23W { + if (a.mag == ONE) { + return FixedTrait::new(HALF_PI, a.sign); + } + + let div = (FixedTrait::ONE() - a * a).sqrt(); // will fail if a > 1 + return atan_fast(a / div); +} + +// Calculates arctan(a) (fixed point) +// See https://stackoverflow.com/a/50894477 for range adjustments +fn atan(a: FP8x23W) -> FP8x23W { + let mut at = a.abs(); + let mut shift = false; + let mut invert = false; + + // Invert value when a > 1 + if (at.mag > ONE) { + at = FixedTrait::ONE() / at; + invert = true; + } + + // Account for lack of precision in polynomaial when a > 0.7 + if (at.mag > 5872026) { + let sqrt3_3 = FixedTrait::new(4843165, false); // sqrt(3) / 3 + at = (at - sqrt3_3) / (FixedTrait::ONE() + at * sqrt3_3); + shift = true; + } + + let r10 = FixedTrait::new(15363, true) * at; + let r9 = (r10 + FixedTrait::new(392482, true)) * at; + let r8 = (r9 + FixedTrait::new(1629064, false)) * at; + let r7 = (r8 + FixedTrait::new(2197820, true)) * at; + let r6 = (r7 + FixedTrait::new(366693, false)) * at; + let r5 = (r6 + FixedTrait::new(1594324, false)) * at; + let r4 = (r5 + FixedTrait::new(11519, false)) * at; + let r3 = (r4 + FixedTrait::new(2797104, true)) * at; + let r2 = (r3 + FixedTrait::new(34, false)) * at; + let mut res = (r2 + FixedTrait::new(8388608, false)) * at; + + // Adjust for sign change, inversion, and shift + if (shift) { + res = res + FixedTrait::new(4392265, false); // pi / 6 + } + + if (invert) { + res = res - FixedTrait::new(HALF_PI, false); + } + + return FixedTrait::new(res.mag, a.sign); +} + +fn atan_fast(a: FP8x23W) -> FP8x23W { + let mut at = a.abs(); + let mut shift = false; + let mut invert = false; + + // Invert value when a > 1 + if (at.mag > ONE) { + at = FixedTrait::ONE() / at; + invert = true; + } + + // Account for lack of precision in polynomaial when a > 0.7 + if (at.mag > 5872026) { + let sqrt3_3 = FixedTrait::new(4843165, false); // sqrt(3) / 3 + at = (at - sqrt3_3) / (FixedTrait::ONE() + at * sqrt3_3); + shift = true; + } + + let (start, low, high) = lut::atan(at.mag); + let partial_step = FixedTrait::new(at.mag - start, false) / FixedTrait::new(58720, false); + let mut res = partial_step * FixedTrait::new(high - low, false) + FixedTrait::new(low, false); + + // Adjust for sign change, inversion, and shift + if (shift) { + res = res + FixedTrait::new(4392265, false); // pi / 6 + } + + if (invert) { + res = res - FixedTrait::::new(HALF_PI, false); + } + + return FixedTrait::new(res.mag, a.sign); +} + +// Calculates cos(a) with a in radians (fixed point) +fn cos(a: FP8x23W) -> FP8x23W { + return sin(FixedTrait::new(HALF_PI, false) - a); +} + +fn cos_fast(a: FP8x23W) -> FP8x23W { + return sin_fast(FixedTrait::new(HALF_PI, false) - a); +} + +fn sin(a: FP8x23W) -> FP8x23W { + let a1 = a.mag % TWO_PI; + let (whole_rem, partial_rem) = u64_safe_divmod(a1, u64_as_non_zero(PI)); + let a2 = FixedTrait::new(partial_rem, false); + let partial_sign = whole_rem == 1; + + let loop_res = a2 * _sin_loop(a2, 7, FixedTrait::ONE()); + return FixedTrait::new(loop_res.mag, a.sign ^ partial_sign && loop_res.mag != 0); +} + +fn sin_fast(a: FP8x23W) -> FP8x23W { + let a1 = a.mag % TWO_PI; + let (whole_rem, mut partial_rem) = u64_safe_divmod(a1, u64_as_non_zero(PI)); + let partial_sign = whole_rem == 1; + + if partial_rem >= HALF_PI { + partial_rem = PI - partial_rem; + } + + let (start, low, high) = lut::sin(partial_rem); + let partial_step = FixedTrait::new(partial_rem - start, false) / FixedTrait::new(51472, false); + let res = partial_step * (FixedTrait::new(high, false) - FixedTrait::new(low, false)) + + FixedTrait::::new(low, false); + + return FixedTrait::new(res.mag, a.sign ^ partial_sign && res.mag != 0); +} + +// Calculates tan(a) with a in radians (fixed point) +fn tan(a: FP8x23W) -> FP8x23W { + let sinx = sin(a); + let cosx = cos(a); + assert(cosx.mag != 0, 'tan undefined'); + return sinx / cosx; +} + +fn tan_fast(a: FP8x23W) -> FP8x23W { + let sinx = sin_fast(a); + let cosx = cos_fast(a); + assert(cosx.mag != 0, 'tan undefined'); + return sinx / cosx; +} + +// Helper function to calculate Taylor series for sin +fn _sin_loop(a: FP8x23W, i: u64, acc: FP8x23W) -> FP8x23W { + let div = (2 * i + 2) * (2 * i + 3); + let term = a * a * acc / FixedTrait::new_unscaled(div, false); + let new_acc = FixedTrait::ONE() - term; + + if (i == 0) { + return new_acc; + } + + return _sin_loop(a, i - 1, new_acc); +} + +// Tests -------------------------------------------------------------------------------------------------------------- + +use traits::Into; + +use orion::numbers::fixed_point::implementations::fp8x23wide::helpers::{ + assert_precise, assert_relative +}; + +#[test] +#[available_gas(3000000)] +fn test_acos() { + let error = Option::Some(84); // 1e-5 + + let a = FixedTrait::ONE(); + assert(acos(a).into() == 0, 'invalid one'); + + let a = FixedTrait::new(ONE / 2, false); + assert_relative(acos(a), 8784530, 'invalid half', error); // 1.0471975506263043 + + let a = FixedTrait::ZERO(); + assert_relative(acos(a), HALF_PI.into(), 'invalid zero', Option::None(())); // PI / 2 + + let a = FixedTrait::new(ONE / 2, true); + assert_relative(acos(a), 17569060, 'invalid neg half', error); // 2.094395102963489 + + let a = FixedTrait::new(ONE, true); + assert_relative(acos(a), PI.into(), 'invalid neg one', Option::None(())); // PI +} + +#[test] +#[available_gas(3000000)] +fn test_acos_fast() { + let error = Option::Some(84); // 1e-5 + + let a = FixedTrait::ONE(); + assert(acos_fast(a).into() == 0, 'invalid one'); + + let a = FixedTrait::new(ONE / 2, false); + assert_relative(acos_fast(a), 8784530, 'invalid half', error); // 1.0471975506263043 + + let a = FixedTrait::ZERO(); + assert_relative(acos_fast(a), HALF_PI.into(), 'invalid zero', Option::None(())); // PI / 2 + + let a = FixedTrait::new(ONE / 2, true); + assert_relative(acos_fast(a), 17569060, 'invalid neg half', error); // 2.094395102963489 + + let a = FixedTrait::new(ONE, true); + assert_relative(acos_fast(a), PI.into(), 'invalid neg one', Option::None(())); // PI +} + +#[test] +#[should_panic] +#[available_gas(1000000)] +fn test_acos_fail() { + let a = FixedTrait::new(2 * ONE, true); + acos(a); +} + +#[test] +#[available_gas(1400000)] +fn test_atan_fast() { + let error = Option::Some(84); // 1e-5 + + let a = FixedTrait::new(2 * ONE, false); + assert_relative(atan_fast(a), 9287437, 'invalid two', error); + + let a = FixedTrait::ONE(); + assert_relative(atan_fast(a), 6588397, 'invalid one', error); + + let a = FixedTrait::new(ONE / 2, false); + assert_relative(atan_fast(a), 3889358, 'invalid half', error); + + let a = FixedTrait::ZERO(); + assert(atan_fast(a).into() == 0, 'invalid zero'); + + let a = FixedTrait::new(ONE / 2, true); + assert_relative(atan_fast(a), -3889358, 'invalid neg half', error); + + let a = FixedTrait::new(ONE, true); + assert_relative(atan_fast(a), -6588397, 'invalid neg one', error); + + let a = FixedTrait::new(2 * ONE, true); + assert_relative(atan_fast(a), -9287437, 'invalid neg two', error); +} + +#[test] +#[available_gas(2600000)] +fn test_atan() { + let a = FixedTrait::new(2 * ONE, false); + assert_relative(atan(a), 9287437, 'invalid two', Option::None(())); + + let a = FixedTrait::ONE(); + assert_relative(atan(a), 6588397, 'invalid one', Option::None(())); + + let a = FixedTrait::new(ONE / 2, false); + assert_relative(atan(a), 3889358, 'invalid half', Option::None(())); + + let a = FixedTrait::ZERO(); + assert(atan(a).into() == 0, 'invalid zero'); + + let a = FixedTrait::new(ONE / 2, true); + assert_relative(atan(a), -3889358, 'invalid neg half', Option::None(())); + + let a = FixedTrait::new(ONE, true); + assert_relative(atan(a), -6588397, 'invalid neg one', Option::None(())); + + let a = FixedTrait::new(2 * ONE, true); + assert_relative(atan(a), -9287437, 'invalid neg two', Option::None(())); +} + +#[test] +#[available_gas(3000000)] +fn test_asin() { + let error = Option::Some(84); // 1e-5 + + let a = FixedTrait::ONE(); + assert_relative(asin(a), HALF_PI.into(), 'invalid one', Option::None(())); // PI / 2 + + let a = FixedTrait::new(ONE / 2, false); + assert_relative(asin(a), 4392265, 'invalid half', error); + + let a = FixedTrait::ZERO(); + assert_precise(asin(a), 0, 'invalid zero', Option::None(())); + + let a = FixedTrait::new(ONE / 2, true); + assert_relative(asin(a), -4392265, 'invalid neg half', error); + + let a = FixedTrait::new(ONE, true); + assert_relative(asin(a), -HALF_PI.into(), 'invalid neg one', Option::None(())); // -PI / 2 +} + +#[test] +#[should_panic] +#[available_gas(1000000)] +fn test_asin_fail() { + let a = FixedTrait::new(2 * ONE, false); + asin(a); +} + +#[test] +#[available_gas(6000000)] +fn test_cos() { + let a = FixedTrait::new(HALF_PI, false); + assert(cos(a).into() == 0, 'invalid half pi'); + + let a = FixedTrait::new(HALF_PI / 2, false); + assert_relative(cos(a), 5931642, 'invalid quarter pi', Option::None(())); // 0.7071067811865475 + + let a = FixedTrait::new(PI, false); + assert_relative(cos(a), -1 * ONE.into(), 'invalid pi', Option::None(())); + + let a = FixedTrait::new(HALF_PI, true); + assert_precise(cos(a), 0, 'invalid neg half pi', Option::None(())); + + let a = FixedTrait::new_unscaled(17, false); + assert_relative(cos(a), -2308239, 'invalid 17', Option::None(())); // -0.2751631780463348 + + let a = FixedTrait::new_unscaled(17, true); + assert_relative(cos(a), -2308236, 'invalid -17', Option::None(())); // -0.2751631780463348 +} + +#[test] +#[available_gas(6000000)] +fn test_cos_fast() { + let error = Option::Some(84); // 1e-5 + + let a = FixedTrait::new(HALF_PI, false); + assert(cos_fast(a).into() == 0, 'invalid half pi'); + + let a = FixedTrait::new(HALF_PI / 2, false); + assert_precise(cos_fast(a), 5931642, 'invalid quarter pi', error); // 0.7071067811865475 + + let a = FixedTrait::new(PI, false); + assert_precise(cos_fast(a), -1 * ONE.into(), 'invalid pi', error); + + let a = FixedTrait::new(HALF_PI, true); + assert_precise(cos(a), 0, 'invalid neg half pi', Option::None(())); + + let a = FixedTrait::new_unscaled(17, false); + assert_precise(cos_fast(a), -2308239, 'invalid 17', error); // -0.2751631780463348 +} + +#[test] +#[available_gas(6000000)] +fn test_sin() { + let a = FixedTrait::new(HALF_PI, false); + assert_precise(sin(a), ONE.into(), 'invalid half pi', Option::None(())); + + let a = FixedTrait::new(HALF_PI / 2, false); + assert_precise(sin(a), 5931642, 'invalid quarter pi', Option::None(())); // 0.7071067811865475 + + let a = FixedTrait::new(PI, false); + assert(sin(a).into() == 0, 'invalid pi'); + + let a = FixedTrait::new(HALF_PI, true); + assert_precise( + sin(a), -ONE.into(), 'invalid neg half pi', Option::None(()) + ); // 0.9999999999939766 + + let a = FixedTrait::new_unscaled(17, false); + assert_precise(sin(a), -8064787, 'invalid 17', Option::None(())); // -0.9613974918793389 + + let a = FixedTrait::new_unscaled(17, true); + assert_precise(sin(a), 8064787, 'invalid -17', Option::None(())); // 0.9613974918793389 +} + +#[test] +#[available_gas(1000000)] +fn test_sin_fast() { + let error = Option::Some(84); // 1e-5 + + let a = FixedTrait::new(HALF_PI, false); + assert_precise(sin_fast(a), ONE.into(), 'invalid half pi', error); + + let a = FixedTrait::new(HALF_PI / 2, false); + assert_precise(sin_fast(a), 5931642, 'invalid quarter pi', error); // 0.7071067811865475 + + let a = FixedTrait::new(PI, false); + assert(sin_fast(a).into() == 0, 'invalid pi'); + + let a = FixedTrait::new(HALF_PI, true); + assert_precise(sin_fast(a), -ONE.into(), 'invalid neg half pi', error); // 0.9999999999939766 + + let a = FixedTrait::new_unscaled(17, false); + assert_precise(sin_fast(a), -8064787, 'invalid 17', error); // -0.9613974918793389 + + let a = FixedTrait::new_unscaled(17, true); + assert_precise(sin_fast(a), 8064787, 'invalid -17', error); // 0.9613974918793389 +} + +#[test] +#[available_gas(8000000)] +fn test_tan() { + let a = FixedTrait::new(HALF_PI / 2, false); + assert_precise(tan(a), ONE.into(), 'invalid quarter pi', Option::None(())); + + let a = FixedTrait::new(PI, false); + assert_precise(tan(a), 0, 'invalid pi', Option::None(())); + + let a = FixedTrait::new_unscaled(17, false); + assert_precise(tan(a), 29309069, 'invalid 17', Option::None(())); // 3.493917677159002 + + let a = FixedTrait::new_unscaled(17, true); + assert_precise(tan(a), -29309106, 'invalid -17', Option::None(())); // -3.493917677159002 +}