1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103

#![feature(const_fn_floating_point_arithmetic)]

#[cfg(test)]
mod tests;

mod table;
use table::{TABLE, MAX_Z};

// Maps some z value to the corresponding table index.
const fn z_to_index(z: f32) -> usize {
    (z / MAX_Z * TABLE.len() as f32) as usize
}

// Maps some table index to the corresponding z value.
const fn index_to_z(i: usize) -> f32 {
    i as f32 * MAX_Z / TABLE.len() as f32
}

/// Lookup function for the Z table.
/// Given a z value, returns the corresponding value of
/// the cumulative distribution function of the standard
/// normal distribution.
/// The input values may be negative.
pub const fn lookup(z: f32) -> f32 {
    if z >= 0.0 {
        lookup_index(z_to_index(z))
    } else {
        1.0 - lookup_index(z_to_index(-z))
    }
}

/// Lookup function for the Z table.
/// Given a z value, returns the corresponding value of
/// the cumulative distribution function of the normal distribution
/// with the provided mean and standard derivation.
/// The input values may be negative.
pub const fn lookup_with(z: f32, mean: f32, standard_derivation: f32) -> f32 {
    lookup((z - mean) / standard_derivation)
}

/// Reverse lookup function for the Z table.
/// Given a value of the cumulative distribution function
/// of the standard normal distribution, returns
/// the corresponding z value.
/// Only inputs between and including 0 and 1 are allowed.
pub const fn reverse_lookup(p: f32) -> f32 {
    assert!(0.0 <= p && p <= 1.0);
    if p >= 0.5 {
        index_to_z(reverse_lookup_index(p))
    } else {
        -index_to_z(reverse_lookup_index(1.0 - p))
    }
}

/// Reverse lookup function for the Z table.
/// Given a value of the cumulative distribution function
/// of the normal distribution with the provided mean and
/// standard derivation, returns the corresponding z value.
/// Only inputs between and including 0 and 1 are allowed.
pub const fn reverse_lookup_with(p: f32, mean: f32, standard_derivation: f32) -> f32 {
    reverse_lookup(p) * standard_derivation + mean
}

// Provides a compile time reverse lookup for the lookup table.
const fn reverse_lookup_index(p: f32) -> usize {
    assert!(0.5 <= p && p <= 1.0);
    let mut prev_abs = std::f32::MAX;
    let mut i = 0;
    loop {
        let curr_abs = abs(p - lookup_index(i));
        if prev_abs < curr_abs {
            if prev_abs < curr_abs {
                return i - 1;
            } else {
                return i;
            }
        }
        prev_abs = curr_abs;
        i += 1;
        if i == TABLE.len() {
            return i;
        }
    }
}

// Computes the absolute value.
const fn abs(x: f32) -> f32 {
    if x >= 0.0 {
        x
    } else {
        -x
    }
}

// Lookup and index in internal lookup table
// without panicking at invalid indices.
const fn lookup_index(i: usize) -> f32 {
    if i >= TABLE.len() {
        1.0
    } else {
        TABLE[i]
    }
}