F32 Normal CDF Instruction and Implementation (#3090)

* F32 Normal CDF Instruction and Implementation

* Linting

* Prop Test and Option Return Removal

* Reverting Change to Cargo.lock

Author: ronald-hood-jr

Cargo.lock

@@ -24,9 +24,11 @@ uint = "0.9"
 proptest = "1.0.0"
 solana-program-test = "1.10.8"
 solana-sdk = "1.10.8"
+libm = "0.2.2"
 
 [lib]
 crate-type = ["cdylib", "lib"]
 
+
 [package.metadata.docs.rs]
 targets = ["x86_64-unknown-linux-gnu"]

libraries/math/Cargo.toml

@@ -41,6 +41,36 @@ pub fn sqrt<T: PrimInt + CheckedShl + CheckedShr>(radicand: T) -> Option<T> {
     Some(result)
 }
 
+/// Calculate the normal cdf of the given number
+///
+/// The approximation is accurate to 3 digits
+///
+/// Code lovingly adapted from the excellent work at:
+///
+/// <https://www.hrpub.org/download/20140305/MS7-13401470.pdf>
+///
+/// The algorithm is based on the implementation in the paper above.
+#[inline(never)]
+pub fn f32_normal_cdf(argument: f32) -> f32 {
+    const PI: f32 = std::f32::consts::PI;
+
+    let mod_argument = if argument < 0.0 {
+        -1.0 * argument
+    } else {
+        argument
+    };
+    let tabulation_numerator: f32 =
+        (1.0 / (1.0 * (2.0 * PI).sqrt())) * (-1.0 * (mod_argument * mod_argument) / 2.0).exp();
+    let tabulation_denominator: f32 =
+        0.226 + 0.64 * mod_argument + 0.33 * (mod_argument * mod_argument + 3.0).sqrt();
+    let y: f32 = 1.0 - tabulation_numerator / tabulation_denominator;
+    if argument < 0.0 {
+        1.0 - y
+    } else {
+        y
+    }
+}
+
 #[cfg(test)]
 mod tests {
     use {super::*, proptest::prelude::*};
@@ -67,4 +97,27 @@ mod tests {
             check_square_root(a as u128);
         }
     }
+
+    fn check_normal_cdf_f32(argument: f32) {
+        let result = f32_normal_cdf(argument);
+        let check_result = 0.5 * (1.0 + libm::erff(argument / std::f32::consts::SQRT_2));
+        let abs_difference: f32 = (result - check_result).abs();
+        assert!(abs_difference <= 0.000_2);
+    }
+
+    #[test]
+    fn test_normal_cdf_f32_min_max() {
+        let test_arguments: [f32; 2] = [f32::MIN, f32::MAX];
+        for i in test_arguments.iter() {
+            check_normal_cdf_f32(*i as f32)
+        }
+    }
+
+    proptest! {
+        #[test]
+        fn test_normal_cdf(a in -1000..1000) {
+
+            check_normal_cdf_f32((a as f32)*0.005);
+        }
+    }
 }

libraries/math/src/approximations.rs

@@ -89,6 +89,14 @@ pub enum MathInstruction {
         argument: f32,
     },
 
+    /// The Normal CDF of a float
+    ///
+    /// No accounts required for this instruction
+    F32NormalCDF {
+        /// The argument
+        argument: f32,
+    },
+
     /// Don't do anything for comparison
     ///
     /// No accounts required for this instruction
@@ -200,6 +208,17 @@ pub fn f32_natural_log(argument: f32) -> Instruction {
     }
 }
 
+/// Create F32 Normal CDF instruction
+pub fn f32_normal_cdf(argument: f32) -> Instruction {
+    Instruction {
+        program_id: id(),
+        accounts: vec![],
+        data: MathInstruction::F32NormalCDF { argument }
+            .try_to_vec()
+            .unwrap(),
+    }
+}
+
 /// Create Noop instruction
 pub fn noop() -> Instruction {
     Instruction {
@@ -347,6 +366,19 @@ mod tests {
         assert_eq!(instruction.program_id, crate::id())
     }
 
+    #[test]
+    fn test_f32_normal_cdf() {
+        let instruction = f32_normal_cdf(f32::MAX);
+        assert_eq!(0, instruction.accounts.len());
+        assert_eq!(
+            instruction.data,
+            MathInstruction::F32NormalCDF { argument: f32::MAX }
+                .try_to_vec()
+                .unwrap()
+        );
+        assert_eq!(instruction.program_id, crate::id())
+    }
+
     #[test]
     fn test_noop() {
         let instruction = noop();

libraries/math/src/instruction.rs

@@ -1,7 +1,11 @@
 //! Program state processor
 
 use {
-    crate::{approximations::sqrt, instruction::MathInstruction, precise_number::PreciseNumber},
+    crate::{
+        approximations::{f32_normal_cdf, sqrt},
+        instruction::MathInstruction,
+        precise_number::PreciseNumber,
+    },
     borsh::BorshDeserialize,
     solana_program::{
         account_info::AccountInfo, entrypoint::ProgramResult, log::sol_log_compute_units, msg,
@@ -33,11 +37,13 @@ fn f32_divide(dividend: f32, divisor: f32) -> f32 {
     dividend / divisor
 }
 
+/// f32_exponentiate
 #[inline(never)]
 fn f32_exponentiate(base: f32, exponent: f32) -> f32 {
     base.powf(exponent)
 }
 
+/// f32_natural_log
 #[inline(never)]
 fn f32_natural_log(argument: f32) -> f32 {
     argument.ln()
@@ -130,6 +136,14 @@ pub fn process_instruction(
             msg!("{}", result as u64);
             Ok(())
         }
+        MathInstruction::F32NormalCDF { argument } => {
+            msg!("Calculating f32 Normal CDF");
+            sol_log_compute_units();
+            let result = f32_normal_cdf(argument);
+            sol_log_compute_units();
+            msg!("{}", result as u64);
+            Ok(())
+        }
         MathInstruction::Noop => {
             msg!("Do nothing");
             msg!("{}", 0_u64);

libraries/math/src/processor.rs

@@ -179,6 +179,21 @@ async fn test_f32_natural_log() {
     banks_client.process_transaction(transaction).await.unwrap();
 }
 
+#[tokio::test]
+async fn test_f32_normal_cdf() {
+    let mut pc = ProgramTest::new("spl_math", id(), processor!(process_instruction));
+
+    // Dial down the BPF compute budget to detect if the operation gets bloated in the future
+    pc.set_compute_max_units(3_100);
+
+    let (mut banks_client, payer, recent_blockhash) = pc.start().await;
+
+    let mut transaction =
+        Transaction::new_with_payer(&[instruction::f32_normal_cdf(0_f32)], Some(&payer.pubkey()));
+    transaction.sign(&[&payer], recent_blockhash);
+    banks_client.process_transaction(transaction).await.unwrap();
+}
+
 #[tokio::test]
 async fn test_noop() {
     let mut pc = ProgramTest::new("spl_math", id(), processor!(process_instruction));