e5f4633c39
which included commits to RCS files with non-trunk default branches. git-svn-id: https://yap.svn.sf.net/svnroot/yap/trunk@5 b08c6af1-5177-4d33-ba66-4b1c6b8b522a
669 lines
17 KiB
Prolog
669 lines
17 KiB
Prolog
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
% clp(q,r) version 1.3.3 %
|
|
% %
|
|
% (c) Copyright 1992,1993,1994,1995 %
|
|
% Austrian Research Institute for Artificial Intelligence (OFAI) %
|
|
% Schottengasse 3 %
|
|
% A-1010 Vienna, Austria %
|
|
% %
|
|
% File: arith.pl %
|
|
% Author: Christian Holzbaur christian@ai.univie.ac.at %
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
|
|
|
|
%
|
|
% common code for R,Q, runtime predicates
|
|
%
|
|
% linearize evaluation, collect vars
|
|
%
|
|
% Todo: +) limited encoding length option
|
|
% +) 2 stage compilation: a) linearization
|
|
% b) specialization to R or Q
|
|
%
|
|
%
|
|
|
|
l2conj( [], true).
|
|
l2conj( [X|Xs], Conj) :-
|
|
( Xs = [], Conj = X
|
|
; Xs = [_|_], Conj = (X,Xc), l2conj( Xs, Xc)
|
|
).
|
|
|
|
% ----------------------------------------------------------------------
|
|
|
|
%
|
|
% float/1 coercion is allowed only at the outermost level in Q
|
|
%
|
|
compile_Q( Term, R, Code) :-
|
|
linearize( Term, Res, Linear),
|
|
specialize_Q( Linear, Code, Ct),
|
|
( Res = boolean, Ct = []
|
|
; Res = float(R), Ct = []
|
|
; Res = rat(N,D), Ct = [ putq(D,N,R) ]
|
|
).
|
|
|
|
%
|
|
% assumes normalized params and puts a normalized result
|
|
%
|
|
compile_Qn( Term, R, Code) :-
|
|
linearize( Term, Res, Linear),
|
|
specialize_Qn( Linear, Code, Ct),
|
|
( Res = boolean, Ct = []
|
|
; Res = float(R), Ct = []
|
|
; Res = rat(N,D), Ct = [ putq(D,N,R) ]
|
|
).
|
|
|
|
|
|
compile_case_signum_Qn( Term, Lt,Z,Gt, Code) :-
|
|
linearize( Term, rat(N,_), Linear),
|
|
specialize_Qn( Linear, Code,
|
|
[
|
|
compare( Rel, N, 0),
|
|
( Rel = <, Lt
|
|
; Rel = =, Z
|
|
; Rel = >, Gt
|
|
)
|
|
]).
|
|
|
|
|
|
specialize_Qn( []) --> [].
|
|
specialize_Qn( [Op|Ops]) -->
|
|
specialize_Qn( Op),
|
|
specialize_Qn( Ops).
|
|
%
|
|
specialize_Qn( op_var(rat(N,D),Var)) --> [ Var=rat(N,D) ]. % <--- here is the difference ---
|
|
specialize_Qn( op_integer(rat(I,1),I)) --> [].
|
|
specialize_Qn( op_rat(rat(N,D),N,D)) --> [].
|
|
specialize_Qn( op_float(rat(N,D),X)) --> [], { float_rat( X, N,D) }.
|
|
specialize_Qn( apply(R,Func)) -->
|
|
specialize_Q_fn( Func, R).
|
|
|
|
|
|
specialize_Q( []) --> [].
|
|
specialize_Q( [Op|Ops]) -->
|
|
specialize_Q( Op),
|
|
specialize_Q( Ops).
|
|
%
|
|
specialize_Q( op_var(rat(N,D),Var)) --> [ getq(Var,N,D) ].
|
|
specialize_Q( op_integer(rat(I,1),I)) --> [].
|
|
specialize_Q( op_rat(rat(N,D),N,D)) --> [], { D > 0 }.
|
|
specialize_Q( op_float(rat(N,D),X)) --> [], { float_rat( X, N,D) }.
|
|
specialize_Q( apply(R,Func)) -->
|
|
specialize_Q_fn( Func, R).
|
|
|
|
specialize_Q_fn( +rat(N,D), rat(N,D)) --> [].
|
|
specialize_Q_fn( numer(rat(N,_)), rat(N,1)) --> [].
|
|
specialize_Q_fn( denom(rat(_,D)), rat(D,1)) --> [].
|
|
specialize_Q_fn( -rat(N0,D), rat(N,D)) --> [ N is -N0 ].
|
|
specialize_Q_fn( abs(rat(Nx,Dx)), rat(N,D)) --> [ N is abs(Nx) ], {D=Dx}.
|
|
specialize_Q_fn( signum(rat(Nx,Dx)), rat(N,D)) --> [ signumq( Nx,Dx, N,D) ].
|
|
specialize_Q_fn( floor(rat(Nx,Dx)), rat(N,D)) --> [ floorq( Nx,Dx, N,D) ].
|
|
specialize_Q_fn( ceiling(rat(Nx,Dx)), rat(N,D)) --> [ ceilingq( Nx,Dx, N,D) ].
|
|
specialize_Q_fn( truncate(rat(Nx,Dx)), rat(N,D)) --> [ truncateq( Nx,Dx, N,D) ].
|
|
specialize_Q_fn( round(rat(Nx,Dx)), rat(N,D)) --> [ roundq( Nx,Dx, N,D) ].
|
|
specialize_Q_fn( log(rat(Nx,Dx)), rat(N,D)) --> [ logq( Nx,Dx, N,D) ].
|
|
specialize_Q_fn( exp(rat(Nx,Dx)), rat(N,D)) --> [ expq( Nx,Dx, N,D) ].
|
|
specialize_Q_fn( sin(rat(Nx,Dx)), rat(N,D)) --> [ sinq( Nx,Dx, N,D) ].
|
|
specialize_Q_fn( cos(rat(Nx,Dx)), rat(N,D)) --> [ cosq( Nx,Dx, N,D) ].
|
|
specialize_Q_fn( tan(rat(Nx,Dx)), rat(N,D)) --> [ tanq( Nx,Dx, N,D) ].
|
|
specialize_Q_fn( asin(rat(Nx,Dx)), rat(N,D)) --> [ asinq( Nx,Dx, N,D) ].
|
|
specialize_Q_fn( acos(rat(Nx,Dx)), rat(N,D)) --> [ acosq( Nx,Dx, N,D) ].
|
|
specialize_Q_fn( atan(rat(Nx,Dx)), rat(N,D)) --> [ atanq( Nx,Dx, N,D) ].
|
|
specialize_Q_fn( float(rat(Nx,Dx)), float(F)) --> [ rat_float( Nx,Dx, F) ].
|
|
%
|
|
specialize_Q_fn( rat(Nx,Dx)+rat(Ny,Dy), rat(N,D)) --> [ addq( Nx,Dx, Ny,Dy, N,D) ].
|
|
specialize_Q_fn( rat(Nx,Dx)-rat(Ny,Dy), rat(N,D)) --> [ subq( Nx,Dx, Ny,Dy, N,D) ].
|
|
specialize_Q_fn( rat(Nx,Dx)*rat(Ny,Dy), rat(N,D)) --> [ mulq( Nx,Dx, Ny,Dy, N,D) ].
|
|
specialize_Q_fn( rat(Nx,Dx)/rat(Ny,Dy), rat(N,D)) --> [ divq( Nx,Dx, Ny,Dy, N,D) ].
|
|
specialize_Q_fn( exp(rat(Nx,Dx),rat(Ny,Dy)), rat(N,D)) --> [ expq( Nx,Dx, Ny,Dy, N,D) ].
|
|
specialize_Q_fn( min(rat(Nx,Dx),rat(Ny,Dy)), rat(N,D)) --> [ minq( Nx,Dx, Ny,Dy, N,D) ].
|
|
specialize_Q_fn( max(rat(Nx,Dx),rat(Ny,Dy)), rat(N,D)) --> [ maxq( Nx,Dx, Ny,Dy, N,D) ].
|
|
%
|
|
specialize_Q_fn( rat(Nx,Dx) < rat(Ny,Dy), boolean) --> [ comq( Nx,Dx, Ny,Dy, <) ].
|
|
specialize_Q_fn( rat(Nx,Dx) > rat(Ny,Dy), boolean) --> [ comq( Ny,Dy, Nx,Dx, <) ].
|
|
specialize_Q_fn( rat(Nx,Dx) =< rat(Ny,Dy), boolean) --> [ comq( Nx,Dx, Ny,Dy, Rel), Rel \== (>) ].
|
|
specialize_Q_fn( rat(Nx,Dx) >= rat(Ny,Dy), boolean) --> [ comq( Ny,Dy, Nx,Dx, Rel), Rel \== (>) ].
|
|
specialize_Q_fn( rat(Nx,Dx) =\= rat(Ny,Dy), boolean) --> [ comq( Nx,Dx, Ny,Dy, Rel), Rel \== (=) ].
|
|
specialize_Q_fn( rat(Nx,Dx) =:= rat(Ny,Dy), boolean) -->
|
|
%
|
|
% *normalized* rationals
|
|
%
|
|
( {Nx = Ny} -> [] ; [ Nx = Ny ] ),
|
|
( {Dx = Dy} -> [] ; [ Dx = Dy ] ).
|
|
|
|
% ----------------------------------------------------------------------
|
|
|
|
compile_R( Term, R, Code) :-
|
|
linearize( Term, Res, Linear),
|
|
specialize_R( Linear, Code, Ct),
|
|
( Res == boolean ->
|
|
Ct = [], R = boolean
|
|
; float(Res) ->
|
|
Ct = [ R=Res ]
|
|
;
|
|
Ct = [ R is Res ]
|
|
).
|
|
|
|
compile_case_signum_R( Term, Lt,Z,Gt, Code) :-
|
|
eps( Eps, NegEps),
|
|
linearize( Term, Res, Linear),
|
|
specialize_R( Linear, Code,
|
|
[
|
|
Rv is Res,
|
|
( Rv < NegEps -> Lt
|
|
; Rv > Eps -> Gt
|
|
; Z
|
|
)
|
|
]).
|
|
|
|
specialize_R( []) --> [].
|
|
specialize_R( [Op|Ops]) -->
|
|
specialize_R( Op),
|
|
specialize_R( Ops).
|
|
%
|
|
specialize_R( op_var(Var,Var)) --> [].
|
|
specialize_R( op_integer(R,I)) --> [], { R is float(I) }.
|
|
specialize_R( op_rat(R,N,D)) --> [], { rat_float( N,D, R) }.
|
|
specialize_R( op_float(F,F)) --> [].
|
|
specialize_R( apply(R,Func)) -->
|
|
specialize_R_fn( Func, R).
|
|
|
|
specialize_R_fn( signum(X), S) -->
|
|
( {var(X)} ->
|
|
{Xe=X}
|
|
;
|
|
[ Xe is X ]
|
|
),
|
|
{
|
|
eps( Eps, NegEps)
|
|
},
|
|
[
|
|
( Xe < NegEps -> S = -1.0
|
|
; Xe > Eps -> S = 1.0
|
|
; S = 0.0
|
|
)
|
|
].
|
|
|
|
specialize_R_fn( +X, X) --> [].
|
|
specialize_R_fn( -X, -X) --> [].
|
|
specialize_R_fn( abs(X), abs(X)) --> [].
|
|
specialize_R_fn( floor(X), float(floor(/*float?*/X))) --> [].
|
|
specialize_R_fn( ceiling(X), float(ceiling(/*float?*/X))) --> [].
|
|
specialize_R_fn( truncate(X), float(truncate(/*float?*/X))) --> [].
|
|
specialize_R_fn( round(X), float(round(/*float?*/X))) --> [].
|
|
specialize_R_fn( log(X), log(X)) --> [].
|
|
specialize_R_fn( exp(X), exp(X)) --> [].
|
|
specialize_R_fn( sin(X), sin(X)) --> [].
|
|
specialize_R_fn( cos(X), cos(X)) --> [].
|
|
specialize_R_fn( tan(X), tan(X)) --> [].
|
|
specialize_R_fn( asin(X), asin(X)) --> [].
|
|
specialize_R_fn( acos(X), acos(X)) --> [].
|
|
specialize_R_fn( atan(X), atan(X)) --> [].
|
|
specialize_R_fn( float(X), float(X)) --> [].
|
|
%
|
|
specialize_R_fn( X+Y, X+Y) --> [].
|
|
specialize_R_fn( X-Y, X-Y) --> [].
|
|
specialize_R_fn( X*Y, X*Y) --> [].
|
|
specialize_R_fn( X/Y, X/Y) --> [].
|
|
specialize_R_fn( exp(X,Y), exp(X,Y)) --> [].
|
|
specialize_R_fn( min(X,Y), min(X,Y)) --> [].
|
|
specialize_R_fn( max(X,Y), max(X,Y)) --> [].
|
|
/**/
|
|
%
|
|
% An absolute eps is of course not very meaningful.
|
|
% An eps scaled by the magnitude of the operands participating
|
|
% in the comparison is too expensive to support in Prolog on the
|
|
% other hand ...
|
|
%
|
|
%
|
|
% -eps 0 +eps
|
|
% ---------------[----|----]----------------
|
|
% < 0 > 0
|
|
% <-----------] [----------->
|
|
% =< 0
|
|
% <---------------------]
|
|
% >= 0
|
|
% [--------------------->
|
|
%
|
|
%
|
|
specialize_R_fn( X < Y, boolean) -->
|
|
{
|
|
eps( Eps, NegEps)
|
|
},
|
|
( {X==0} ->
|
|
[ Y > Eps ]
|
|
; {Y==0} ->
|
|
[ X < NegEps ]
|
|
;
|
|
[ X-Y < NegEps ]
|
|
).
|
|
specialize_R_fn( X > Y, boolean) --> specialize_R_fn( Y < X, boolean).
|
|
specialize_R_fn( X =< Y, boolean) -->
|
|
{
|
|
eps( Eps, _)
|
|
},
|
|
[ X-Y < Eps ].
|
|
specialize_R_fn( X >= Y, boolean) --> specialize_R_fn( Y =< X, boolean).
|
|
specialize_R_fn( X =:= Y, boolean) -->
|
|
{
|
|
eps( Eps, NegEps)
|
|
},
|
|
( {X==0} ->
|
|
[ Y >= NegEps, Y =< Eps ]
|
|
; {Y==0} ->
|
|
[ X >= NegEps, X =< Eps ]
|
|
;
|
|
[
|
|
Diff is X-Y,
|
|
Diff =< Eps,
|
|
Diff >= NegEps
|
|
]
|
|
).
|
|
specialize_R_fn( X =\= Y, boolean) -->
|
|
{
|
|
eps( Eps, NegEps)
|
|
},
|
|
[
|
|
Diff is X-Y,
|
|
( Diff < NegEps -> true ; Diff > Eps )
|
|
].
|
|
/**/
|
|
|
|
/**
|
|
%
|
|
% b30427, pp.218
|
|
%
|
|
specialize_R_fn( X > Y, boolean) --> specialize_R_fn( Y < X, boolean).
|
|
specialize_R_fn( X < Y, boolean) -->
|
|
[ scaled_eps(X,Y,E), Y-X > E ].
|
|
|
|
specialize_R_fn( X >= Y, boolean) --> specialize_R_fn( Y =< X, boolean).
|
|
specialize_R_fn( X =< Y, boolean) -->
|
|
[ scaled_eps(X,Y,E), X-Y =< E ]. % \+ >
|
|
|
|
specialize_R_fn( X =:= Y, boolean) -->
|
|
[ scaled_eps(X,Y,E), abs(X-Y) =< E ].
|
|
|
|
specialize_R_fn( X =\= Y, boolean) -->
|
|
[ scaled_eps(X,Y,E), abs(X-Y) > E ].
|
|
|
|
|
|
scaled_eps( X, Y, Eps) :-
|
|
exponent( X, Ex),
|
|
exponent( Y, Ey),
|
|
arith_eps( E),
|
|
Max is max(Ex,Ey),
|
|
( Max < 0 ->
|
|
Eps is E/(1<<Max)
|
|
;
|
|
Eps is E*(1<<Max)
|
|
).
|
|
|
|
exponent( X, E) :-
|
|
A is abs(X),
|
|
float_rat( A, N, D),
|
|
E is msb(N+1)-msb(D).
|
|
**/
|
|
|
|
% ----------------------------------------------------------------------
|
|
|
|
linearize( Term, Res, Linear) :-
|
|
linearize( Term, Res, Vs,[], Lin, []),
|
|
keysort( Vs, Vss),
|
|
( Vss = [], Linear = Lin
|
|
; Vss = [V|Vt], join_vars( Vt, V, Linear, Lin)
|
|
).
|
|
|
|
%
|
|
% flatten the evaluation, collect variables, shared by Q,R,...
|
|
%
|
|
linearize( X, R, [X-R|Vs],Vs) --> {var(X)}, !, [ ].
|
|
linearize( X, R, Vs,Vs) --> {integer(X)}, !, [ op_integer(R,X) ].
|
|
linearize( X, R, Vs,Vs) --> {float(X)}, !, [ op_float(R,X) ].
|
|
linearize( rat(N,D), R, Vs,Vs) --> !, [ op_rat(R,N,D) ].
|
|
linearize( Term, R, V0,V1) -->
|
|
{
|
|
functor( Term, N, A),
|
|
functor( Skeleton, N, A)
|
|
},
|
|
linearize_args( A, Term, Skeleton, V0,V1), [ apply(R,Skeleton) ].
|
|
|
|
linearize_args( 0, _, _, Vs,Vs) --> [].
|
|
linearize_args( N, T, S, V0,V2) -->
|
|
{
|
|
arg( N, T, Arg),
|
|
arg( N, S, Res),
|
|
N1 is N-1
|
|
},
|
|
linearize( Arg, Res, V0,V1),
|
|
linearize_args( N1, T, S, V1,V2).
|
|
|
|
join_vars( [], Y-Ry) --> [ op_var(Ry,Y) ].
|
|
join_vars( [X-Rx|Xs], Y-Ry) -->
|
|
( {X==Y} ->
|
|
{Rx=Ry},
|
|
join_vars( Xs, Y-Ry)
|
|
;
|
|
[ op_var(Ry,Y) ],
|
|
join_vars( Xs, X-Rx)
|
|
).
|
|
|
|
% ---------------------------------- runtime system ---------------------------
|
|
|
|
%
|
|
% C candidate
|
|
%
|
|
limit_encoding_length( 0,D, _, 0,D) :- !. % msb ...
|
|
limit_encoding_length( N,D, Bits, Nl,Dl) :-
|
|
Shift is min(max(msb(abs(N)),msb(D))-Bits,
|
|
min(msb(abs(N)),msb(D))),
|
|
Shift > 0,
|
|
!,
|
|
Ns is N>>Shift,
|
|
Ds is D>>Shift,
|
|
Gcd is gcd(Ns,Ds),
|
|
Nl is Ns//Gcd,
|
|
Dl is Ds//Gcd.
|
|
limit_encoding_length( N,D, _, N,D).
|
|
|
|
|
|
%
|
|
% No longer backconvert to integer
|
|
%
|
|
% putq( 1, N, N) :- !.
|
|
putq( D, N, rat(N,D)).
|
|
|
|
getq( Exp, N,D) :- var( Exp), !,
|
|
raise_exception( instantiation_error(getq(Exp,N,D),1)).
|
|
getq( I, I,1) :- integer(I), !.
|
|
getq( F, N,D) :- float( F), !, float_rat( F, N,D).
|
|
getq( rat(N,D), N,D) :-
|
|
integer( N),
|
|
integer( D),
|
|
D > 0,
|
|
1 =:= gcd(N,D).
|
|
|
|
%
|
|
% actually just a joke to have this stuff in Q ...
|
|
%
|
|
expq( N,D, N1,D1) :- rat_float( N,D, X), F is exp(X), float_rat( F, N1,D1).
|
|
logq( N,D, N1,D1) :- rat_float( N,D, X), F is log(X), float_rat( F, N1,D1).
|
|
sinq( N,D, N1,D1) :- rat_float( N,D, X), F is sin(X), float_rat( F, N1,D1).
|
|
cosq( N,D, N1,D1) :- rat_float( N,D, X), F is cos(X), float_rat( F, N1,D1).
|
|
tanq( N,D, N1,D1) :- rat_float( N,D, X), F is tan(X), float_rat( F, N1,D1).
|
|
asinq( N,D, N1,D1) :- rat_float( N,D, X), F is asin(X), float_rat( F, N1,D1).
|
|
acosq( N,D, N1,D1) :- rat_float( N,D, X), F is acos(X), float_rat( F, N1,D1).
|
|
atanq( N,D, N1,D1) :- rat_float( N,D, X), F is atan(X), float_rat( F, N1,D1).
|
|
|
|
%
|
|
% for integer powers we can do it in Q
|
|
%
|
|
expq( Nx,Dx, Ny,Dy, N,D) :-
|
|
( Dy =:= 1 ->
|
|
( Ny >= 0 ->
|
|
powq( Ny, Nx,Dx, 1,1, N,D)
|
|
;
|
|
Nabs is -Ny,
|
|
powq( Nabs, Nx,Dx, 1,1, N1,D1),
|
|
( N1 < 0 ->
|
|
N is -D1, D is -N1
|
|
;
|
|
N = D1, D = N1
|
|
)
|
|
)
|
|
;
|
|
rat_float( Nx,Dx, Fx),
|
|
rat_float( Ny,Dy, Fy),
|
|
F is exp(Fx,Fy),
|
|
float_rat( F, N, D)
|
|
).
|
|
|
|
%
|
|
% positive integer powers of rational
|
|
%
|
|
powq( 0, _, _, Nt,Dt, Nt,Dt) :- !.
|
|
powq( 1, Nx,Dx, Nt,Dt, Nr,Dr) :- !, mulq( Nx,Dx, Nt,Dt, Nr,Dr).
|
|
powq( N, Nx,Dx, Nt,Dt, Nr,Dr) :-
|
|
N1 is N >> 1,
|
|
( N /\ 1 =:= 0 ->
|
|
Nt1 = Nt, Dt1 = Dt
|
|
;
|
|
mulq( Nx,Dx, Nt,Dt, Nt1,Dt1)
|
|
),
|
|
mulq( Nx,Dx, Nx,Dx, Nxx,Dxx),
|
|
powq( N1, Nxx,Dxx, Nt1,Dt1, Nr,Dr).
|
|
|
|
|
|
/*
|
|
%
|
|
% the choicepoint ruins the party ...
|
|
%
|
|
mulq( Na,Da, Nb,Db, Nc,Dc) :-
|
|
Gcd1 is gcd(Na,Db),
|
|
( Gcd1 =:= 1 -> Na1=Na,Db1=Db; Na1 is Na//Gcd1,Db1 is Db//Gcd1 ),
|
|
Gcd2 is gcd(Nb,Da),
|
|
( Gcd2 =:= 1 -> Nb1=Nb,Da1=Da; Nb1 is Nb//Gcd2,Da1 is Da//Gcd2 ),
|
|
Nc is Na1 * Nb1,
|
|
Dc is Da1 * Db1.
|
|
*/
|
|
mulq( Na,Da, Nb,Db, Nc,Dc) :-
|
|
Gcd1 is gcd(Na,Db),
|
|
Na1 is Na//Gcd1,
|
|
Db1 is Db//Gcd1,
|
|
Gcd2 is gcd(Nb,Da),
|
|
Nb1 is Nb//Gcd2,
|
|
Da1 is Da//Gcd2,
|
|
Nc is Na1 * Nb1,
|
|
Dc is Da1 * Db1.
|
|
|
|
/*
|
|
divq( Na,Da, Nb,Db, Nc,Dc) :-
|
|
Gcd1 is gcd(Na,Nb),
|
|
( Gcd1 =:= 1 -> Na1=Na,Nb1=Nb; Na1 is Na//Gcd1,Nb1 is Nb//Gcd1 ),
|
|
Gcd2 is gcd(Da,Db),
|
|
( Gcd2 =:= 1 -> Da1=Da,Db1=Db; Da1 is Da//Gcd2,Db1 is Db//Gcd2 ),
|
|
( Nb1 < 0 -> % keep denom positive !!!
|
|
Nc is -(Na1 * Db1),
|
|
Dc is Da1 * (-Nb1)
|
|
;
|
|
Nc is Na1 * Db1,
|
|
Dc is Da1 * Nb1
|
|
).
|
|
*/
|
|
divq( Na,Da, Nb,Db, Nc,Dc) :-
|
|
Gcd1 is gcd(Na,Nb),
|
|
Na1 is Na//Gcd1,
|
|
Nb1 is Nb//Gcd1,
|
|
Gcd2 is gcd(Da,Db),
|
|
Da1 is Da//Gcd2,
|
|
Db1 is Db//Gcd2,
|
|
( Nb1 < 0 -> % keep denom positive !!!
|
|
Nc is -(Na1 * Db1),
|
|
Dc is Da1 * (-Nb1)
|
|
;
|
|
Nc is Na1 * Db1,
|
|
Dc is Da1 * Nb1
|
|
).
|
|
|
|
%
|
|
% divq_11( Nb,Db, Nc,Dc) :- divq( 1,1, Nb,Db, Nc,Dc).
|
|
%
|
|
divq_11( Nb,Db, Nc,Dc) :-
|
|
( Nb < 0 -> % keep denom positive !!!
|
|
Nc is -Db,
|
|
Dc is -Nb
|
|
;
|
|
Nc is Db,
|
|
Dc is Nb
|
|
).
|
|
|
|
'divq_-11'( Nb,Db, Nc,Dc) :-
|
|
( Nb < 0 -> % keep denom positive !!!
|
|
Nc is Db,
|
|
Dc is -Nb
|
|
;
|
|
Nc is -Db,
|
|
Dc is Nb
|
|
).
|
|
|
|
/*
|
|
addq( Na,Da, Nb,Db, Nc,Dc) :-
|
|
Gcd1 is gcd(Da,Db),
|
|
( Gcd1 =:= 1 -> % This is the case (for random input) with
|
|
% probability 6/(pi**2).
|
|
Nc is Na*Db + Nb*Da,
|
|
Dc is Da*Db
|
|
;
|
|
T is Na*(Db//Gcd1) + Nb*(Da//Gcd1),
|
|
Gcd2 is gcd(T,Gcd1),
|
|
Nc is T//Gcd2,
|
|
Dc is (Da//Gcd1) * (Db//Gcd2)
|
|
).
|
|
*/
|
|
addq( Na,Da, Nb,Db, Nc,Dc) :-
|
|
Gcd1 is gcd(Da,Db),
|
|
T is Na*(Db//Gcd1) + Nb*(Da//Gcd1),
|
|
Gcd2 is gcd(T,Gcd1),
|
|
Nc is T//Gcd2,
|
|
Dc is (Da//Gcd1) * (Db//Gcd2).
|
|
|
|
/*
|
|
subq( Na,Da, Nb,Db, Nc,Dc) :-
|
|
Gcd1 is gcd(Da,Db),
|
|
( Gcd1 =:= 1 -> % This is the case (for random input) with
|
|
% probability 6/(pi**2).
|
|
Nc is Na*Db - Nb*Da,
|
|
Dc is Da*Db
|
|
;
|
|
T is Na*(Db//Gcd1) - Nb*(Da//Gcd1),
|
|
Gcd2 is gcd(T,Gcd1),
|
|
Nc is T//Gcd2,
|
|
Dc is (Da//Gcd1) * (Db//Gcd2)
|
|
).
|
|
*/
|
|
subq( Na,Da, Nb,Db, Nc,Dc) :-
|
|
Gcd1 is gcd(Da,Db),
|
|
T is Na*(Db//Gcd1) - Nb*(Da//Gcd1),
|
|
Gcd2 is gcd(T,Gcd1),
|
|
Nc is T//Gcd2,
|
|
Dc is (Da//Gcd1) * (Db//Gcd2).
|
|
|
|
comq( Na,Da, Nb,Db, S) :- % todo: avoid multiplication by looking a signs first !!!
|
|
Xa is Na * Db,
|
|
Xb is Nb * Da,
|
|
compare( S, Xa, Xb).
|
|
|
|
minq( Na,Da, Nb,Db, N,D) :-
|
|
comq( Na,Da, Nb,Db, Rel),
|
|
( Rel = =, N=Na, D=Da
|
|
; Rel = <, N=Na, D=Da
|
|
; Rel = >, N=Nb, D=Db
|
|
).
|
|
|
|
maxq( Na,Da, Nb,Db, N,D) :-
|
|
comq( Na,Da, Nb,Db, Rel),
|
|
( Rel = =, N=Nb, D=Db
|
|
; Rel = <, N=Nb, D=Db
|
|
; Rel = >, N=Na, D=Da
|
|
).
|
|
|
|
signumq( N,_, S,1) :-
|
|
compare( Rel, N, 0),
|
|
rel2sig( Rel, S).
|
|
|
|
rel2sig( <, -1).
|
|
rel2sig( >, 1).
|
|
rel2sig( =, 0).
|
|
|
|
|
|
|
|
% -----------------------------------------------------------------------------
|
|
|
|
truncateq( N,D, R,1) :-
|
|
R is N // D.
|
|
|
|
%
|
|
% returns the greatest integral value less than or
|
|
% equal to x. This corresponds to IEEE rounding toward nega-
|
|
% tive infinity
|
|
%
|
|
floorq( N,1, N,1) :- !.
|
|
floorq( N,D, R,1) :-
|
|
( N < 0 ->
|
|
R is N // D - 1
|
|
;
|
|
R is N // D
|
|
).
|
|
|
|
%
|
|
% returns the least integral value greater than or
|
|
% equal to x. This corresponds to IEEE rounding toward posi-
|
|
% tive infinity
|
|
%
|
|
ceilingq( N,1, N,1) :- !.
|
|
ceilingq( N,D, R,1) :-
|
|
( N > 0 ->
|
|
R is N // D + 1
|
|
;
|
|
R is N // D
|
|
).
|
|
|
|
%
|
|
% rounding towards zero
|
|
%
|
|
roundq( N,D, R,1) :-
|
|
% rat_float( N,D, F), % cheating, can do that in Q
|
|
% R is integer(round(F)).
|
|
I is N//D,
|
|
subq( N,D, I,1, Rn,Rd),
|
|
Rna is abs(Rn),
|
|
( comq( Rna,Rd, 1,2, <) ->
|
|
R = I
|
|
; I >= 0 ->
|
|
R is I+1
|
|
;
|
|
R is I-1
|
|
).
|
|
|
|
% ------------------------------- rational -> float -------------------------------
|
|
%
|
|
% The problem here is that SICStus converts BIG fractions N/D into +-nan
|
|
% if it does not fit into a float
|
|
%
|
|
% | ?- X is msb(integer(1.0e+308)).
|
|
% X = 1023
|
|
%
|
|
|
|
rat_float( Nx,Dx, F) :-
|
|
limit_encoding_length( Nx,Dx, 1023, Nxl,Dxl),
|
|
F is Nxl / Dxl.
|
|
|
|
% ------------------------------- float -> rational -------------------------------
|
|
|
|
float_rat( F, N, D) :-
|
|
float_rat( 100, F, F, 1,0,0,1, N0,D0), % at most 100 iterations
|
|
( D0 < 0 -> % sign normalization
|
|
D is -D0,
|
|
N is -N0
|
|
;
|
|
D = D0,
|
|
N = N0
|
|
).
|
|
|
|
float_rat( 0, _, _, Na,_,Da,_, Na,Da) :- !.
|
|
float_rat( _, _, X, Na,_,Da,_, Na,Da) :-
|
|
0.0 =:= abs(X-Na/Da),
|
|
!.
|
|
float_rat( N, F, X, Na,Nb,Da,Db, Nar,Dar) :-
|
|
I is integer(F),
|
|
( I =:= F -> % guard against zero division
|
|
Nar is Na*I+Nb, % 1.0 -> 1/1 and not 0/1 (first iter.) !!!
|
|
Dar is Da*I+Db
|
|
;
|
|
Na1 is Na*I+Nb,
|
|
Da1 is Da*I+Db,
|
|
F1 is 1/(F-I),
|
|
N1 is N-1,
|
|
float_rat( N1, F1, X, Na1,Na,Da1,Da, Nar,Dar)
|
|
).
|
|
|