library(numerics), special funs from crate puruspe
This commit is contained in:
David C. Norris
187
src/lib/numerics/.ediprolog-temp
Normal file
187
src/lib/numerics/.ediprolog-temp
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)].
|
||||
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)].
|
||||
265
src/lib/numerics/special_functions.pl
Normal file
265
src/lib/numerics/special_functions.pl
Normal file
@@ -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).
|
||||
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]).
|
||||
Reference in New Issue
Block a user