From: David C. Norris Date: Sun, 10 Aug 2025 14:01:01 +0000 (-0400) Subject: library(numerics), special funs from crate puruspe X-Git-Tag: v0.10.0~28^2~4 X-Git-Url: https://git.sagredo.dev/?a=commitdiff_plain;h=af3792f79a4827eb26157b6511a883334303fc60;p=scryer-prolog.git library(numerics), special funs from crate puruspe --- diff --git a/Cargo.lock b/Cargo.lock index 0feaf573..d3452c8d 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -1528,6 +1528,16 @@ dependencies = [ "cpufeatures", ] +[[package]] +name = "lambert_w" +version = "1.2.24" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "05dcaea43deea2259ce41aa33abf6e9ec8bc3e7c65bb222f2de618e10191b132" +dependencies = [ + "num-complex", + "num-traits", +] + [[package]] name = "lazy_static" version = "1.5.0" @@ -1867,6 +1877,15 @@ version = "0.3.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "61807f77802ff30975e01f4f071c8ba10c022052f98b3294119f3e615d13e5be" +[[package]] +name = "num-complex" +version = "0.4.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "73f88a1307638156682bada9d7604135552957b7818057dcef22705b4d509495" +dependencies = [ + "num-traits", +] + [[package]] name = "num-format" version = "0.4.4" @@ -2254,6 +2273,16 @@ dependencies = [ "yansi", ] +[[package]] +name = "puruspe" +version = "0.4.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e2692b0d5427c05a258862ef7e4ee9bbdb03b97554fae3b7fcf6b4ac4a47e034" +dependencies = [ + "lambert_w", + "num-complex", +] + [[package]] name = "quick-xml" version = "0.26.0" @@ -2716,6 +2745,7 @@ dependencies = [ "phf", "pprof", "proc-macro2", + "puruspe", "quote", "rand", "reqwest", diff --git a/Cargo.toml b/Cargo.toml index 67f65bc9..fb4bde0a 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -62,6 +62,7 @@ scryer-modular-bitfield = "0.11.4" num-order = { version = "1.2.0" } ordered-float = "5.0.0" phf = { version = "0.11", features = ["macros"] } +puruspe = "0.4.1" rand = "0.8.5" ring = { version = "0.17.8", features = ["wasm32_unknown_unknown_js"] } ripemd = "0.1.3" diff --git a/build/instructions_template.rs b/build/instructions_template.rs index 1cf004b1..91144ea9 100644 --- a/build/instructions_template.rs +++ b/build/instructions_template.rs @@ -517,6 +517,30 @@ enum SystemClauseType { #[cfg(feature = "crypto-full")] #[strum_discriminants(strum(props(Arity = "6", Name = "$crypto_data_decrypt")))] CryptoDataDecrypt, + #[strum_discriminants(strum(props(Arity = "3", Name = "$beta")))] + Beta, + #[strum_discriminants(strum(props(Arity = "4", Name = "$betai")))] + BetaI, + #[strum_discriminants(strum(props(Arity = "4", Name = "$invbetai")))] + InvBetaI, + #[strum_discriminants(strum(props(Arity = "2", Name = "$gamma")))] + Gamma, + #[strum_discriminants(strum(props(Arity = "2", Name = "$ln_gamma")))] + LnGamma, + #[strum_discriminants(strum(props(Arity = "3", Name = "$gammp")))] + GammP, + #[strum_discriminants(strum(props(Arity = "3", Name = "$gammq")))] + GammQ, + #[strum_discriminants(strum(props(Arity = "3", Name = "$invgammp")))] + InvGammP, + #[strum_discriminants(strum(props(Arity = "2", Name = "$erf")))] + Erf, + #[strum_discriminants(strum(props(Arity = "2", Name = "$erfc")))] + Erfc, + #[strum_discriminants(strum(props(Arity = "2", Name = "$inverf")))] + InvErf, + #[strum_discriminants(strum(props(Arity = "2", Name = "$inverfc")))] + InvErfc, #[strum_discriminants(strum(props(Arity = "4", Name = "$ed25519_sign_raw")))] Ed25519SignRaw, #[strum_discriminants(strum(props(Arity = "4", Name = "$ed25519_verify_raw")))] @@ -1899,6 +1923,21 @@ fn generate_instruction_preface() -> TokenStream { functor!(atom!("call"), [atom_as_cell(name), fixnum(arity)]) } // + &Instruction::CallBeta | + &Instruction::CallBetaI | + &Instruction::CallInvBetaI | + &Instruction::CallGamma | + &Instruction::CallGammP | + &Instruction::CallGammQ | + &Instruction::CallInvGammP | + &Instruction::CallLnGamma | + &Instruction::CallErf | + &Instruction::CallErfc | + &Instruction::CallInvErf | + &Instruction::CallInvErfc => { + let (name, arity) = self.to_name_and_arity(); + functor!(atom!("call"), [atom_as_cell(name), fixnum(arity)]) + } &Instruction::ExecuteAtomChars | &Instruction::ExecuteAtomCodes | &Instruction::ExecuteAtomLength | @@ -2142,6 +2181,21 @@ fn generate_instruction_preface() -> TokenStream { functor!(atom!("execute"), [atom_as_cell(name), fixnum(arity)]) } // + &Instruction::ExecuteBeta | + &Instruction::ExecuteBetaI | + &Instruction::ExecuteInvBetaI | + &Instruction::ExecuteGamma | + &Instruction::ExecuteGammP | + &Instruction::ExecuteGammQ | + &Instruction::ExecuteInvGammP | + &Instruction::ExecuteLnGamma | + &Instruction::ExecuteErf | + &Instruction::ExecuteErfc | + &Instruction::ExecuteInvErf | + &Instruction::ExecuteInvErfc => { + let (name, arity) = self.to_name_and_arity(); + functor!(atom!("execute"), [atom_as_cell(name), fixnum(arity)]) + } &Instruction::Deallocate => { functor!(atom!("deallocate")) } diff --git a/src/lib/numerics/.ediprolog-temp b/src/lib/numerics/.ediprolog-temp new file mode 100644 index 00000000..0f8e739f --- /dev/null +++ b/src/lib/numerics/.ediprolog-temp @@ -0,0 +1,187 @@ +% Efforts toward literate tests with quads + +:- use_module(library(iso_ext)). +:- use_module(library(pio)). +:- use_module(library(lists)). +:- use_module(library(dcgs)). +:- use_module(library(format)). +:- use_module(library(reif)). +:- use_module(library(debug)). +:- use_module(library(lambda)). +:- use_module(library(error)). +:- use_module(library(time)). + +:- use_module(testutils). +:- use_module(special_functions). + +portray_term(Stream) :- + read_term(Stream, Term, []), + portray_clause(Term). + +?- check_module_quads(special_functions, _). +% Checking 11 quads .. +% CHECKING.. (?-A=0.6174468790806071,erf(A,A),B is-A,erf(B,B)). +% CHECKING.. (?-try_falsify(odd_t(erf,real(A)))). +% CHECKING.. (?-witness(odd_t(erf,real(A),false))). +% CHECKING.. (?-witness((real(A),erf(A,B),erf(-A,C),abs(B+C)>0))). +% CHECKING.. (?-length(A,B)). +% CHECKING.. (?-real(A),erf(A,B),erfc(A,C),abs(B+C-1)>epsilon). +% CHECKING.. (?-try_falsify(δ_inverses_t(40*epsilon,erf,inverf,interval(-2,2,A)))). +% CHECKING.. (?-try_falsify(δ_inverses_t(40*epsilon,erfc,inverfc,interval(-2,2,A)))). +% CHECKING.. (?-A=10,B is A+1,gamma(B,C),int_realfact(A,D)). +% CHECKING.. (?-gamma_P_Q(1.2,2.3,A,B),abs(A+B-1) terms_quads(Terms, Quads) + ; Quads = [Term|Quads_], + terms_quads(Terms, Quads_) + ). +terms_quads([], []). + +term_type(Term-_, Type) :- + ( Term = (?- _) -> Type = query + ; Term = (_,_) -> Type = answer_description + ; Term = (_;_) -> Type = answer_description + ; Term == true -> Type = answer_description + ; Term == false -> Type = answer_description + ; Type = clause + ). + +?- term_type(test("erf is odd",try_falsify(odd_t(erf,real(_L))))-_, Type). + Type = clause. + +read_terms(Stream, Terms) :- + read_terms_(Stream, [], Terms). +read_terms_(Stream, Terms0, Terms) :- + Options = [variable_names(VarNames)], + read_term(Stream, Term, Options), + ( Term = end_of_file -> reverse(Terms0, Terms) + ; read_terms_(Stream, [Term-VarNames|Terms0], Terms) + ). + +%% zip(+Xs, +Ys, ?XYs) +%% zip(?Xs, ?Ys, +XYs) +% +% List XYs interleaves same-length lists Xs and Ys. +zip([X|Xs], [Y|Ys], [X,Y|XYs]) :- + zip(Xs, Ys, XYs). +zip([], [], []). + +?- zip(Xs, Ys, XYs). % MGQ does not terminate + error('$interrupt_thrown',repl/0). + +% The following suggested by Ulrich via Quad Works chat +?- zip(Xs, Ys, XYs), false. % loops + +?- zip(X, [4,5,6], [1,4,2,5,3,6]). + X = [1,2,3]. + +?- zip([1,2,3], Y, [1,4,2,5,3,6]). + Y = [4,5,6]. + +?- zip([1,2,3], [4,5,6], Z). + Z = [1,4,2,5,3,6]. + +?- zip(Xs, Ys, [1,4,2,5,3,6]). + Xs = [1,2,3], Ys = [4,5,6]. + +% 3. Demonstrate checking 1 quad, the top two elements of a QAs list. +check_qu_ad(Q-QVN, A-AVN) :- + Q = ?-(G), + phrase(portray_clause_(Q), LitQ), % NB: LitQ terminates w/ newline + format("% CHECKING.. ",[]), + ( A == true -> call(G) + ; A == false -> ( call(G) -> false + ; true + ) + ; phrase(unconj(A), As) -> + ( length(As, N), + n_answers(N, A, AVN, ADs), + n_answers(N, G, QVN, Answers), + maplist(contains, ADs, Answers) + ) + ; % Otherwise, we have the ',' case of a solitary answer + call(G), + call(A), + QVN == AVN + ), + format("~s", [LitQ]). + +% Answer-description AD (qua set-of-bindings) contains Answer. +contains(AD, Answer) :- append(Answer, _, AD). + +?- contains(['Xs'=[C],'L'=1,'_A'=C,'_B'=D], ['Xs'=[A],'L'=1]). + C = A. + +?- check_qu_ad((?-length(_F,_G))-['Xs'=_F,'L'=_G],(_H=[],_I=0;_H=[_J],_I=1;_H=[_J,_K],_I=2;...)-['Xs'=_H,'L'=_I,'_A'=_J,'_B'=_K]). +% CHECKING.. (?-length(A,B)). + _F = [_A,_B], _G = 2, _H = [_J,_K], _I = 2. + +% Unravel the nested (;)/1 applications of multiple-AD structures. +unconj(Conj) --> { Conj = (Elt;Conj_) }, + [Elt], + unconj(Conj_). +unconj(...) --> []. + +?- phrase(unconj((_H=[],_I=0;_H=[_J],_I=1;_H=[_J,_K],_I=2;...)), List). + List = [(_H=[],_I=0),(_H=[_J],_I=1),(_H=[_J,_K],_I=2)]. + +empty_anstack :- + ( retract('$anstack'(_)), fail + ; asserta('$anstack'([])) + ). + +push(VN) :- + retract('$anstack'(As)), + asserta('$anstack'([VN|As])). + +backtrack(N) :- + ( '$anstack'(Ans), + length(Ans, N) -> true + ; fail + ). + +n_answers(N, G, VN, ADs) :- + must_be(integer, N), + ( N > 0 -> n_answers_(N, G, VN, ADs) + ; domain_error(not_less_than_zero, N, n_answers/4) + ). + +n_answers_(N, G, VN, ADs) :- + empty_anstack, + call(G), push(VN), + backtrack(N), + !, + retract('$anstack'(As)), + reverse(As, ADs). + +?- n_answers(3, length(Xs, L), ('Xs'=Xs,'Len'=L), ADs). + Xs = [_A,_B], L = 2, ADs = [('Xs'=[],'Len'=0),('Xs'=[_C],'Len'=1),('Xs'=[_D,_E],'Len'=2)]. diff --git a/src/lib/numerics/quadtests.pl b/src/lib/numerics/quadtests.pl new file mode 100644 index 00000000..0f8e739f --- /dev/null +++ b/src/lib/numerics/quadtests.pl @@ -0,0 +1,187 @@ +% Efforts toward literate tests with quads + +:- use_module(library(iso_ext)). +:- use_module(library(pio)). +:- use_module(library(lists)). +:- use_module(library(dcgs)). +:- use_module(library(format)). +:- use_module(library(reif)). +:- use_module(library(debug)). +:- use_module(library(lambda)). +:- use_module(library(error)). +:- use_module(library(time)). + +:- use_module(testutils). +:- use_module(special_functions). + +portray_term(Stream) :- + read_term(Stream, Term, []), + portray_clause(Term). + +?- check_module_quads(special_functions, _). +% Checking 11 quads .. +% CHECKING.. (?-A=0.6174468790806071,erf(A,A),B is-A,erf(B,B)). +% CHECKING.. (?-try_falsify(odd_t(erf,real(A)))). +% CHECKING.. (?-witness(odd_t(erf,real(A),false))). +% CHECKING.. (?-witness((real(A),erf(A,B),erf(-A,C),abs(B+C)>0))). +% CHECKING.. (?-length(A,B)). +% CHECKING.. (?-real(A),erf(A,B),erfc(A,C),abs(B+C-1)>epsilon). +% CHECKING.. (?-try_falsify(δ_inverses_t(40*epsilon,erf,inverf,interval(-2,2,A)))). +% CHECKING.. (?-try_falsify(δ_inverses_t(40*epsilon,erfc,inverfc,interval(-2,2,A)))). +% CHECKING.. (?-A=10,B is A+1,gamma(B,C),int_realfact(A,D)). +% CHECKING.. (?-gamma_P_Q(1.2,2.3,A,B),abs(A+B-1) terms_quads(Terms, Quads) + ; Quads = [Term|Quads_], + terms_quads(Terms, Quads_) + ). +terms_quads([], []). + +term_type(Term-_, Type) :- + ( Term = (?- _) -> Type = query + ; Term = (_,_) -> Type = answer_description + ; Term = (_;_) -> Type = answer_description + ; Term == true -> Type = answer_description + ; Term == false -> Type = answer_description + ; Type = clause + ). + +?- term_type(test("erf is odd",try_falsify(odd_t(erf,real(_L))))-_, Type). + Type = clause. + +read_terms(Stream, Terms) :- + read_terms_(Stream, [], Terms). +read_terms_(Stream, Terms0, Terms) :- + Options = [variable_names(VarNames)], + read_term(Stream, Term, Options), + ( Term = end_of_file -> reverse(Terms0, Terms) + ; read_terms_(Stream, [Term-VarNames|Terms0], Terms) + ). + +%% zip(+Xs, +Ys, ?XYs) +%% zip(?Xs, ?Ys, +XYs) +% +% List XYs interleaves same-length lists Xs and Ys. +zip([X|Xs], [Y|Ys], [X,Y|XYs]) :- + zip(Xs, Ys, XYs). +zip([], [], []). + +?- zip(Xs, Ys, XYs). % MGQ does not terminate + error('$interrupt_thrown',repl/0). + +% The following suggested by Ulrich via Quad Works chat +?- zip(Xs, Ys, XYs), false. % loops + +?- zip(X, [4,5,6], [1,4,2,5,3,6]). + X = [1,2,3]. + +?- zip([1,2,3], Y, [1,4,2,5,3,6]). + Y = [4,5,6]. + +?- zip([1,2,3], [4,5,6], Z). + Z = [1,4,2,5,3,6]. + +?- zip(Xs, Ys, [1,4,2,5,3,6]). + Xs = [1,2,3], Ys = [4,5,6]. + +% 3. Demonstrate checking 1 quad, the top two elements of a QAs list. +check_qu_ad(Q-QVN, A-AVN) :- + Q = ?-(G), + phrase(portray_clause_(Q), LitQ), % NB: LitQ terminates w/ newline + format("% CHECKING.. ",[]), + ( A == true -> call(G) + ; A == false -> ( call(G) -> false + ; true + ) + ; phrase(unconj(A), As) -> + ( length(As, N), + n_answers(N, A, AVN, ADs), + n_answers(N, G, QVN, Answers), + maplist(contains, ADs, Answers) + ) + ; % Otherwise, we have the ',' case of a solitary answer + call(G), + call(A), + QVN == AVN + ), + format("~s", [LitQ]). + +% Answer-description AD (qua set-of-bindings) contains Answer. +contains(AD, Answer) :- append(Answer, _, AD). + +?- contains(['Xs'=[C],'L'=1,'_A'=C,'_B'=D], ['Xs'=[A],'L'=1]). + C = A. + +?- check_qu_ad((?-length(_F,_G))-['Xs'=_F,'L'=_G],(_H=[],_I=0;_H=[_J],_I=1;_H=[_J,_K],_I=2;...)-['Xs'=_H,'L'=_I,'_A'=_J,'_B'=_K]). +% CHECKING.. (?-length(A,B)). + _F = [_A,_B], _G = 2, _H = [_J,_K], _I = 2. + +% Unravel the nested (;)/1 applications of multiple-AD structures. +unconj(Conj) --> { Conj = (Elt;Conj_) }, + [Elt], + unconj(Conj_). +unconj(...) --> []. + +?- phrase(unconj((_H=[],_I=0;_H=[_J],_I=1;_H=[_J,_K],_I=2;...)), List). + List = [(_H=[],_I=0),(_H=[_J],_I=1),(_H=[_J,_K],_I=2)]. + +empty_anstack :- + ( retract('$anstack'(_)), fail + ; asserta('$anstack'([])) + ). + +push(VN) :- + retract('$anstack'(As)), + asserta('$anstack'([VN|As])). + +backtrack(N) :- + ( '$anstack'(Ans), + length(Ans, N) -> true + ; fail + ). + +n_answers(N, G, VN, ADs) :- + must_be(integer, N), + ( N > 0 -> n_answers_(N, G, VN, ADs) + ; domain_error(not_less_than_zero, N, n_answers/4) + ). + +n_answers_(N, G, VN, ADs) :- + empty_anstack, + call(G), push(VN), + backtrack(N), + !, + retract('$anstack'(As)), + reverse(As, ADs). + +?- n_answers(3, length(Xs, L), ('Xs'=Xs,'Len'=L), ADs). + Xs = [_A,_B], L = 2, ADs = [('Xs'=[],'Len'=0),('Xs'=[_C],'Len'=1),('Xs'=[_D,_E],'Len'=2)]. diff --git a/src/lib/numerics/special_functions.pl b/src/lib/numerics/special_functions.pl new file mode 100644 index 00000000..c6e488d3 --- /dev/null +++ b/src/lib/numerics/special_functions.pl @@ -0,0 +1,265 @@ +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + Written 2025 by David C. Norris (david@precisionmethods.guru) + As with all things floating-point, use at your own risk. + Part of Scryer Prolog. +- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + +/** Special math functions in the Error, Gamma and Beta families + +The underlying Rust implementations come from the +[puruspe](https://docs.rs/puruspe/latest/puruspe/) crate. +*/ + +:- module(special_functions, [ + erf/2 + ,erfc/2 + ,inverf/2 + ,inverfc/2 + ,gamma/2 + ,gamma/3 + ,gamma_P_Q/4 + ,invgammp/3 + ,log_gamma/2 + ,beta/3 + ,betai/4 + ,invbetai/4 + ,test/2 + ,test_special_functions/0 + ,try_falsify/1 + ,witness/1 + ]). + +:- use_module(library(numerics/testutils)). + +%% erf(+Xexpr, -Erf) +% +% X is Xexpr ∈ ℝ, Erf is erf(X). +% +% [DLMF §7.2.1](https://dlmf.nist.gov/7.2#E1), +% [`puruspe::error::erf`](https://docs.rs/puruspe/latest/puruspe/error/fn.erf.html) +erf(Xexpr, Erf) :- + X is Xexpr, + builtins:must_be_number(X, erf/2), + '$erf'(X, Erf). + +% Demonstrate the roots of x - erf(x)) +?- X0 = 0.6174468790806071, erf(X0, X0), _X0 is -X0, erf(_X0, _X0). + X0 = 0.6174468790806071, _X0 = -0.6174468790806071. + +% erf is an odd function: +?- try_falsify(odd_t(erf, real(_))). + false. + +% Another way to say the same thing.. +?- witness(odd_t(erf, real(_), false)). + false. + +% ..and yet one more: +?- witness((real(X), erf(X,Erf), erf(-X,_Erf), abs(Erf+_Erf) > 0)). + false. + +% TODO: Remove this general answer description, +% added merely as a quad-check test case: +?- length(Xs, L). + Xs = [], L = 0 +; Xs = [_A], L = 1 +; Xs = [_A,_B], L = 2 +; ... . + +%% erfc(+X, -Erfc) +% +% Erfc is erfc(X) for X ∈ ℝ. +% +% [DLMF §7.2.2](https://dlmf.nist.gov/7.2#E2), +% [`puruspe::error::erfc`](https://docs.rs/puruspe/latest/puruspe/error/fn.erfc.html) +erfc(X, Erfc) :- + builtins:must_be_number(X, erfc/2), + '$erfc'(X, Erfc). + +?- real(X), erf(X, Erf), erfc(X, Erfc), abs(Erf+Erfc-1) > epsilon. + false. + +%% inverf(+ErfX, -X) +% +% X is erf⁻¹(ErfX) for ErfX ∈ (-1,1). +inverf(ErfX, X) :- + builtins:must_be_number(ErfX, inverf/2), + '$inverf'(ErfX, X). + +?- try_falsify(δ_inverses_t(40*epsilon, erf, inverf, interval(-2,2,_))). + false. + +%% inverfc(+ErfcX, -X) +% +% X is erfc⁻¹(ErfcX) for ErfcX ∈ (0,2). +inverfc(ErfcX, X) :- + builtins:must_be_number(ErfcX, inverfc/2), + '$inverfc'(ErfcX, X). + +?- try_falsify(δ_inverses_t(40*epsilon, erfc, inverfc, interval(-2,2,_))). + false. + +%% gamma(+X, -Gamma) +% +% Gamma is Γ(X), the [ordinary] gamma function evaluated at X ∈ ℝ. +% +% [DLMF §5.2.1](https://dlmf.nist.gov/5.2#E1) +% [`puruspe::gamma::gamma`](https://docs.rs/puruspe/latest/puruspe/gamma/fn.gamma.html) +gamma(X, Gamma) :- + builtins:must_be_number(X, gamma/2), + '$gamma'(X, Gamma). + +% Γ(n+1) ≡ n! +?- N = 10, N1 is N+1, gamma(N1, ΓN1), int_realfact(N, Γ11). + N = 10, N1 = 11, ΓN1 = 3628800.0, Γ11 = 3628800.0. + + +%% gamma(+A, +X, -Gamma) +% +% Gamma is Γ(A,X), the upper incomplete gamma function, where A > 0 +% is the shape parameter and X ≥ 0 is the lower limit of integration. +% +% [DLMF §8.2.2](https://dlmf.nist.gov/8.2#E2), +gamma(A, X, Gamma) :- + builtins:must_be_number(A, gammq/3), + builtins:must_be_number(X, gammq/3), + '$gammq'(A, X, Q), + gamma(A, GammaA), + Gamma is Q*GammaA. + +%% gamma_P_Q(+A, +X, -P, -Q) +% +% For shape parameter A > 0 and lower limit of integration X ≥ 0, +% +% * P is γ(A,X)/Γ(X), the regularized _lower_ incomplete gamma function, and +% +% * Q is Γ(A,X)/Γ(X), the regularized _upper_ incomplete gamma function. +% +% [DLMF §8.2.4](https://dlmf.nist.gov/8.2#E4), +% [`puruspe::gammp::gammp`](https://docs.rs/puruspe/latest/puruspe/gamma/fn.gammp.html), +% [`puruspe::gammp::gammq`](https://docs.rs/puruspe/latest/puruspe/gamma/fn.gammq.html) +gamma_P_Q(A, X, P, Q) :- + builtins:must_be_number(A, gamma_P_Q/4), + builtins:must_be_number(X, gamma_P_Q/4), + '$gammp'(A, X, P), + '$gammq'(A, X, Q). + +% P + Q ≈ 1 +?- gamma_P_Q(1.2, 2.3, P, Q), abs(P + Q - 1) < epsilon. + P = 0.8621845438106976, Q = 0.1378154561893024. + +%% invgammp(+A, +P, -X) +% +% Given shape parameter A > 0 and probability P ∈ [0,1), +% +% X is the unique solution of P = P(A,X), where P(-,-) is the +% regularized lower incomplete gamma function. +% +% [`puruspe::gamma::invgammp`](https://docs.rs/puruspe/latest/puruspe/gamma/fn.invgammp.html) +invgammp(A, P, X) :- + builtins:must_be_number(A, invgammp/3), + builtins:must_be_number(P, invgammp/3), + '$invgammp'(P, A, X). + +?- A = 1.5, P = 0.7, invgammp(A, P, X), gamma_P_Q(A, X, P_, _), abs(P-P_) < epsilon. + A = 1.5, P = 0.7, X = 1.8324353915624363, P_ = 0.7000000000000001. + +%% log_gamma(+X, -LogGamma) +% +% LogGamma is ln(Γ(X)), the natural logarithm of Γ(X). +% +% [`puruspe::gamma::ln_gamma`](https://docs.rs/puruspe/latest/puruspe/gamma/fn.ln_gamma.html) +log_gamma(X, LnGamma) :- + builtins:must_be_number(X, log_gamma/2), + '$ln_gamma'(X, LnGamma). + +%% beta(+X, +Y, -B) +% +% B is B(X,Y) ≡ Γ(X)*Γ(Y)/Γ(X+Y) +% +% [DLMF §5.12.1](https://dlmf.nist.gov/5.12#E1) +% [`puruspe::beta::beta`](https://docs.rs/puruspe/latest/puruspe/beta/fn.beta.html) +beta(X, Y, B) :- + builtins:must_be_number(X, beta/3), + builtins:must_be_number(Y, beta/3), + '$beta'(X, Y, B). + +%% betai(+A, +B, +X, -Ix) +% +% Given: +% +% * shape parameters A > 0 and B > 0, +% * upper limit of integration X ∈ [0,1], +% +% Ix is Iₓ(A,B) ≡ B(X;A,B)/B(A,B), the regularized incomplete beta function; +% +% [DLMF §8.17.2](https://dlmf.nist.gov/8.17#E2), +% [`puruspe::beta::betai`](https://docs.rs/puruspe/latest/puruspe/beta/fn.betai.html) +betai(A, B, X, Ix) :- + builtins:must_be_number(A, betai/4), + builtins:must_be_number(B, betai/4), + builtins:must_be_number(X, betai/4), + '$betai'(A, B, X, Ix). + +%% invbetai(+A, +B, +P, -X) +% +% Given: +% +% * shape parameters A > 0 and B > 0, +% * probability P ∈ [0,1], +% +% X ∈ [0,1] is the unique solution of P = Iₓ(A,B) ≡ B(X;A,B)/B(A,B). +% +% [`puruspe::beta::invbetai`](https://docs.rs/puruspe/latest/puruspe/beta/fn.invbetai.html) +invbetai(A, B, P, X) :- + builtins:must_be_number(A, invbetai/4), + builtins:must_be_number(B, invbetai/4), + builtins:must_be_number(P, invbetai/4), + '$invbetai'(P, A, B, X). + + +% ============================== TESTS ============================== + +%% test_special_functions +% +% Run all tests defined in this module. (These tests _succeed_ when +% they find counterexamples, so the 'desirable' result is `false`.) +test_special_functions :- + format("Seeking counterexamples to assertions:~n", []), + test(T, G), format("% ~s ~n", [T]), + call(G). + +% We default to 1M falsification attempts per assertion, and -- more +% importantly -- use the (unexported) testutils:try_falsify_/2, to +% avoid testutils:reproducibly/0 fixing an RNG seed. Thus we obtain +% truly pseudorandom tests untainted by seed-hacking impropriety. +:- meta_predicate(try_falsify(1)). +try_falsify(G) :- testutils:try_falsify_(10^6, G). + +% A unary witness/1 predicate similarly renders queries more concise. +:- meta_predicate(witness(0)). +witness(G) :- testutils:witness(10^6, G). + +:- discontiguous(test/2). + +%% test(+Name, ?Goal) +% +% Tests have the signature established by @bakaq's test_framework, +% each with a user-facing string Name, and a Goal which serves as an +% _assertion_ by succeeding iff the Name'd desirable property holds. +test("erf is odd", try_falsify(odd_t(erf, real(_)))). + +test("pos root of erf(x)-x", \+ (X0 = 0.6174468790806071, erf(X0, X0))). + +test("erfc ≈ 1 - erf", try_falsify(erf_plus_erfc_unity_t(real(_)))). + +erf_plus_erfc_unity_t(Any, T) :- + call_free(Any, X), erf(X, Erf), erfc(X, Erfc), + ( abs(Erf + Erfc - 1) < epsilon -> T = true + ; T = false + ). + +test("inverf ≈ erf⁻¹", + try_falsify(δ_inverses_t(40*epsilon, erf, inverf, interval(-2,2,_)))). + +test("('false' is good)", false). diff --git a/src/lib/numerics/testutils.pl b/src/lib/numerics/testutils.pl new file mode 100644 index 00000000..7238fb14 --- /dev/null +++ b/src/lib/numerics/testutils.pl @@ -0,0 +1,147 @@ +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + Written 2025 by David C. Norris (david@precisionmethods.guru) + As with all things floating-point, use at your own risk. + Part of Scryer Prolog. +- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + +/** Utility predicates for testing numerics + +*/ + +:- module(testutils, [ + try_falsify/2 + ,witness/2 + ,real/1 + ,posreal/1 + ,interval/3 + ,call_free/2 + ,odd_t/3 + ,δ_inverses_t/5 + ,int_realfact/2 + ]). + +:- use_module(library(lists)). +:- use_module(library(random)). +:- use_module(library(error)). +:- use_module(library(format)). + +%% real(-X) +% +% X = tan(U) for uniform U ~ U[-π,π). +real(X) :- + random(U), + X is tan(pi*(U - 0.5)). + +%% posreal(-X) +% +% X = 1/U for uniform U ~ U[0,1). +posreal(X) :- + random(U), + X is 1/U - 1. + +%% interval(+A, +B, -X) +% +% X ~ U[A,B). +interval(A, B, X) :- + random(U), + X is A + (B-A)*U. + +:- meta_predicate(try_falsify(+, 1)). + +% NOTE: NON-reproducible tests would in fact be MORE STRINGENT, +% creating opportunities to detect rare error cases over time. +reproducibly :- set_random(seed(2025)). + +%% try_falsify(+IntExpr, ?G_1) +% +% Make (N is IntExpr) attempts to falsify the partial goal G/1, +% reporting the first counterexample found. +try_falsify(N, G_1) :- reproducibly -> try_falsify_(N, G_1). +try_falsify_(IntExpr, G_1) :- N is IntExpr, must_be(integer, N), N > 0, + ( call(G_1, false) -> counterexample(G_1) + ; N_ is N - 1, + try_falsify_(N_, G_1) + ). + +%% call_free(?G, -V) +% +% Call goal G, the single free variable of which is bound to V. +call_free(G, V) :- term_variables(G, [V]), call(G). + +:- meta_predicate(odd_t(2, 0, ?)). + +:- meta_predicate(witness(+, 0)). + +%% witness(+N, ?OhNo) +% +% Make N attempts to satisfy the goal OhNo, reporting the first +% counterexample found. +witness(N, OhNo) :- N > 0, + ( call(OhNo) -> counterexample(OhNo) + ; N_ is N - 1, + witness(N_, OhNo) + ). + +% The goal below is obviously _designed_ to succeed 10% of the time, +% giving us plenty of chances to see 'counterexamples': +?- witness(10, (random(X), format("% X = ~f~n", [X]), X > 0.9)). +% X = 0.8084510175379878 +% X = 0.3173976152322231 +% X = 0.6016479707924276 +% X = 0.2782828276608216 +% X = 0.5737059731916494 +% X = 0.4300520177737992 +% X = 0.5411038305190763 +% X = 0.6901983097300066 +% X = 0.05013429563996663 +% X = 0.7321078902421636 + false. +% X = 0.6401599325865659 +% X = 0.7024957773981413 +% X = 0.5798797162672757 +% X = 0.3107013295892864 +% X = 0.3792925200660049 +% X = 0.3424593598278811 +% X = 0.0589899862019303 +% X = 0.9692529829158145 +% COUNTEREXAMPLE: user:random(0.9692529829158145),(current_output(user_output),pio:phrase_to_stream(format:format_([%, ,X, ,=, ,~,f,~,n],[0.9692529829158145]),user_output),flush_output(user_output)),0.9692529829158145>0.9 + X = 0.9692529829158145. + +%% odd_t(+F_2, +Any, ?T) +% +% T is the truth-value from testing that function F_2 is +% [odd](https://en.wikipedia.org/wiki/Even_and_odd_functions) at a +% value X obtained via call_free(Any, X). +odd_t(F, Any, T) :- + ( call_free(Any, X), + _X is -X, + call(F, X, Fx), + call(F, _X, _Fx), + _Fx is -Fx -> T = true + ; T = false + ). + +:- meta_predicate(δ_inverses_t(?, 2, 2, 0, ?)). + +%% δ_inverses_t(+Δ, +F_2, +Finv_2, +Any, ?T) +% +% For given functions F_2 and Finv_2, T is the truth-value from +% testing that Finv_2 ∘ F_2 is within Δ of the identity at a value X +% obtained via call_free(Any, X). +δ_inverses_t(Δ, F, Finv, Any, T) :- + ( call_free(Any, X), + call(F, X, Fx), + call(Finv, Fx, X_), + ( abs(X - X_) < Δ -> T = true + ; T = false + ) + ). + +%% int_realfact(+N, -FactN) +% +% FactN is floating-point N! +int_realfact(N, FactN) :- + N > 0, N_ is N - 1, int_realfact(N_, FactN_), FactN is N*FactN_. +int_realfact(0, 1.0). + +counterexample(G) :- format("% COUNTEREXAMPLE: ~w~n", [G]). diff --git a/src/machine/dispatch.rs b/src/machine/dispatch.rs index 8a8921e6..59db4e1c 100644 --- a/src/machine/dispatch.rs +++ b/src/machine/dispatch.rs @@ -4691,6 +4691,102 @@ impl Machine { self.crypto_data_decrypt(); step_or_fail!(self, self.machine_st.p = self.machine_st.cp); } + &Instruction::CallBeta => { + self.beta(); + step_or_fail!(self, self.machine_st.p += 1); + } + &Instruction::ExecuteBeta => { + self.beta(); + step_or_fail!(self, self.machine_st.p = self.machine_st.cp); + } + &Instruction::CallBetaI => { + self.betai(); + step_or_fail!(self, self.machine_st.p += 1); + } + &Instruction::ExecuteBetaI => { + self.betai(); + step_or_fail!(self, self.machine_st.p = self.machine_st.cp); + } + &Instruction::CallInvBetaI => { + self.invbetai(); + step_or_fail!(self, self.machine_st.p += 1); + } + &Instruction::ExecuteInvBetaI => { + self.invbetai(); + step_or_fail!(self, self.machine_st.p = self.machine_st.cp); + } + &Instruction::CallGamma => { + self.gamma(); + step_or_fail!(self, self.machine_st.p += 1); + } + &Instruction::ExecuteGamma => { + self.gamma(); + step_or_fail!(self, self.machine_st.p = self.machine_st.cp); + } + &Instruction::CallLnGamma => { + self.ln_gamma(); + step_or_fail!(self, self.machine_st.p += 1); + } + &Instruction::ExecuteLnGamma => { + self.ln_gamma(); + step_or_fail!(self, self.machine_st.p = self.machine_st.cp); + } + &Instruction::CallGammP => { + self.gammp(); + step_or_fail!(self, self.machine_st.p += 1); + } + &Instruction::ExecuteGammP => { + self.gammp(); + step_or_fail!(self, self.machine_st.p = self.machine_st.cp); + } + &Instruction::CallInvGammP => { + self.invgammp(); + step_or_fail!(self, self.machine_st.p += 1); + } + &Instruction::ExecuteInvGammP => { + self.invgammp(); + step_or_fail!(self, self.machine_st.p = self.machine_st.cp); + } + &Instruction::CallGammQ => { + self.gammq(); + step_or_fail!(self, self.machine_st.p += 1); + } + &Instruction::ExecuteGammQ => { + self.gammq(); + step_or_fail!(self, self.machine_st.p = self.machine_st.cp); + } + &Instruction::CallErf => { + self.erf(); + step_or_fail!(self, self.machine_st.p += 1); + } + &Instruction::ExecuteErf => { + self.erf(); + step_or_fail!(self, self.machine_st.p = self.machine_st.cp); + } + &Instruction::CallErfc => { + self.erfc(); + step_or_fail!(self, self.machine_st.p += 1); + } + &Instruction::ExecuteErfc => { + self.erfc(); + step_or_fail!(self, self.machine_st.p = self.machine_st.cp); + } + &Instruction::CallInvErf => { + self.inverf(); + step_or_fail!(self, self.machine_st.p += 1); + } + &Instruction::ExecuteInvErf => { + self.inverf(); + step_or_fail!(self, self.machine_st.p = self.machine_st.cp); + } + &Instruction::CallInvErfc => { + self.inverfc(); + step_or_fail!(self, self.machine_st.p += 1); + } + &Instruction::ExecuteInvErfc => { + self.inverfc(); + 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); diff --git a/src/machine/system_calls.rs b/src/machine/system_calls.rs index 3f367a39..12094c36 100644 --- a/src/machine/system_calls.rs +++ b/src/machine/system_calls.rs @@ -84,6 +84,8 @@ use sha3::{Sha3_224, Sha3_256, Sha3_384, Sha3_512}; use crrl::{ed25519, secp256k1, x25519}; +pub(crate) mod special_math; + #[cfg(feature = "tls")] use native_tls::{Identity, TlsAcceptor, TlsConnector}; diff --git a/src/machine/system_calls/special_math.rs b/src/machine/system_calls/special_math.rs new file mode 100644 index 00000000..907c9748 --- /dev/null +++ b/src/machine/system_calls/special_math.rs @@ -0,0 +1,123 @@ +use crate::machine::Number; +use crate::Machine; +use crate::offset_table::OffsetTable; // These appear now +use ordered_float::OrderedFloat; // to be required. +use puruspe::error::*; +use puruspe::gamma::*; +use puruspe::beta::*; + +macro_rules! number_as_f64 { + ($self: ident, $reg: literal) => {{ + match Number::try_from(($self.deref_register($reg), &$self.machine_st.arena.f64_tbl)) { + 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!() + } + } + }}; +} + +macro_rules! return_f64_reg { + ($self: ident, $val: ident, $reg: literal) => {{ + let return_value = $self.deref_register($reg); + $self.machine_st.unify_f64($val, return_value); + }}; +} + +impl Machine { + + #[inline(always)] + pub(crate) fn erf(&mut self) { + let x = number_as_f64!(self, 1); + let erf_x = float_alloc!(erf(x), self.machine_st.arena); + return_f64_reg!(self, erf_x, 2); + } + + #[inline(always)] + pub(crate) fn erfc(&mut self) { + let x = number_as_f64!(self, 1); + let erfc_x = float_alloc!(erfc(x), self.machine_st.arena); + return_f64_reg!(self, erfc_x, 2); + } + + #[inline(always)] + pub(crate) fn inverf(&mut self) { + let erf_x = number_as_f64!(self, 1); + let x = float_alloc!(inverf(erf_x), self.machine_st.arena); + return_f64_reg!(self, x, 2); + } + + #[inline(always)] + pub(crate) fn inverfc(&mut self) { + let erfc_x = number_as_f64!(self, 1); + let x = float_alloc!(inverfc(erfc_x), self.machine_st.arena); + return_f64_reg!(self, x, 2); + } + + #[inline(always)] + pub(crate) fn gamma(&mut self) { + let x = number_as_f64!(self, 1); + let gamma_x = float_alloc!(gamma(x), self.machine_st.arena); + return_f64_reg!(self, gamma_x, 2); + } + + #[inline(always)] + pub(crate) fn gammp(&mut self) { + let a = number_as_f64!(self, 1); + let x = number_as_f64!(self, 2); + let gammp_a_x = float_alloc!(gammp(a, x), self.machine_st.arena); + return_f64_reg!(self, gammp_a_x, 3); + } + + #[inline(always)] + pub(crate) fn gammq(&mut self) { + let a = number_as_f64!(self, 1); + let x = number_as_f64!(self, 2); + let gammq_a_x = float_alloc!(gammq(a, x), self.machine_st.arena); + return_f64_reg!(self, gammq_a_x, 3); + } + + #[inline(always)] + pub(crate) fn invgammp(&mut self) { + let p = number_as_f64!(self, 1); + let a = number_as_f64!(self, 2); + let x = float_alloc!(invgammp(p, a), self.machine_st.arena); + return_f64_reg!(self, x, 3); + } + + #[inline(always)] + pub(crate) fn ln_gamma(&mut self) { + let x = number_as_f64!(self, 1); + let ln_gamma_x = float_alloc!(ln_gamma(x), self.machine_st.arena); + return_f64_reg!(self, ln_gamma_x, 2); + } + + #[inline(always)] + pub(crate) fn beta(&mut self) { + let x = number_as_f64!(self, 1); + let y = number_as_f64!(self, 2); + let beta_x_y = float_alloc!(beta(x,y), self.machine_st.arena); + return_f64_reg!(self, beta_x_y, 3); + } + + #[inline(always)] + pub(crate) fn betai(&mut self) { + let a = number_as_f64!(self, 1); + let b = number_as_f64!(self, 2); + let x = number_as_f64!(self, 3); + let betai_a_b_x = float_alloc!(betai(a,b,x), self.machine_st.arena); + return_f64_reg!(self, betai_a_b_x, 4); + } + + #[inline(always)] + pub(crate) fn invbetai(&mut self) { + let a = number_as_f64!(self, 1); + let b = number_as_f64!(self, 2); + let p = number_as_f64!(self, 3); + let x = float_alloc!(invbetai(a,b,p), self.machine_st.arena); + return_f64_reg!(self, x, 4); + } + +}