Refactor normal distribution for f64 constants

CONTRIBUTING.md

@@ -0,0 +1,8 @@
+# Contributing
+
+The tests in this project require unstable Rust. To switch to unstable in
+this repository, run e.g.
+
+```shell
+rustup override set nightly
+```

Cargo.toml

@@ -23,3 +23,7 @@ lto = true
 
 [badges]
 travis-ci = { repository = "danielhstahl/black_scholes_rust" }
+
+[[example]]
+name = "profile"
+src = "examples/profile.rs"

benches/benches.rs

@@ -1,14 +1,18 @@
 #![feature(test)]
 extern crate test;
-use test::Bencher;
+use test::{black_box, Bencher};
+
+/// Convenience alias for [`black_box`].
+const BB: fn(f64) -> f64 = black_box;
+
 #[bench]
 fn bench_call_price(b: &mut Bencher) {
     let r = 0.05;
     let sig = 0.3;
     let t = 1.0;
     let asset = 50.0;
     let k = 50.0;
-    b.iter(|| black_scholes::call(asset, k, r, sig, t))
+    b.iter(|| black_scholes::call(BB(asset), BB(k), BB(r), BB(sig), BB(t)))
 }
 
 #[bench]
@@ -18,7 +22,7 @@ fn bench_all_price(b: &mut Bencher) {
     let t = 1.0;
     let asset = 50.0;
     let k = 50.0;
-    b.iter(|| black_scholes::compute_all(asset, k, r, sig, t))
+    b.iter(|| black_scholes::compute_all(BB(asset), BB(k), BB(r), BB(sig), BB(t)))
 }
 
 #[bench]
@@ -29,18 +33,18 @@ fn bench_all_price_no_cache(b: &mut Bencher) {
     let asset = 50.0;
     let k = 50.0;
     b.iter(|| {
-        black_scholes::call(asset, k, r, sig, t);
-        black_scholes::call_delta(asset, k, r, sig, t);
-        black_scholes::call_gamma(asset, k, r, sig, t);
-        black_scholes::call_vega(asset, k, r, sig, t);
-        black_scholes::call_theta(asset, k, r, sig, t);
-        black_scholes::call_rho(asset, k, r, sig, t);
+        black_scholes::call(BB(asset), BB(k), BB(r), BB(sig), BB(t));
+        black_scholes::call_delta(BB(asset), BB(k), BB(r), BB(sig), BB(t));
+        black_scholes::call_gamma(BB(asset), BB(k), BB(r), BB(sig), BB(t));
+        black_scholes::call_vega(BB(asset), BB(k), BB(r), BB(sig), BB(t));
+        black_scholes::call_theta(BB(asset), BB(k), BB(r), BB(sig), BB(t));
+        black_scholes::call_rho(BB(asset), BB(k), BB(r), BB(sig), BB(t));
 
-        black_scholes::put(asset, k, r, sig, t);
-        black_scholes::put_delta(asset, k, r, sig, t);
-        black_scholes::put_gamma(asset, k, r, sig, t);
-        black_scholes::put_vega(asset, k, r, sig, t);
-        black_scholes::put_theta(asset, k, r, sig, t);
-        black_scholes::put_rho(asset, k, r, sig, t);
+        black_scholes::put(BB(asset), BB(k), BB(r), BB(sig), BB(t));
+        black_scholes::put_delta(BB(asset), BB(k), BB(r), BB(sig), BB(t));
+        black_scholes::put_gamma(BB(asset), BB(k), BB(r), BB(sig), BB(t));
+        black_scholes::put_vega(BB(asset), BB(k), BB(r), BB(sig), BB(t));
+        black_scholes::put_theta(BB(asset), BB(k), BB(r), BB(sig), BB(t));
+        black_scholes::put_rho(BB(asset), BB(k), BB(r), BB(sig), BB(t));
     })
 }

examples/profile.rs

