From 5016341f0baf9f63d0add6a8bee1906d8b7c7391 Mon Sep 17 00:00:00 2001 From: hnakashima Date: Sun, 21 Jan 2024 18:54:38 +0900 Subject: [PATCH] fix: Fix to use f64::MIN_POSITIVE refactor: Change to use more f64 methods --- src/bachelier.rs | 30 ++++++------ src/constants.rs | 7 ++- src/erf_cody.rs | 76 ++++++++++++++---------------- src/lets_be_rational.rs | 95 +++++++++++++++++++++----------------- src/normal_distribution.rs | 4 +- src/rational_cubic.rs | 4 +- 6 files changed, 109 insertions(+), 107 deletions(-) diff --git a/src/bachelier.rs b/src/bachelier.rs index 435bb9f..f60c873 100644 --- a/src/bachelier.rs +++ b/src/bachelier.rs @@ -1,3 +1,4 @@ +use std::cmp::Ordering; use crate::constants::{ONE_OVER_SQRT_TWO_PI, SQRT_TWO_PI}; use crate::normal_distribution::norm_pdf; @@ -13,7 +14,7 @@ fn phi_tilde_times_x(x: f64) -> f64 { return ONE_OVER_SQRT_TWO_PI + x * (0.5 + x * g); } - if x > 0.0 { + if x.is_sign_positive() { return phi_tilde_times_x(-x) + x; } @@ -22,7 +23,7 @@ fn phi_tilde_times_x(x: f64) -> f64 { return (-0.5 * (x * x)).exp() * g; } - let w = 1.0 / (x * x); + let w = (x * x).recip(); let g = (2.999_999_999_999_991 + w * (2.365_455_662_782_315E2 + w * (6.812_677_344_935_879E3 + w * (8.969_794_159_836_079E4 + w * (5.516_392_059_126_862E5 + w * (1.434_506_112_333_566_2E6 + w * (1.150_498_824_634_488_2E6 + 1.186_760_040_099_769_1E4 * w))))))) / (1.0 + w * (8.384_852_209_273_714E1 + w * (2.655_135_058_780_958E3 + w * (4.055_529_088_467_379E4 + w * (3.166_737_476_299_376_6E5 + w * (1.232_979_595_802_432_2E6 + w * (2.140_981_054_061_905E6 + 1.214_566_780_409_316E6 * w))))))); ONE_OVER_SQRT_TWO_PI * (-0.5 * (x * x)).exp() * w * (1.0 - g * w) } @@ -35,12 +36,12 @@ fn inv_phi_tilde(phi_tilde_star: f64) -> f64 { if phi_tilde_star > 1.0 { return -inv_phi_tilde(1.0 - phi_tilde_star); } - if phi_tilde_star >= 0.0 { + if !phi_tilde_star.is_sign_negative() { return f64::NAN; } let x_bar = if phi_tilde_star < -0.00188203927 { // Equation (2.1) - let g = 1.0 / (phi_tilde_star - 0.5); + let g = (phi_tilde_star - 0.5).recip(); let g2 = g * g; // Equation (2.2) let xi_bar = (0.032114372355 - g2 * (0.016969777977 - g2 * (0.002620733246 - 0.000096066952861 * g2))) / @@ -80,8 +81,7 @@ pub(crate) fn bachelier(forward: f64, strike: f64, sigma: f64, t: f64, q: bool) if s < f64::MIN_POSITIVE { return intrinsic_value(forward, strike, q); } - let theta = if q { 1.0 } else { -1.0 }; - let moneyness = theta * (forward - strike); + let moneyness = if q { forward - strike } else { strike - forward }; let x = moneyness / s; s * phi_tilde_times_x(x) } @@ -91,16 +91,16 @@ pub(crate) fn implied_normal_volatility(price: f64, forward: f64, strike: f64, t return price * SQRT_TWO_PI / t.sqrt(); } let intrinsic = intrinsic_value(forward, strike, q); - if price == intrinsic { - return 0.0; - } - if price < intrinsic { - return f64::NEG_INFINITY; + match price.total_cmp(&intrinsic) { + Ordering::Less => { f64::NEG_INFINITY } + Ordering::Equal => { 0.0 } + Ordering::Greater => { + let absolute_moneyness = (forward - strike).abs(); + let phi_tilde_star = (intrinsic - price) / absolute_moneyness; + let x_star = inv_phi_tilde(phi_tilde_star); + absolute_moneyness / (x_star * t.sqrt()).abs() + } } - let absolute_moneyness = (forward - strike).abs(); - let phi_tilde_star = (intrinsic - price) / absolute_moneyness; - let x_star = inv_phi_tilde(phi_tilde_star); - absolute_moneyness / (x_star * t.sqrt()).abs() } diff --git a/src/constants.rs b/src/constants.rs index 8dd886d..b79a7aa 100644 --- a/src/constants.rs +++ b/src/constants.rs @@ -1,19 +1,18 @@ -pub(crate) const TWO_PI: f64 = std::f64::consts::TAU; pub(crate) const SQRT_PI_OVER_TWO: f64 = 1.253_314_137_315_500_3; pub(crate) const SQRT_TWO_PI: f64 = 2.506_628_274_631_000_7; pub(crate) const SQRT_THREE: f64 = 1.732_050_807_568_877_2; pub(crate) const SQRT_ONE_OVER_THREE: f64 = 0.577_350_269_189_625_7; pub(crate) const TWO_PI_OVER_SQRT_TWENTY_SEVEN: f64 = 1.209_199_576_156_145_2; pub(crate) const SQRT_THREE_OVER_THIRD_ROOT_TWO_PI: f64 = 0.938_643_487_427_383_6; -pub(crate) const PI_OVER_SIX: f64 = std::f64::consts::FRAC_PI_6; pub(crate) const FOURTH_ROOT_DBL_EPSILON: f64 = 0.0001220703125; pub(crate) const SIXTEENTH_ROOT_DBL_EPSILON: f64 = 0.10511205190671433; -pub(crate) const SQRT_DBL_MIN: f64 = 1.4916681462400413e-154; +pub(crate) const SQRT_MIN_POSITIVE: f64 = 1.4916681462400413e-154; pub(crate) const SQRT_DBL_MAX: f64 = 1.3407807929942596e154; +// pub(crate) const SQRT_DBL_MIN: f64 = 1.3407807929942596e154; // Set this to 0 if you want positive results for (positive) denormalised inputs, else to DBL_MIN. // Note that you cannot achieve full machine accuracy from denormalised inputs! pub(crate) const DENORMALISATION_CUTOFF: f64 = 0.0; pub(crate) const ONE_OVER_SQRT_TWO_PI: f64 = 0.3989422804014327; pub(crate) const VOLATILITY_VALUE_TO_SIGNAL_PRICE_IS_BELOW_INTRINSIC: f64 = f64::NEG_INFINITY; pub(crate) const VOLATILITY_VALUE_TO_SIGNAL_PRICE_IS_ABOVE_MAXIMUM: f64 = f64::INFINITY; -pub(crate) const LN_TWO_PI: f64 = 1.8378770664093453; \ No newline at end of file +pub(crate) const HALF_OF_LN_TWO_PI: f64 = 0.918_938_533_204_672_8; diff --git a/src/erf_cody.rs b/src/erf_cody.rs index db75b74..a231fe1 100644 --- a/src/erf_cody.rs +++ b/src/erf_cody.rs @@ -48,14 +48,8 @@ const Q: [f64; 5] = [ 0.002_335_204_976_268_691_8, ]; -const ZERO: f64 = 0.0; -// const HALF: f64 = 0.5; -const ONE: f64 = 1.0; -const TWO: f64 = 2.0; -const FOUR: f64 = 4.0; const SQRPI: f64 = 0.564_189_583_547_756_3; const THRESH: f64 = 0.46875; -const SIXTEN: f64 = 16.0; const XINF: f64 = f64::MAX; const XNEG: f64 = -26.628; const XSMALL: f64 = 1.11e-16; @@ -70,10 +64,10 @@ pub(crate) fn erfc_cody(x: f64) -> f64 { /* Author/date: W. J. Cody, January 8, 1985 */ /* -------------------------------------------------------------------- */ let y = x.abs(); - let mut ysq = ZERO; + let mut ysq = 0.0; let mut xden; let mut xnum; - let mut result = ZERO; + let mut result = 0.0; if y <= THRESH { if y > XSMALL { @@ -88,9 +82,9 @@ pub(crate) fn erfc_cody(x: f64) -> f64 { } result = x * (xnum + A[3]) / (xden + B[3]); - result = ONE - result; + result = 1.0 - result; return result; - } else if y <= FOUR { + } else if y <= 4.0 { xnum = C[8] * y; xden = y; @@ -100,16 +94,16 @@ pub(crate) fn erfc_cody(x: f64) -> f64 { } result = (xnum + C[7]) / (xden + D[7]); - ysq = (y * SIXTEN).trunc() / SIXTEN; + ysq = (y * 16.0).trunc() / 16.0; let del = (y - ysq) * (y + ysq); result *= (-ysq * ysq).exp() * (-del).exp(); } else if y >= XBIG { - if x < ZERO { - result = TWO - result; + if x.is_sign_negative() { + result = 2.0 - result; } return result; } else { - ysq = ONE / (y * y); + ysq = (y * y).recip(); xnum = P[5] * ysq; xden = ysq; @@ -120,12 +114,12 @@ pub(crate) fn erfc_cody(x: f64) -> f64 { result = ysq * (xnum + P[4]) / (xden + Q[4]); result = (SQRPI - result) / y; - ysq = (y * SIXTEN).trunc() / SIXTEN; + ysq = (y * 16.0).trunc() / 16.0; let del = (y - ysq) * (y + ysq); result *= (-ysq * ysq).exp() * (-del).exp(); } - if x < ZERO { - result = TWO - result; + if x.is_sign_negative() { + result = 2.0 - result; } result } @@ -137,10 +131,10 @@ pub(crate) fn erfcx_cody(x: f64) -> f64 { /* Author/date: W. J. Cody, March 30, 1987 */ /* ------------------------------------------------------------------ */ let y = x.abs(); - let mut ysq = ZERO; + let mut ysq = 0.0; let mut xden; let mut xnum; - let mut result = ZERO; + let mut result = 0.0; if y <= THRESH { if y > XSMALL { @@ -155,11 +149,11 @@ pub(crate) fn erfcx_cody(x: f64) -> f64 { } result = x * (xnum + A[3]) / (xden + B[3]); - result = ONE - result; + result = 1.0 - result; result *= ysq.exp(); return result; - } else if y <= FOUR { + } else if y <= 4.0 { xnum = C[8] * y; xden = y; @@ -170,11 +164,11 @@ pub(crate) fn erfcx_cody(x: f64) -> f64 { result = (xnum + C[7]) / (xden + D[7]); } else if y >= XBIG { if y >= XMAX { - if x < ZERO { + if x.is_sign_negative() { if x < XNEG { result = XINF; } else { - let ysq = (x * SIXTEN).trunc() / SIXTEN; + let ysq = (x * 16.0).trunc() / 16.0; let del = (x - ysq) * (x + ysq); let y = (ysq * ysq).exp() * del.exp(); result = (y + y) - result; @@ -183,11 +177,11 @@ pub(crate) fn erfcx_cody(x: f64) -> f64 { return result; } else if y >= XHUGE { result = SQRPI / y; - if x < ZERO { + if x.is_sign_negative() { if x < XNEG { result = XINF; } else { - let ysq = (x * SIXTEN).trunc() / SIXTEN; + let ysq = (x * 16.0).trunc() / 16.0; let del = (x - ysq) * (x + ysq); let y = (ysq * ysq).exp() * del.exp(); result = (y + y) - result; @@ -196,7 +190,7 @@ pub(crate) fn erfcx_cody(x: f64) -> f64 { return result; } } else { - ysq = ONE / (y * y); + ysq = (y * y).recip(); xnum = P[5] * ysq; xden = ysq; @@ -207,11 +201,11 @@ pub(crate) fn erfcx_cody(x: f64) -> f64 { result = ysq * (xnum + P[4]) / (xden + Q[4]); result = (SQRPI - result) / y; } - if x < ZERO { + if x.is_sign_negative() { if x < XNEG { result = XINF; } else { - let ysq = (x * SIXTEN).trunc() / SIXTEN; + let ysq = (x * 16.0).trunc() / 16.0; let del = (x - ysq) * (x + ysq); let y = (ysq * ysq).exp() * del.exp(); result = (y + y) - result; @@ -223,7 +217,7 @@ pub(crate) fn erfcx_cody(x: f64) -> f64 { #[cfg(test)] mod tests { use crate::erf_cody::{ - erfc_cody, erfcx_cody, FOUR, THRESH, XBIG, XHUGE, XMAX, XNEG, ZERO, + erfc_cody, erfcx_cody, THRESH, XBIG, XHUGE, XMAX, XNEG, }; #[test] @@ -237,13 +231,13 @@ mod tests { let x = erfc_cody(-THRESH + f64::EPSILON); assert_eq!(x, 1.4926134732179377); - let x = erfc_cody(FOUR + f64::EPSILON); + let x = erfc_cody(4.0 + f64::EPSILON); assert_eq!(x, 1.541725790028002e-8); - let x = erfc_cody(FOUR - f64::EPSILON); + let x = erfc_cody(4.0 - f64::EPSILON); assert_eq!(x, 1.541725790028002e-8); - let x = erfc_cody(-FOUR - f64::EPSILON); + let x = erfc_cody(-4.0 - f64::EPSILON); assert_eq!(x, 1.999999984582742); - let x = erfc_cody(-FOUR + f64::EPSILON); + let x = erfc_cody(-4.0 + f64::EPSILON); assert_eq!(x, 1.999999984582742); let x = erfc_cody(XBIG + f64::EPSILON); @@ -273,9 +267,9 @@ mod tests { let x = erfc_cody(-XHUGE + f64::EPSILON); assert_eq!(x, 2.0); - let x = erfc_cody(ZERO + f64::EPSILON); + let x = erfc_cody(0.0 + f64::EPSILON); assert_eq!(x, 0.9999999999999998); - let x = erfc_cody(ZERO - f64::EPSILON); + let x = erfc_cody(0.0 - f64::EPSILON); assert_eq!(x, 1.0000000000000002); let x = erfc_cody(XNEG + f64::EPSILON); @@ -295,13 +289,13 @@ mod tests { let x = erfcx_cody(-THRESH + f64::EPSILON); assert_eq!(x, 1.8594024168714214); - let x = erfcx_cody(FOUR + f64::EPSILON); + let x = erfcx_cody(4.0 + f64::EPSILON); assert_eq!(x, 0.1369994576250614); - let x = erfcx_cody(FOUR - f64::EPSILON); + let x = erfcx_cody(4.0 - f64::EPSILON); assert_eq!(x, 0.1369994576250614); - let x = erfcx_cody(-FOUR - f64::EPSILON); + let x = erfcx_cody(-4.0 - f64::EPSILON); assert_eq!(x, 17772220.904016286); - let x = erfcx_cody(-FOUR + f64::EPSILON); + let x = erfcx_cody(-4.0 + f64::EPSILON); assert_eq!(x, 17772220.904016286); let x = erfcx_cody(XBIG + f64::EPSILON); @@ -331,9 +325,9 @@ mod tests { let x = erfcx_cody(-XHUGE + f64::EPSILON); assert_eq!(x, 1.7976931348623157e308); - let x = erfcx_cody(ZERO + f64::EPSILON); + let x = erfcx_cody(0.0 + f64::EPSILON); assert_eq!(x, 0.9999999999999998); - let x = erfcx_cody(ZERO - f64::EPSILON); + let x = erfcx_cody(0.0 - f64::EPSILON); assert_eq!(x, 1.0000000000000002); let x = erfcx_cody(XNEG + f64::EPSILON); diff --git a/src/lets_be_rational.rs b/src/lets_be_rational.rs index cb20689..ddad641 100644 --- a/src/lets_be_rational.rs +++ b/src/lets_be_rational.rs @@ -1,5 +1,5 @@ use std::f64::consts::SQRT_2; -use crate::constants::{DENORMALISATION_CUTOFF, FOURTH_ROOT_DBL_EPSILON, LN_TWO_PI, PI_OVER_SIX, SIXTEENTH_ROOT_DBL_EPSILON, SQRT_DBL_MAX, SQRT_DBL_MIN, SQRT_ONE_OVER_THREE, SQRT_PI_OVER_TWO, SQRT_THREE, SQRT_THREE_OVER_THIRD_ROOT_TWO_PI, SQRT_TWO_PI, TWO_PI, TWO_PI_OVER_SQRT_TWENTY_SEVEN, VOLATILITY_VALUE_TO_SIGNAL_PRICE_IS_ABOVE_MAXIMUM, VOLATILITY_VALUE_TO_SIGNAL_PRICE_IS_BELOW_INTRINSIC}; +use crate::constants::{DENORMALISATION_CUTOFF, FOURTH_ROOT_DBL_EPSILON, HALF_OF_LN_TWO_PI, SIXTEENTH_ROOT_DBL_EPSILON, SQRT_DBL_MAX, SQRT_MIN_POSITIVE, SQRT_ONE_OVER_THREE, SQRT_PI_OVER_TWO, SQRT_THREE, SQRT_THREE_OVER_THIRD_ROOT_TWO_PI, SQRT_TWO_PI, TWO_PI_OVER_SQRT_TWENTY_SEVEN, VOLATILITY_VALUE_TO_SIGNAL_PRICE_IS_ABOVE_MAXIMUM, VOLATILITY_VALUE_TO_SIGNAL_PRICE_IS_BELOW_INTRINSIC}; use crate::erf_cody::{erfc_cody, erfcx_cody}; use crate::normal_distribution::{inverse_norm_cdf, norm_cdf, norm_pdf}; use crate::rational_cubic::{convex_rational_cubic_control_parameter_to_fit_second_derivative_at_left_side, convex_rational_cubic_control_parameter_to_fit_second_derivative_at_right_side, rational_cubic_interpolation}; @@ -12,16 +12,26 @@ fn householder3_factor(v: f64, h2: f64, h3: f64) -> f64 { (1.0 + 0.5 * h2 * v) / fn householder4_factor(v: f64, h2: f64, h3: f64, h4: f64) -> f64 { (1.0 + v * (h2 + v * h3 / 6.0)) / (1.0 + v * (1.5 * h2 + v * (h2 * h2 / 4.0 + h3 / 3.0 + v * h4 / 24.0))) } fn normalised_intrinsic(x: f64, q: bool) -> f64 { - if (q && (x <= 0.0)) || (!q && (x >= 0.0)) { + if (q && !x.is_sign_positive()) || (!q && !x.is_sign_negative()) { return 0.0; } let x2 = x * x; if x2 < 98.0 * FOURTH_ROOT_DBL_EPSILON { - return (if q { 1.0 } else { -1.0 }) * x * (1.0 + x2 * ((1.0 / 24.0) + x2 * ((1.0 / 1920.0) + x2 * ((1.0 / 322560.0) + (1.0 / 92897280.0) * x2)))).abs().max(0.0); + let ret = x * (1.0 + x2 * ((1.0 / 24.0) + x2 * ((1.0 / 1920.0) + x2 * ((1.0 / 322560.0) + (1.0 / 92897280.0) * x2)))).abs().max(0.0); + return if q { + ret + } else { + -ret + }; } let b_max = (0.5 * x).exp(); let one_over_b_max = b_max.recip(); - (if q { 1.0 } else { -1.0 }) * (b_max - one_over_b_max).abs().max(0.0) + let ret = (b_max - one_over_b_max).abs().max(0.0); + if q { + ret + } else { + -ret + } } #[inline] @@ -38,7 +48,7 @@ const ASYMPTOTIC_EXPANSION_ACCURACY_THRESHOLD: f64 = -10.0; const SMALL_T_EXPANSION_OF_NORMALISED_BLACK_THRESHOLD: f64 = 2.0 * SIXTEENTH_ROOT_DBL_EPSILON; fn asymptotic_expansion_of_normalised_black_call_over_vega(h: f64, t: f64) -> f64 { - assert!((h < -f64::abs(ASYMPTOTIC_EXPANSION_ACCURACY_THRESHOLD)) && (h + t < -f64::abs(SMALL_T_EXPANSION_OF_NORMALISED_BLACK_THRESHOLD + ASYMPTOTIC_EXPANSION_ACCURACY_THRESHOLD))); + assert!((h < -ASYMPTOTIC_EXPANSION_ACCURACY_THRESHOLD.abs()) && (h + t < -(SMALL_T_EXPANSION_OF_NORMALISED_BLACK_THRESHOLD + ASYMPTOTIC_EXPANSION_ACCURACY_THRESHOLD).abs())); let e = square(t / h); let r = (h + t) * (h - t); let q = square(h / r); @@ -46,35 +56,35 @@ fn asymptotic_expansion_of_normalised_black_call_over_vega(h: f64, t: f64) -> f6 let asymptotic_expansion_sum = 2.0 + q * (-6.0E0 - 2.0 * e + 3.0 * q * (1.0E1 + e * (2.0E1 + 2.0 * e) + 5.0 * q * (-1.4E1 + e * (-7.0E1 + e * (-4.2E1 - 2.0 * e)) + 7.0 * q * (1.8E1 + e * (1.68E2 + e * (2.52E2 + e * (7.2E1 + 2.0 * e))) + 9.0 * q * (-2.2E1 + e * (-3.3E2 + e * (-9.24E2 + e * (-6.6E2 + e * (-1.1E2 - 2.0 * e)))) + 1.1E1 * q * (2.6E1 + e * (5.72E2 + e * (2.574E3 + e * (3.432E3 + e * (1.43E3 + e * (1.56E2 + 2.0 * e))))) + 1.3E1 * q * (-3.0E1 + e * (-9.1E2 + e * (-6.006E3 + e * (-1.287E4 + e * (-1.001E4 + e * (-2.73E3 + e * (-2.1E2 - 2.0 * e)))))) + 1.5E1 * q * (3.4E1 + e * (1.36E3 + e * (1.2376E4 + e * (3.8896E4 + e * (4.862E4 + e * (2.4752E4 + e * (4.76E3 + e * (2.72E2 + 2.0 * e))))))) + 1.7E1 * q * (-3.8E1 + e * (-1.938E3 + e * (-2.3256E4 + e * (-1.00776E5 + e * (-1.84756E5 + e * (-1.51164E5 + e * (-5.4264E4 + e * (-7.752E3 + e * (-3.42E2 - 2.0 * e)))))))) + 1.9E1 * q * (4.2E1 + e * (2.66E3 + e * (4.0698E4 + e * (2.3256E5 + e * (5.8786E5 + e * (7.05432E5 + e * (4.0698E5 + e * (1.08528E5 + e * (1.197E4 + e * (4.2E2 + 2.0 * e))))))))) + 2.1E1 * q * (-4.6E1 + e * (-3.542E3 + e * (-6.7298E4 + e * (-4.90314E5 + e * (-1.63438E6 + e * (-2.704156E6 + e * (-2.288132E6 + e * (-9.80628E5 + e * (-2.01894E5 + e * (-1.771E4 + e * (-5.06E2 - 2.0 * e)))))))))) + 2.3E1 * q * (5.0E1 + e * (4.6E3 + e * (1.0626E5 + e * (9.614E5 + e * (4.08595E6 + e * (8.9148E6 + e * (1.04006E7 + e * (6.53752E6 + e * (2.16315E6 + e * (3.542E5 + e * (2.53E4 + e * (6.0E2 + 2.0 * e))))))))))) + 2.5E1 * q * (-5.4E1 + e * (-5.85E3 + e * (-1.6146E5 + e * (-1.77606E6 + e * (-9.37365E6 + e * (-2.607579E7 + e * (-4.01166E7 + e * (-3.476772E7 + e * (-1.687257E7 + e * (-4.44015E6 + e * (-5.9202E5 + e * (-3.51E4 + e * (-7.02E2 - 2.0 * e)))))))))))) + 2.7E1 * q * (5.8E1 + e * (7.308E3 + e * (2.3751E5 + e * (3.12156E6 + e * (2.003001E7 + e * (6.919458E7 + e * (1.3572783E8 + e * (1.5511752E8 + e * (1.0379187E8 + e * (4.006002E7 + e * (8.58429E6 + e * (9.5004E5 + e * (4.7502E4 + e * (8.12E2 + 2.0 * e))))))))))))) + 2.9E1 * q * (-6.2E1 + e * (-8.99E3 + e * (-3.39822E5 + e * (-5.25915E6 + e * (-4.032015E7 + e * (-1.6934463E8 + e * (-4.1250615E8 + e * (-6.0108039E8 + e * (-5.3036505E8 + e * (-2.8224105E8 + e * (-8.870433E7 + e * (-1.577745E7 + e * (-1.472562E6 + e * (-6.293E4 + e * (-9.3E2 - 2.0 * e)))))))))))))) + 3.1E1 * q * (6.6E1 + e * (1.0912E4 + e * (4.74672E5 + e * (8.544096E6 + e * (7.71342E7 + e * (3.8707344E8 + e * (1.14633288E9 + e * (2.07431664E9 + e * (2.33360622E9 + e * (1.6376184E9 + e * (7.0963464E8 + e * (1.8512208E8 + e * (2.7768312E7 + e * (2.215136E6 + e * (8.184E4 + e * (1.056E3 + 2.0 * e))))))))))))))) + 3.3E1 * (-7.0E1 + e * (-1.309E4 + e * (-6.49264E5 + e * (-1.344904E7 + e * (-1.4121492E8 + e * (-8.344518E8 + e * (-2.9526756E9 + e * (-6.49588632E9 + e * (-9.0751353E9 + e * (-8.1198579E9 + e * (-4.6399188E9 + e * (-1.6689036E9 + e * (-3.67158792E8 + e * (-4.707164E7 + e * (-3.24632E6 + e * (-1.0472E5 + e * (-1.19E3 - 2.0 * e))))))))))))))))) * q)))))))))))))))); let b_over_vega = (t / r) * asymptotic_expansion_sum; - f64::abs(f64::max(b_over_vega, 0.)) + b_over_vega.max(0.0).abs() } fn small_t_expansion_of_normalised_black_call_over_vega(h: f64, t: f64) -> f64 { let w = t * t; let h2 = h * h; - let a = 1_f64 + h * SQRT_PI_OVER_TWO * erfcx_cody(-(1.0 / SQRT_2) * h); + let a = 1_f64 + h * SQRT_PI_OVER_TWO * erfcx_cody(-SQRT_2.recip() * h); let b_over_vega = 2.0 * t * (a + w * ((-1.0 + 3.0 * a + a * h2) / 6.0 + w * ((-7.0 + 15.0 * a + h2 * (-1.0 + 10.0 * a + a * h2)) / 120.0 + w * ((-57.0 + 105.0 * a + h2 * (-18.0 + 105.0 * a + h2 * (-1.0 + 21.0 * a + a * h2))) / 5040.0 + w * ((-561.0 + 945.0 * a + h2 * (-285.0 + 1260.0 * a + h2 * (-33.0 + 378.0 * a + h2 * (-1.0 + 36.0 * a + a * h2)))) / 362880.0 + w * ((-6555.0 + 10395.0 * a + h2 * (-4680.0 + 17325.0 * a + h2 * (-840.0 + 6930.0 * a + h2 * (-52.0 + 990.0 * a + h2 * (-1.0 + 55.0 * a + a * h2))))) / 39916800.0 + ((-89055.0 + 135135.0 * a + h2 * (-82845.0 + 270270.0 * a + h2 * (-20370.0 + 135135.0 * a + h2 * (-1926.0 + 25740.0 * a + h2 * (-75.0 + 2145.0 * a + h2 * (-1.0 + 78.0 * a + a * h2)))))) * w) / 6227020800.0)))))); - f64::abs(f64::max(b_over_vega, 0_f64)) + b_over_vega.max(0.0).abs() } fn normalised_black_call_with_optimal_use_of_codys_functions(x: f64, s: f64) -> f64 { const CODYS_THRESHOLD: f64 = 0.46875; let h = x / s; let t = 0.5 * s; - let q1 = -(1f64 / SQRT_2) * (h + t); - let q2 = -(1f64 / SQRT_2) * (h - t); - let two_b: f64; - if q1 < CODYS_THRESHOLD { - if q2 < CODYS_THRESHOLD { - two_b = 0.5 * ((0.5 * x).exp() * erfc_cody(q1) - (-0.5 * x).exp() * erfc_cody(q2)); + let q1 = -SQRT_2.recip() * (h + t); + let q2 = -SQRT_2.recip() * (h - t); + let two_b: f64 = + if q1 < CODYS_THRESHOLD { + if q2 < CODYS_THRESHOLD { + 0.5 * ((0.5 * x).exp() * erfc_cody(q1) - (-0.5 * x).exp() * erfc_cody(q2)) + } else { + 0.5 * ((0.5 * x).exp() * erfc_cody(q1) - (-0.5 * (h * h + t * t)).exp() * erfcx_cody(q2)) + } + } else if q2 < CODYS_THRESHOLD { + 0.5 * ((-0.5 * (h * h + t * t)).exp() * erfcx_cody(q1) - (-0.5 * x).exp() * erfc_cody(q2)) } else { - two_b = 0.5 * ((0.5 * x).exp() * erfc_cody(q1) - (-0.5 * (h * h + t * t)).exp() * erfcx_cody(q2)); - } - } else if q2 < CODYS_THRESHOLD { - two_b = 0.5 * ((-0.5 * (h * h + t * t)).exp() * erfcx_cody(q1) - (-0.5 * x).exp() * erfc_cody(q2)); - } else { - two_b = 0.5 * ((-0.5 * (h * h + t * t)).exp() * (erfcx_cody(q1) - erfcx_cody(q2))); - } + 0.5 * ((-0.5 * (h * h + t * t)).exp() * (erfcx_cody(q1) - erfcx_cody(q2))) + }; two_b.abs().max(0.0) } @@ -82,7 +92,7 @@ fn normalised_vega(x: f64, s: f64) -> f64 { let ax = x.abs(); if ax <= 0.0 { (1.0 / SQRT_TWO_PI) * (-0.125 * s * s).exp() - } else if s <= 0.0 || s <= ax * SQRT_DBL_MIN { + } else if s <= 0.0 || s <= ax * SQRT_MIN_POSITIVE { 0.0 } else { (1.0 / SQRT_TWO_PI) * (-0.5 * (square(x / s) + square(0.5 * s))).exp() @@ -92,16 +102,16 @@ fn normalised_vega(x: f64, s: f64) -> f64 { fn ln_normalised_vega(x: f64, s: f64) -> f64 { let ax = x.abs(); if ax <= 0.0 { - -(LN_TWO_PI / 2.0) - 0.125 * s * s + -HALF_OF_LN_TWO_PI - 0.125 * s * s } else if s <= 0.0 { f64::MIN } else { - -(LN_TWO_PI / 2.0) - 0.5 * (square(x / s) + square(0.5 * s)) + -HALF_OF_LN_TWO_PI - 0.5 * (square(x / s) + 0.25 * s * s) } } fn normalised_black_call(x: f64, s: f64) -> f64 { - if x > 0.0 { + if x.is_sign_positive() { return normalised_intrinsic_call(x) + normalised_black_call(-x, s); } if s <= x.abs() * DENORMALISATION_CUTOFF { @@ -118,7 +128,7 @@ fn normalised_black_call(x: f64, s: f64) -> f64 { fn normalised_black_call_over_vega_and_ln_vega(x: f64, s: f64) -> (f64, f64) { - if x > 0.0 { + if x.is_sign_positive() { let (bx, ln_vega) = normalised_black_call_over_vega_and_ln_vega(-x, s); return (normalised_intrinsic_call(x) * (-ln_vega).exp() + bx, ln_vega); } @@ -142,7 +152,7 @@ fn normalised_black(x: f64, s: f64, theta: bool) -> f64 { pub(crate) fn black(f: f64, k: f64, sigma: f64, t: f64, q: bool) -> f64 { let intrinsic = if !q { k - f } else { f - k }.max(0f64).abs(); - if (q && ((f - k) > 0.0)) || (!q && ((f - k) < 0.0)) { + if (q && ((f - k).is_sign_positive())) || (!q && ((f - k).is_sign_negative())) { return intrinsic + black(f, k, sigma, t, !q); } intrinsic.max((f.sqrt() * k.sqrt()) * normalised_black((f / k).ln(), sigma * t.sqrt(), q)) @@ -155,12 +165,12 @@ fn compute_f_lower_map_and_first_two_derivatives(x: f64, s: f64) -> (f64, f64, f let s2 = s * s; let phi_m = norm_cdf(-z); let phi = norm_pdf(z); - let fpp = PI_OVER_SIX * y / (s2 * s) * phi_m * (8.0 * SQRT_THREE * s * ax + (3.0 * s2 * (s2 - 8.0) - 8.0 * x * x) * phi_m / phi) * (2.0 * y + 0.25 * s2).exp(); + let fpp = std::f64::consts::FRAC_PI_6 * y / (s2 * s) * phi_m * (8.0 * SQRT_THREE * s * ax + (3.0 * s2 * (s2 - 8.0) - 8.0 * x * x) * phi_m / phi) * (2.0 * y + 0.25 * s2).exp(); let (fp, f) = if s.is_subnormal() { (1.0, 0.0) } else { let phi2 = phi_m * phi_m; - let fp_val = TWO_PI * y * phi2 * (y + 0.125 * s * s).exp(); + let fp_val = std::f64::consts::TAU * y * phi2 * (y + 0.125 * s * s).exp(); let f_val = if x.is_subnormal() { 0.0 } else { @@ -207,15 +217,14 @@ fn take_step(x_min: f64, x_max: f64, x: f64, dx: f64) -> (f64, f64) { } fn unchecked_normalised_implied_volatility_from_a_transformed_rational_guess_with_limited_iterations( - mut beta: f64, mut x: f64, q: bool, n: i32, + mut beta: f64, mut x: f64, q: bool, n: u8, ) -> f64 { - if (q && (x > 0.0)) || (!q && (x < 0.0)) { + if (q && (x.is_sign_positive())) || (!q && (x.is_sign_negative())) { beta = (beta - normalised_intrinsic(x, q)).max(0.).abs(); } if !q { x = -x; } - assert!(x <= 0.); if beta <= 0. || beta < DENORMALISATION_CUTOFF { return 0.0; } @@ -227,9 +236,9 @@ fn unchecked_normalised_implied_volatility_from_a_transformed_rational_guess_wit let mut f = f64::MIN; let mut s; let mut ds = f64::MIN; - let mut s_left = f64::MIN; + let mut s_left = f64::MIN_POSITIVE; let mut s_right = f64::MAX; - let s_c = (2. * x.abs()).sqrt(); + let s_c = SQRT_2 * x.abs().sqrt(); let b_c = normalised_black_call(x, s_c); let v_c = normalised_vega(x, s_c); if beta < b_c { @@ -274,19 +283,19 @@ fn unchecked_normalised_implied_volatility_from_a_transformed_rational_guess_wit return s; } else { let v1 = normalised_vega(x, s1); - let r_im = convex_rational_cubic_control_parameter_to_fit_second_derivative_at_right_side(b1, b_c, s1, s_c, 1.0 / v1, 1.0 / v_c, 0.0, false); - s = rational_cubic_interpolation(beta, b1, b_c, s1, s_c, 1.0 / v1, 1.0 / v_c, r_im); + let r_im = convex_rational_cubic_control_parameter_to_fit_second_derivative_at_right_side(b1, b_c, s1, s_c, v1.recip(), v_c.recip(), 0.0, false); + s = rational_cubic_interpolation(beta, b1, b_c, s1, s_c, v1.recip(), v_c.recip(), r_im); s_left = s1; s_right = s_c; } } else { - let s_u = if v_c > f64::MIN { s_c + (b_max - b_c) / v_c } else { s_c }; + let s_u = if v_c > f64::MIN_POSITIVE { s_c + (b_max - b_c) / v_c } else { s_c }; let b_u = normalised_black_call(x, s_u); if beta <= b_u { let v_u = normalised_vega(x, s_u); let r_u_m = convex_rational_cubic_control_parameter_to_fit_second_derivative_at_left_side( - b_c, b_u, s_c, s_u, 1.0 / v_c, 1.0 / v_u, 0.0, false); - s = rational_cubic_interpolation(beta, b_c, b_u, s_c, s_u, 1.0 / v_c, 1.0 / v_u, r_u_m); + b_c, b_u, s_c, s_u, v_c.recip(), v_u.recip(), 0.0, false); + s = rational_cubic_interpolation(beta, b_c, b_u, s_c, s_u, v_c.recip(), v_u.recip(), r_u_m); s_left = s_c; s_right = s_u; } else { @@ -296,7 +305,7 @@ fn unchecked_normalised_implied_volatility_from_a_transformed_rational_guess_wit let r_uu = convex_rational_cubic_control_parameter_to_fit_second_derivative_at_left_side(b_u, b_max, f_upper_map_h, 0.0, d_f_upper_map_h_d_beta, -0.5, d2_f_upper_map_h_d_beta2, true); f = rational_cubic_interpolation(beta, b_u, b_max, f_upper_map_h, 0.0, d_f_upper_map_h_d_beta, -0.5, r_uu); } - if f <= 0.0 { + if !f.is_sign_positive() { let h = b_max - b_u; let t = (beta - b_u) / h; f = (f_upper_map_h * (1.0 - t) + 0.5 * h * t) * (1.0 - t); @@ -307,7 +316,7 @@ fn unchecked_normalised_implied_volatility_from_a_transformed_rational_guess_wit while iterations < n && ds.abs() > f64::EPSILON * s { let h = x / s; let t = s / 2.0; - let gp = (2.0 / SQRT_TWO_PI) / (erfcx_cody((t + h) * (1.0 / SQRT_2)) + erfcx_cody((t - h) * (1.0 / SQRT_2))); + let gp = (2.0 / SQRT_TWO_PI) / (erfcx_cody((t + h) * SQRT_2.recip()) + erfcx_cody((t - h) * SQRT_2.recip())); let b_bar = normalised_vega(x, s) / gp; let g = (beta_bar / b_bar).ln(); let x_over_s_square = (h * h) / s; @@ -353,7 +362,7 @@ fn implied_volatility_from_a_transformed_rational_guess_with_limited_iterations( k: f64, t: f64, mut q: bool, - n: i32, + n: u8, ) -> f64 { let intrinsic = (if !q { k - f } else { f - k }).max(0.0).abs(); if price < intrinsic { @@ -367,7 +376,7 @@ fn implied_volatility_from_a_transformed_rational_guess_with_limited_iterations( } let x = (f / k).ln(); // Map in-the-money to out-of-the-money - if (q && (x > 0.0)) || (!q && (x < 0.0)) { // q && x > 0.0 or !q && x < 0.0 + if (q && (x.is_sign_positive())) || (!q && (x.is_sign_negative())) { price = (price - intrinsic).max(0.0).abs(); q = !q; } diff --git a/src/normal_distribution.rs b/src/normal_distribution.rs index 6e3299d..81e0b66 100644 --- a/src/normal_distribution.rs +++ b/src/normal_distribution.rs @@ -115,7 +115,7 @@ pub(crate) fn inverse_norm_cdf(u: f64) -> f64 { q * (((((((A7 * r + A6) * r + A5) * r + A4) * r + A3) * r + A2) * r + A1) * r + A0) / (((((((B7 * r + B6) * r + B5) * r + B4) * r + B3) * r + B2) * r + B1) * r + 1.0) } else { - let mut r = if q < 0.0 { u } else { 1.0 - u }; + let mut r = if q.is_sign_negative() { u } else { 1.0 - u }; r = (-r.ln()).sqrt(); let ret = if r < SPLIT2 { @@ -127,6 +127,6 @@ pub(crate) fn inverse_norm_cdf(u: f64) -> f64 { (((((((E7 * r + E6) * r + E5) * r + E4) * r + E3) * r + E2) * r + E1) * r + E0) / (((((((F7 * r + F6) * r + F5) * r + F4) * r + F3) * r + F2) * r + F1) * r + 1.0) }; - if q < 0.0 { -ret } else { ret } + if q.is_sign_negative() { -ret } else { ret } } } \ No newline at end of file diff --git a/src/rational_cubic.rs b/src/rational_cubic.rs index 15e1f8a..34160d5 100644 --- a/src/rational_cubic.rs +++ b/src/rational_cubic.rs @@ -51,7 +51,7 @@ pub(crate) fn rational_cubic_control_parameter_to_fit_second_derivative_at_left_ } let denominator = (y_r - y_l) / h - d_l; if is_zero(denominator) { - if numerator > 0.0 { + if numerator.is_sign_positive() { MAXIMUM_RATIONAL_CUBIC_CONTROL_PARAMETER_VALUE } else { MINIMUM_RATIONAL_CUBIC_CONTROL_PARAMETER_VALUE @@ -77,7 +77,7 @@ pub(crate) fn rational_cubic_control_parameter_to_fit_second_derivative_at_right } let denominator = d_r - (y_r - y_l) / h; if is_zero(denominator) { - return if numerator > 0.0 { + return if numerator.is_sign_positive() { MAXIMUM_RATIONAL_CUBIC_CONTROL_PARAMETER_VALUE } else { MINIMUM_RATIONAL_CUBIC_CONTROL_PARAMETER_VALUE