"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"
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"
"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"
"phf",
"pprof",
"proc-macro2",
+ "puruspe",
"quote",
"rand",
"reqwest",
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"
#[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")))]
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 |
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"))
}
--- /dev/null
+% 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)<epsilon).
+% CHECKING.. (?-A=1.5,B=0.7,invgammp(A,B,C),gamma_P_Q(A,C,D,E),abs(B-D)<epsilon).
+ true.
+
+check_module_quads(Module, Quads) :-
+ use_module(Module),
+ read_quads(Module, Quads),
+ zip(Qs, ADs, Quads),
+ length(Qs, NQ),
+ format("% Checking ~d quads ..~n", [NQ]),
+ maplist(check_qu_ad, Qs, ADs).
+
+read_quads(Module, Quads) :-
+ module_terms(Module, Terms),
+ terms_quads(Terms, Quads).
+
+module_terms(Module, Terms) :-
+ module_file(Module, File),
+ setup_call_cleanup(
+ open(File, read, Stream, [type(text)]),
+ read_terms(Stream, Terms),
+ close(Stream)
+ ).
+
+module_file(Module, File) :- atom_concat(Module, '.pl', File).
+
+% Given a list of terms, filter out the predicate clauses.
+% TODO: Arg 1 is really a list of Term-VarNames _pairs_;
+% it would be very nice to find a less unsightly
+% name than 'TermVN' for these!
+terms_quads([Term|Terms], Quads) :-
+ ( term_type(Term, clause) -> 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)].
--- /dev/null
+% 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)<epsilon).
+% CHECKING.. (?-A=1.5,B=0.7,invgammp(A,B,C),gamma_P_Q(A,C,D,E),abs(B-D)<epsilon).
+ true.
+
+check_module_quads(Module, Quads) :-
+ use_module(Module),
+ read_quads(Module, Quads),
+ zip(Qs, ADs, Quads),
+ length(Qs, NQ),
+ format("% Checking ~d quads ..~n", [NQ]),
+ maplist(check_qu_ad, Qs, ADs).
+
+read_quads(Module, Quads) :-
+ module_terms(Module, Terms),
+ terms_quads(Terms, Quads).
+
+module_terms(Module, Terms) :-
+ module_file(Module, File),
+ setup_call_cleanup(
+ open(File, read, Stream, [type(text)]),
+ read_terms(Stream, Terms),
+ close(Stream)
+ ).
+
+module_file(Module, File) :- atom_concat(Module, '.pl', File).
+
+% Given a list of terms, filter out the predicate clauses.
+% TODO: Arg 1 is really a list of Term-VarNames _pairs_;
+% it would be very nice to find a less unsightly
+% name than 'TermVN' for these!
+terms_quads([Term|Terms], Quads) :-
+ ( term_type(Term, clause) -> 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)].
--- /dev/null
+/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+ 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).
--- /dev/null
+/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+ 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]).
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);
use crrl::{ed25519, secp256k1, x25519};
+pub(crate) mod special_math;
+
#[cfg(feature = "tls")]
use native_tls::{Identity, TlsAcceptor, TlsConnector};
--- /dev/null
+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);
+ }
+
+}
projects / scryer-prolog.git / commitdiff