Add normal distribution functions

Author: timschmidt

README.md

@@ -320,7 +320,8 @@ assert_eq!(value.exact_rational(), Some(Rational::fraction(7, 8).unwrap()));
 ```
 
 `Simple` supports arithmetic, roots, powers, logs, exponentials, trig, inverse trig,
-inverse hyperbolic functions, integers, decimals, fractions, `pi`, and `e`.
+inverse hyperbolic functions, normal distribution helpers (`erf`, `dnorm`, `pnorm`,
+`qnorm`), integers, decimals, fractions, `pi`, and `e`.
 
 ## Conversions
 

src/computable/approximation.rs

@@ -18,3 +18,4 @@ include!("approximation/trig.rs");
 include!("approximation/logarithms.rs");
 include!("approximation/inverse_trig.rs");
 include!("approximation/inverse_hyperbolic.rs");
+include!("approximation/statistics.rs");

src/computable/approximation/dispatch.rs

@@ -79,6 +79,12 @@ impl Approximation {
             PrescaledTan(c) => tan(signal, c, p),
             PrescaledTanRational(r) => tan_rational(signal, r, p),
             PrescaledCot(c) => cot(signal, c, p),
+            ErfSeries(c) => erf_series(signal, c, p),
+            NormalQuantile {
+                p: prob,
+                seed,
+                seed_prec,
+            } => normal_quantile(signal, prob, seed, *seed_prec, p),
         }
     }
 }
@@ -122,4 +128,3 @@ fn raw(kind: Approximation) -> Computable {
         signal: None,
     }
 }
-

src/computable/approximation/representation.rs

@@ -143,6 +143,12 @@ pub(super) enum Approximation {
     // the same local quotient kernel once digits are requested.
     PrescaledTanRational(Rational),
     PrescaledCot(Computable),
+    ErfSeries(Computable),
+    NormalQuantile {
+        p: Computable,
+        seed: BigInt,
+        seed_prec: Precision,
+    },
 }
 
 #[derive(Clone, Copy, Debug, PartialEq, Eq)]
@@ -165,4 +171,3 @@ pub(super) enum SharedConstant {
     AtanInv2,
     AtanInv5,
 }
-

src/computable/approximation/statistics.rs

@@ -0,0 +1,75 @@
+fn to_prec(n: &BigInt) -> Precision {
+    n.to_i32().unwrap_or(Precision::MAX)
+}
+
+fn erf_series(signal: &Option<Signal>, op: &Computable, p: Precision) -> BigInt {
+    let rough_x = op.approx_signal(signal, -10);
+    let x_sq_approx = (&rough_x * &rough_x) >> 20;
+
+    let n_estimate = {
+        let estimate = &x_sq_approx + BigInt::from(-p) + BigInt::from(10);
+        if estimate < BigInt::one() {
+            BigInt::one()
+        } else {
+            estimate
+        }
+    };
+    let magnitude_bits = to_prec(&((&x_sq_approx * BigInt::from(3)) / BigInt::from(2))) + 2;
+    let guard_bits = (n_estimate.magnitude().bits() as Precision) + magnitude_bits + 4;
+    let calc_precision = p - guard_bits;
+    let op_prec = calc_precision - 8;
+    let op_appr = op.approx_signal(signal, op_prec);
+    let max_trunc_error = signed::ONE.deref() << (p - 4 - calc_precision);
+
+    let mut n: i64 = 0;
+    let mut current_term = scale(op_appr.clone(), op_prec - calc_precision);
+    let mut current_sum = current_term.clone();
+    loop {
+        if should_stop(signal) {
+            break;
+        }
+        let prev = current_term;
+        let mut t = scale(&prev * &op_appr, op_prec);
+        t = scale(t * &op_appr, op_prec);
+        t = (t * signed::TWO.deref()) / BigInt::from(2 * n + 3);
+        n += 1;
+        if t.is_zero() {
+            break;
+        }
+
+        let prev_abs = prev.abs();
+        let cur_abs = t.abs();
+        if cur_abs < prev_abs {
+            let denom = &prev_abs - &cur_abs;
+            if &cur_abs * &prev_abs < &max_trunc_error * &denom {
+                break;
+            }
+        }
+        current_sum += &t;
+        current_term = t;
+    }
+    scale(current_sum, calc_precision - p)
+}
+
+const NORMAL_QUANTILE_RESULT_MSD_BOUND: Precision = 5;
+
+fn normal_quantile(
+    signal: &Option<Signal>,
+    p: &Computable,
+    seed: &BigInt,
+    seed_prec: Precision,
+    prec: Precision,
+) -> BigInt {
+    if prec >= seed_prec || should_stop(signal) {
+        return scale(seed.clone(), seed_prec - prec);
+    }
+
+    let appr_prec = (NORMAL_QUANTILE_RESULT_MSD_BOUND + prec) / 2 - 6;
+    let xn_int = normal_quantile(signal, p, seed, seed_prec, appr_prec);
+    let xn = Computable::integer(xn_int).shift_left(appr_prec);
+    let fx = xn.clone().pnorm();
+    let phi_xn = xn.clone().dnorm();
+    let x_next = xn.add(fx.add(p.clone().negate()).multiply(phi_xn.inverse()).negate());
+
+    (x_next.approx_signal(signal, prec - 2) + signed::TWO.deref()) >> 2
+}

src/computable/node/algebra.rs

@@ -483,6 +483,55 @@ impl Computable {
         }
     }
 
+    /// Error function, erf(x).
+    pub fn erf(self) -> Computable {
+        if let Some(rational) = self.exact_rational()
+            && rational.sign() == Sign::NoSign
+        {
+            return Self::zero();
+        }
+        let series = Self {
+            internal: Box::new(Approximation::ErfSeries(self.clone())),
+            cache: RefCell::new(Cache::Invalid),
+            bound: RefCell::new(BoundCache::Invalid),
+            exact_sign: RefCell::new(ExactSignCache::Invalid),
+            signal: None,
+        };
+        let gaussian = self.square().negate().exp();
+        let two_over_sqrt_pi = Self::pi().sqrt().inverse().shift_left(1);
+        two_over_sqrt_pi.multiply(gaussian).multiply(series)
+    }
+
+    /// Standard normal CDF.
+    pub fn pnorm(self) -> Computable {
+        if let Some(rational) = self.exact_rational()
+            && rational.sign() == Sign::NoSign
+        {
+            return Self::rational(HALF_RATIONAL.clone());
+        }
+        let sqrt2 = Self::sqrt_constant(2).unwrap_or_else(|| Self::integer(BigInt::from(2)).sqrt());
+        let z = self.multiply(sqrt2.inverse());
+        Self::one().add(z.erf()).shift_right(1)
+    }
+
+    /// Standard normal density.
+    pub fn dnorm(self) -> Computable {
+        let neg_half_x_sq = self.square().shift_right(1).negate();
+        let sqrt_2pi = Self::pi().shift_left(1).sqrt();
+        neg_half_x_sq.exp().multiply(sqrt_2pi.inverse())
+    }
+
+    /// Standard normal quantile by Newton iteration with the analytic density.
+    pub fn normal_quantile(p: Computable, seed: BigInt, seed_prec: Precision) -> Computable {
+        Self {
+            internal: Box::new(Approximation::NormalQuantile { p, seed, seed_prec }),
+            cache: RefCell::new(Cache::Invalid),
+            bound: RefCell::new(BoundCache::Invalid),
+            exact_sign: RefCell::new(ExactSignCache::Invalid),
+            signal: None,
+        }
+    }
+
     /// Attach an abort signal checked by long-running approximation routines.
     pub fn abort(&mut self, s: Signal) {
         self.signal = Some(s);

src/computable/node/primitive_constructors.rs

@@ -165,6 +165,25 @@ impl Computable {
                 (Approximation::PrescaledCot(left), Approximation::PrescaledCot(right)) => {
                     Computable::internal_structural_eq(left, right)
                 }
+                (Approximation::ErfSeries(left), Approximation::ErfSeries(right)) => {
+                    Computable::internal_structural_eq(left, right)
+                }
+                (
+                    Approximation::NormalQuantile {
+                        p: left_p,
+                        seed: left_seed,
+                        seed_prec: left_seed_prec,
+                    },
+                    Approximation::NormalQuantile {
+                        p: right_p,
+                        seed: right_seed,
+                        seed_prec: right_seed_prec,
+                    },
+                ) => {
+                    left_seed == right_seed
+                        && left_seed_prec == right_seed_prec
+                        && Computable::internal_structural_eq(left_p, right_p)
+                }
                 _ => false,
             }
         }
@@ -749,4 +768,3 @@ impl Computable {
         }
     }
 }
-

src/computable/node/tests.rs

@@ -1192,6 +1192,53 @@ mod tests {
         assert_close(value.clone().ln().exp(), value, -40, 2);
     }
 
+    #[test]
+    fn erf_known_values() {
+        assert_close(
+            Computable::zero().erf(),
+            Computable::zero(),
+            -160,
+            2,
+        );
+        assert_close(
+            Computable::rational(Rational::fraction(1, 2).unwrap()).erf(),
+            Computable::rational("0.5204998778130465376827466538919645287364".parse().unwrap()),
+            -90,
+            2,
+        );
+        assert_close(
+            Computable::one().erf(),
+            Computable::rational("0.8427007929497148693412206350826092592960".parse().unwrap()),
+            -90,
+            2,
+        );
+    }
+
+    #[test]
+    fn normal_density_and_cdf_known_values() {
+        assert_close(
+            Computable::zero().dnorm(),
+            Computable::rational("0.39894228040143267793994605993438186847585863".parse().unwrap()),
+            -120,
+            2,
+        );
+        assert_close(
+            Computable::one().pnorm(),
+            Computable::rational("0.8413447460685429485852325456320379224779".parse().unwrap()),
+            -120,
+            2,
+        );
+    }
+
+    #[test]
+    fn normal_quantile_inverts_cdf() {
+        let two = Computable::rational(Rational::new(2));
+        let p = two.clone().pnorm();
+        let seed = BigInt::from((1.9999_f64 * f64::from(1_u32 << 13)).round() as i64);
+        let q = Computable::normal_quantile(p, seed, -13);
+        assert_close(q, two, -120, 2);
+    }
+
     #[test]
     fn add() {
         let three: BigInt = "3".parse().unwrap();

src/real/arithmetic.rs

@@ -1,3 +1,4 @@
+use crate::computable::Precision;
 use crate::{
     CertifiedRealEquality, CertifiedRealOrdering, CertifiedRealSign, Computable, DomainFacts,
     DomainStatus, ExpressionDegree, IdentityFacts, MagnitudeBits, OrderingFacts, PrimitiveFacts,