Merge pull request #3042 from dcnorris/numerics
library(numerics), special funs from crate puruspe
This commit is contained in:
Mark Thom
30
Cargo.lock
generated
30
Cargo.lock
generated
@@ -1534,6 +1534,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"
|
||||
@@ -1873,6 +1883,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"
|
||||
@@ -2260,6 +2279,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"
|
||||
@@ -2723,6 +2752,7 @@ dependencies = [
|
||||
"phf",
|
||||
"pprof",
|
||||
"proc-macro2",
|
||||
"puruspe",
|
||||
"quote",
|
||||
"rand",
|
||||
"reqwest",
|
||||
|
||||
@@ -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"
|
||||
|
||||
@@ -766,6 +766,8 @@ The modules that ship with Scryer Prolog are also called
|
||||
* [`uuid`](src/lib/uuid.pl) UUIDv4 generation and hex representation
|
||||
* [`tls`](src/lib/tls.pl)
|
||||
Predicates for negotiating TLS connections explicitly.
|
||||
* [`numerics/special_functions`](src/lib/numerics/special_functions.pl)
|
||||
Predicates for erf, gamma, beta, and related special functions.
|
||||
* [`ugraphs`](src/lib/ugraphs.pl) Graph manipulation library
|
||||
* [`simplex`](src/lib/simplex.pl) Providing `assignment/2`,
|
||||
`transportation/4` and other predicates for solving linear
|
||||
|
||||
@@ -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"))
|
||||
}
|
||||
|
||||
187
src/lib/numerics/quadtests.pl
Normal file
187
src/lib/numerics/quadtests.pl
Normal file
@@ -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)<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)].
|
||||
257
src/lib/numerics/special_functions.pl
Normal file
257
src/lib/numerics/special_functions.pl
Normal file
@@ -0,0 +1,257 @@
|
||||
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
||||
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.
|
||||
|
||||
%% 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).
|
||||
147
src/lib/numerics/testutils.pl
Normal file
147
src/lib/numerics/testutils.pl
Normal file
@@ -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]).
|
||||
@@ -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);
|
||||
|
||||
@@ -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};
|
||||
|
||||
|
||||
121
src/machine/system_calls/special_math.rs
Normal file
121
src/machine/system_calls/special_math.rs
Normal file
@@ -0,0 +1,121 @@
|
||||
use crate::machine::Number;
|
||||
use crate::offset_table::OffsetTable;
|
||||
use crate::Machine;
|
||||
use ordered_float::OrderedFloat;
|
||||
use puruspe::beta::*;
|
||||
use puruspe::error::*;
|
||||
use puruspe::gamma::*;
|
||||
|
||||
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);
|
||||
}
|
||||
}
|
||||
Reference in New Issue
Block a user