@@ -0,0 +1,28 @@
+use std::hint::black_box;
+
+/// Convenience alias for [`black_box`].
+const BB: fn(f64) -> f64 = black_box;
+
+fn main() {
+    let r = 0.05;
+    let sig = 0.3;
+    let t = 1.0;
+    let asset = 50.0;
+    let k = 50.0;
+
+    for _ in 0..100_000 {
+        black_scholes::call(BB(asset), BB(k), BB(r), BB(sig), BB(t));
+        black_scholes::call_delta(BB(asset), BB(k), BB(r), BB(sig), BB(t));
+        black_scholes::call_gamma(BB(asset), BB(k), BB(r), BB(sig), BB(t));
+        black_scholes::call_vega(BB(asset), BB(k), BB(r), BB(sig), BB(t));
+        black_scholes::call_theta(BB(asset), BB(k), BB(r), BB(sig), BB(t));
+        black_scholes::call_rho(BB(asset), BB(k), BB(r), BB(sig), BB(t));
+
+        black_scholes::put(BB(asset), BB(k), BB(r), BB(sig), BB(t));
+        black_scholes::put_delta(BB(asset), BB(k), BB(r), BB(sig), BB(t));
+        black_scholes::put_gamma(BB(asset), BB(k), BB(r), BB(sig), BB(t));
+        black_scholes::put_vega(BB(asset), BB(k), BB(r), BB(sig), BB(t));
+        black_scholes::put_theta(BB(asset), BB(k), BB(r), BB(sig), BB(t));
+        black_scholes::put_rho(BB(asset), BB(k), BB(r), BB(sig), BB(t));
+    }
+}

src/lib.rs

@@ -2,24 +2,27 @@
 //! A Black Scholes option pricing library.
 use serde::Serialize;
 use special::Error;
-use std::f64::consts::{FRAC_2_SQRT_PI, PI, SQRT_2};
+use std::f64::consts::{FRAC_1_PI, FRAC_1_SQRT_2, FRAC_2_SQRT_PI, SQRT_2};
+
+/// 1/sqrt(2π)
+#[allow(clippy::excessive_precision)]
+const FRAC_1_SQRT_2PI: f64 = 0.3989422804014326779399460599343818684758586311649346576659258296;
 
 fn cum_norm(x: f64) -> f64 {
-    (x / SQRT_2).error() * 0.5 + 0.5
+    (x * FRAC_1_SQRT_2).error() * 0.5 + 0.5
 }
+
 fn inc_norm(x: f64) -> f64 {
-    (-x.powi(2) / 2.0).exp() / (PI.sqrt() * SQRT_2)
+    (-x.powi(2) * 0.5).exp() * FRAC_1_SQRT_2PI
 }
 
 fn d1(s: f64, k: f64, discount: f64, sqrt_maturity_sigma: f64) -> f64 {
     (s / (k * discount)).ln() / sqrt_maturity_sigma + 0.5 * sqrt_maturity_sigma
 }
+
+#[inline(always)]
 fn max_or_zero(v: f64) -> f64 {
-    if v > 0.0 {
-        v
-    } else {
-        0.0
-    }
+    v.max(0.0)
 }
 
 /// Returns BS call option formula with discount and volatility already computed.
@@ -490,7 +493,7 @@ fn approximate_vol(price: f64, s: f64, k: f64, rate: f64, maturity: f64) -> f64
     let helper_1 = s - x;
     let c1 = price - helper_1 * 0.5;
     let c2 = c1.powi(2);
-    let c3 = helper_1.powi(2) / PI;
+    let c3 = helper_1.powi(2) * FRAC_1_PI;
     let bridge_1 = c2 - c3;
     let bridge_m = if bridge_1 > 0.0 { bridge_1.sqrt() } else { 0.0 };
     coef * (c1 + bridge_m) / maturity.sqrt()
@@ -731,6 +734,8 @@ mod tests {
     use rand::distributions::{Distribution, Uniform};
     use rand::rngs::StdRng;
     use rand::SeedableRng;
+    use std::f64::consts::PI;
+
     fn get_rng_seed(seed: [u8; 32]) -> StdRng {
         SeedableRng::from_seed(seed)
     }
@@ -1038,4 +1043,9 @@ mod tests {
         assert_approx_eq!(result.put_vega, put_vega(s, k, rate, sigma, maturity));
         assert_approx_eq!(result.put_rho, put_rho(s, k, rate, sigma, maturity));
     }
+
+    #[test]
+    fn constants_are_correct() {
+        assert_approx_eq!(FRAC_1_SQRT_2PI, (2.0 * PI).sqrt().recip());
+    }
 }