Add math/special:erf/2

Cargo.toml

@@ -16,13 +16,14 @@ rust-version = "1.77"
 crate-type = ["cdylib", "rlib"]
 
 [features]
-default = ["ffi", "repl", "hostname", "tls", "http", "crypto-full"]
+default = ["ffi", "repl", "hostname", "tls", "http", "crypto-full", "special-math"]
 ffi = ["dep:libffi"]
 repl = ["dep:crossterm", "dep:ctrlc", "dep:rustyline"]
 hostname = ["dep:hostname"]
 tls = ["dep:native-tls"]
 http = ["dep:warp", "dep:reqwest"]
 crypto-full = []
+special-math = ["dep:puruspe"]
 "rust-version-1.80" = []
 
 [build-dependencies]
@@ -62,6 +63,7 @@ scryer-modular-bitfield = "0.11.4"
 num-order = { version = "1.2.0" }
 ordered-float = "4.2.2"
 phf = { version = "0.11", features = ["macros"] }
+puruspe = { version = "0.3.0", optional = true }
 rand = "0.8.5"
 regex = "1.10.6"
 ring = { version = "0.17.8", features = ["wasm32_unknown_unknown_js"] }

build/instructions_template.rs

@@ -515,6 +515,9 @@ enum SystemClauseType {
     #[cfg(feature = "crypto-full")]
     #[strum_discriminants(strum(props(Arity = "6", Name = "$crypto_data_decrypt")))]
     CryptoDataDecrypt,
+    #[cfg(feature = "special-math")]
+    #[strum_discriminants(strum(props(Arity = "2", Name = "$erf")))]
+    Erf,
     #[strum_discriminants(strum(props(Arity = "4", Name = "$ed25519_sign_raw")))]
     Ed25519SignRaw,
     #[strum_discriminants(strum(props(Arity = "4", Name = "$ed25519_verify_raw")))]
@@ -1930,6 +1933,12 @@ fn generate_instruction_preface() -> TokenStream {
                         functor!(atom!("call"), [atom(name), fixnum(arity)])
                     }
                     //
+                    #[cfg(feature = "special-math")]
+                    &Instruction::CallErf => {
+                        let (name, arity) = self.to_name_and_arity();
+                        functor!(atom!("call"), [atom(name), fixnum(arity)])
+                    }
+                    //
                     &Instruction::ExecuteAtomChars |
                     &Instruction::ExecuteAtomCodes |
                     &Instruction::ExecuteAtomLength |
@@ -2168,6 +2177,12 @@ fn generate_instruction_preface() -> TokenStream {
                         functor!(atom!("execute"), [atom(name), fixnum(arity)])
                     }
                     //
+                    #[cfg(feature = "special-math")]
+                    &Instruction::ExecuteErf => {
+                        let (name, arity) = self.to_name_and_arity();
+                        functor!(atom!("execute"), [atom(name), fixnum(arity)])
+                    }
+                    //
                     &Instruction::Deallocate => {
                         functor!(atom!("deallocate"))
                     }

src/lib/math/special.pl

@@ -0,0 +1,80 @@
+/** Special math functions in the Error, Gamma and Beta families
+
+The underlying Rust implementations come from the puruspe crate.
+*/
+
+:- module(special, [
+              erf/2
+                  %,erfc/2
+                  %,inverf/2
+                  %,inverfc/2
+                  %,gamma/2
+                  %,gammq/3
+                  %,invgammp/3
+                  %,ln_gamma/2
+                  %,beta/3
+                  %,betai/4
+                  %,invbetai/4
+          ]).
+
+%% erf(+X, -Erf)
+%
+% Erf is erf(X).
+erf(X, Erf) :- '$erf'(X, Erf).
+
+%% TODO: erfc(+X, -Erfc)
+%
+% Erfc is erfc(X).
+
+%% TODO: inverf(+X, -InvErf)
+%
+% InvErf is erf⁻¹(X). % TODO: Check this; puruspe docs have a mistake.
+
+%% TODO: inverfc(+X, -InvErfc)
+%
+% InvErfc is erfc⁻¹(X). % TODO: Check this; puruspe docs have a mistake.
+
+%% TODO: gamma(+X, -Gamma)
+%
+% Gamma is Γ(X).
+
+%% TODO: gammp(+A, +X, -P)
+%
+% P is P(A,X), the regularized _lower_ incomplete gamma function.
+% X ≥ 0 is the upper limit of integration
+% A > 0 is the shape parameter
+
+%% TODO: gammq(+A, +X, -Q)
+%
+% Q is Q(A,X) := Γ(A,X)/Γ(A), the regularized _upper_ incomplete gamma
+% function.
+% X ≥ 0 is the lower limit of integration
+% A > 0 is the shape parameter
+
+%% TODO: invgammp(+P, +A, -X)
+%
+% X is the unique solution of P = P(A,X), where P(-,-) is the
+% regularized lower incomplete gamma function.
+% P ∈ [0,1] is a probability
+% A > 0 is the _shape parameter_
+
+%% TODO: ln_gamma(+Z, -LnGamma)
+%
+% LnGamma is ln(Γ(Z)), the natural logarithm of Γ(Z).
+
+%% TODO: beta(+Z, +W, -B)
+%
+% B is B(Z,W) := Γ(Z)*Γ(W)/Γ(Z+W)
+
+%% TODO: betai(+A, +B, +X, -I)
+%
+% I is Iₓ(A,B) := B(X;A,B)/B(A,B), the regularized incomplete beta function
+% A > 0 and B > 0 are the first and second shape parameters
+% X ∈ [0,1] is the upper limit of integration.
+
+%% TODO: invbetai(+P, +A, +B, -X)
+%
+% X is the unique solution of P = Iₓ(A,B) := B(X;A,B)/B(A,B)
+% P ∈ [0,1] is a probability
+% A > 0 and B > 0 are the first and second shape parameters
+% X ∈ [0,1] is the upper limit of integration.

src/machine/dispatch.rs

@@ -4523,6 +4523,16 @@ impl Machine {
                         self.crypto_data_decrypt();
                         step_or_fail!(self, self.machine_st.p = self.machine_st.cp);
                     }
+                    #[cfg(feature = "special-math")]
+                    &Instruction::CallErf => {
+                        self.erf();
+                        step_or_fail!(self, self.machine_st.p += 1);
+                    }
+                    #[cfg(feature = "special-math")]
+                    &Instruction::ExecuteErf => {
+                        self.erf();
+                        step_or_fail!(self, self.machine_st.p = self.machine_st.cp);
+                    }
                     &Instruction::CallCryptoCurveScalarMult => {
                         self.crypto_curve_scalar_mult();
                         step_or_fail!(self, self.machine_st.p += 1);

src/machine/system_calls.rs

@@ -83,6 +83,9 @@ use sha3::{Sha3_224, Sha3_256, Sha3_384, Sha3_512};
 
 use crrl::{ed25519, secp256k1, x25519};
 
+#[cfg(feature = "special-math")]
+use puruspe::error::erf;
+
 #[cfg(feature = "tls")]
 use native_tls::{Identity, TlsAcceptor, TlsConnector};
 
@@ -7748,6 +7751,23 @@ impl Machine {
         );
     }
 
+    #[cfg(feature = "special-math")]
+    #[inline(always)]
+    pub(crate) fn erf(&mut self) {
+        let x = self.deref_register(1);
+        let x = match Number::try_from(x) {
+            Ok(Number::Float(n)) => n.into_inner(),
+            Ok(Number::Fixnum(n)) => n.get_num() as f64,
+            Ok(Number::Integer(n)) => n.to_f64().value(),
+            _ => {
+                unreachable!()
+            }
+        };
+        let erf_x = float_alloc!(erf(x), self.machine_st.arena);
+        let return_value = self.deref_register(2);
+        self.machine_st.unify_f64(erf_x, return_value);
+    }
+
     #[inline(always)]
     pub(crate) fn crypto_curve_scalar_mult(&mut self) {
         let stub_gen = || functor_stub(atom!("crypto_curve_scalar_mult"), 4);