/* LONG.PL : Arbitrary precision rational arithmetic package.
Copyright (C) 1981 - R.A.O'Keefe. Updated: 30 August 82
Designed and written by Richard O'Keefe.
Scenery by Lawrence Byrd.
As "number" is an immutable name in NIP I have replaced it by
"ratnum" throughout. (Ken Johnson 6-5-87)
This package provides arithmetic for arbitrary precision rational
numbers. The normal domain of prolog 'integers' is extended to
full rational 'numbers'. This domain includes all Prolog integers.
The predicate:
ratnum(N)
will recognise any number in this extended domain.
Rational numbers are produced by using the predicates
eval(Command)
eval(Expression,Answer)
Expression can involve any form of rational number, whether such
numbers can be represented by Prolog integer or not. Any form of
number produced as output by "eval" is acceptable as input to it.
For convenience the Answer produced by eval is normalised as follows:
a) Integers X (where |X| <= 99999) are represented as Prolog integers;
b) 1/0, 0/0, -1/0 are represented as infinity, undefined, neginfinity;
c) All other numbers are represented as full rationals in reduced form
i.e. numerator and denominator are relatively prime.
In the current representation, one normalised number will unify with
another (including an integer) iff the two numbers are equal. But it
is better to test for equality between arbitrary numbers by calling
eval(N1=:=N2)
which also handles infinity & undefined, and is guaranteed to work.
Once created, representations of rational numbers can be passed round
your program, used with eval, or printed. The predicate
portray_ratnum(Number)
will pretty-print arbitrary numbers, and will fail for anything
else. In particular, it will not evaluate an expression. (But
eval(write(Expr)) combines evaluation and printing if you want.)
If this is connected up to your general "portray" mechanism, you
will never have to see the internal representation of rationals.
It is ill-advised to write procedures which assume knowledge of
this internal representation as it is subject to change (rarely),
not to mention that such activities are against all the principles
of abstraction and structured programming.
NB Note that eval/1 and eval/2 will only evaluate fully numeric
expressions. If there is some garbage in the expression (such
as an atom) then no evaluation at all occurs and the whole
input expression is returned untouched. If you want to evaluate
mixed symbolic and numeric expressions then use tidy/2 (from
TIDY.PL) which is designed for this purpose.
FIXES
[3 April 81]
Added the functions numer(_), denom(_) and sign(_) to the
evaluator (ie eva2).
[8 April 81]
removed choice-points from comq, and corrected sign(.).
replaced the log routine completely.
[14 April 81]
changed all XXXr routines to XXXn (for Natural or zero)
changed all XXXs routines to XXXm (for Modified (Natural routines))
changed all XXXl routines to XXXz (for the ring Z of integers)
replaced "digits" by "conn" as I've meant to for some time.
removed experimental 'xwd' code which doesn't work compiled. Eheu.
changed estq,chkq,gest to estd,chkd,estg (estimate division digit,
check digit, estimate Gcd) to avoid confusion; they don't use rationals.
rewrote norm_number and renamed it to standardise.
laid the trig routines out in MY style not Lawrence's.
Increased the radix from 10,000 to 100,000 after fixing addn to use
unsigned numbers.
[21 April 81]
Continued tidying things up.
made 0^(1/N) = 0; this was an oversight.
added new xwd(,) code in eva2, and beautified portray_number.
[8 July 81]
fixed mode error bug in eva2(abs(_),_). Foolish oversight.
[9 Sept 81]
fixed negative number bugs in arccos and arcsin. How long have
these been around without anyone (except Bernard) noticing?
Also shunted some cuts around in the same general area.
[13 Sept 81]
corrected typo {da=>Da} in gcdq/4.
[2 Dec 81]
corrected a benign bug in ratnum/5 (100 had been written
where R should have been), and some minor cosmetic changes.
Unified error reporting into long_error[_message]. Added a
few mode declarations for trig functions.
[9 Dec 81]
when writing eval up for EuroSam, discovered that logs aren't
handled properly. Rewrote absq and logq to return 'undefined'
in more cases, instead of failing.
[27 July 82]
changed prnq/3 to portray_ratnum/3 and laid it out properly.
changed prin/1 to putn/1 (put Natural) to avoid conflict elsewhere.
Made this stuff call put/1 where it made sense, and used ASCII
codes instead of strings. Don't know if it matters, really.
Also rewrote arctaneval completely, so that it should succeed in a
few more cases. Really, the trig stuff is PITIFUL. Please, will
someone do a proper job of it (preferrably someone PAID to do it).
[30 August 82]
fixed bug in gcdq/4 so that gcd(1/2,1/4) = 1/4, not 1/2!
[12 July 1983]
arctaneval used to call addn/4, and there isn't any such
predicate. Made it call addn/5.
[6 May 1987 Ken Johnson]
For transportability, Radix set to 100 (will fit inside 16
bits even allowing for a sign bit). Replaced / by "div" in
is-expressions.
*/
/* OPERATORS */
:- op(300, xfx, div). % integer quotient A div B = fix(A/B).
:- op(500, yfx, ++ ).
:- op(500, yfx, -- ).
/* MODES and types */
% The comments at the right give the argument types for each predicate.
% The predicates can of course be called with any arguments, but these
% are the only types they are supposed to work on or deliver.
% ? = any Prolog term, possibly including variables.
% E = an arithmetic Expression, a term.
% A = a Prolog atom (but not an integer).
% I = a Prolog integer. Generally positive, but not always.
% T = a Truth-value, 'true' or 'false'.
% S = a Sign, '+' or '-'. {Sometimes can be 0 or *.}
% R = a Relational operator, {<, =, >; sometimes =<, >=, =/=}
% N = a long positive (Natural) number.
% Q = a rational number.
%% Top Level %%
% ratnum(+) % Q
% eval(+) % E
% eval(+, -) % E Q
% eva2(+, -) % E Q
% relational_op(+, -, -) % R R T
% combine_ops(+, +, +, -) % R R T T
% portray_ratnum(+) % Q
% portray_ratnum(+, +, +) % S N N
% putn(+) % N
%% Conversions %%
% ratnum(+, +, ?, ?, ?), % Q I S N N
% binrad(+, +, -), % I I N
% standardise(+, ?), % Q Q
%
%%% Low Level %%
%
% add(+, +, ?), % Q Q Q
% multiply(+, +, ?), % Q Q Q
% power(+, +, ?), % Q Q Q
%
%%% Rational Arithmetic %%
%
% mod2(+, ?), % Q I
% intq(+, +, -), % Q I Q
%% gcdq(+, +, +, -), % Q Q I Q
%/* invq(+, +, -), % Q I Q */
% mulq(+, +, +, -), % Q Q I Q
% divq(+, +, +, -), % Q Q I Q
% divo(+, +, +, -, -), % Q Q I Q Q
% powq(+, +, +, -), % Q Q I Q
% negq(+, +, -), % Q I Q
% addq(+, +, +, -), % Q Q I Q
% subq(+, +, +, -), % Q Q I Q
% comq(+, +, +, ?), % Q Q I R
% nthq(+, +, +, -), % I Q I Q
% nthn(+, +, +, -), % I N I N
% newton(+, +, +, +, -), % I N N I N
% newton(+, +, +, +, +, -), % R I N N I N
%
% %% Long Arithmetic %%
%
% addz(+,+, +,+, +, -,-), % S N S N I S N
% addn(+, +, +, +, -), % N N I I N
% add1(+, +, -), % N I N
% comz(+,+, +,+, ?), % S N S N R
% comn(+, +, +, ?), % N N R R
% com1(+, +, +, -), % I I R R
% subz(+,+, +,+, +, -,-), % S N S N I S N
% subn(+, +, +, -,-), % N N I S N
% subn(+, +, +, +, -,-), % R N N I S N
% prune(+, -), % N N
% subp(+, +, +, +, -), % N N I I N
% sub1(+, +, -), % N I N
% sign(+, +, -), % S S S
%/* mulz(+,+, +,+, +, -,-), % S N S N I S N */
% muln(+, +, +, -), % N N I N
% muln(+, +, +, +, -), % N N N I N
% mul1(+, +, +, -), % N I I N
% mul1(+, +, +, +, -), % N I I I N
% powz(+,+, +, +, -,-), % S N I I S N
% pown(+, +, +, +, -), % I N N I N
% divz(+,+, +,+, +, -,-, -,-), % S N S N I S N S N
% divn(+, +, +, -, -), % N N I N N
% conn(?, ?, ?), % I N N
% % both +, +, - and -, -, + are used.
% div1(+, +, +, -, -), % N I I N N
% divm(+, +, +, -, -), % N N I N N
% div2(+, +, +, -, -), % N N I N N
% estd(+, +, +, -), % N N I I
% chkd(+, +, +, +, +, -, -), % N N I I I I N
%/* gcdz(+,+, +,+, +, -, -,-, -,-), % S N S N I N S N S N */
% gcdn(+, +, +, -, -, -), % N N I N N N
% gcdn(+, +, +, -), % N N I N
% gcdn(+, +, +, +, -), % R N N I N
% estg(+, +, +, -), % N N I I
%
%%% Logarithms %%
%
% logq(+, +, +, -), % Q Q I Q
% logq(+,+, +,+,+,-), % R R Q Q I Q
% absq(+, -, -), % Q S Q
% logq(+, -,-), % S S N
% oneq(+, -, -), % Q R Q
% ratlog(+, +, +, -), % Q Q I Q
% ratlog(+,+, +,+,+, -), % S S Q Q I Q
% lograt(+,+,+, -,-), % Q Q I N N
% loop(+, +, +, -), % N N I N
% loop(+,+,+,+,-), % N N N I N
% logn(+,+,+,+,-), % Q I Q Q I I
%
%%% Trigonometry %%
%
% sineval(+, -), % Q Q
% coseval(+, -), % Q Q
% taneval(+, -), % Q Q
% arcsineval(+, -), % Q Q
% arccoseval(+, -), % Q Q
% arctaneval(+, -), % Q Q
% arctaneval(+, +, -, -), % N N N N
% sineval1(?, ?), % Q Q
%
%%% Error handing %%
%
% long_error(+, ?), % A ?
% long_error_message(+, -). % A A
%
/* Implementation
The internal representation for rationals is of the form:
ratnum(Sign, Numerator, Denominator)
where
Sign is in {+,-}
Numerator is a list of (Prolog) integers
Denominator is a list of (Prolog) integers
The lists of Prolog integers represent arbitrary precision unsigned
long integers
eg [n0,n1,....,nz]
is n0+R*(n1+R*(....R*nz)...)
where R is the Radix.
The Radix used in the current version is 100. Most of the code
in this module is completely independent of the radix - it all
uses the value passed in by the top level procedures. However the
printing routine currently assumes that the radix is a power of
10 as this makes things easier. In general the radix must be such
that both:
Radix^2 - 1
and Radix*2 + 1
are representable as Prolog integers (which
are 18 bit quantities on the DEC10). This is a little restrictive,
however, and this implementation only assumes that Radix^2 - 1 is
"obtainable" as an intermediate during Prolog arithmetic. On the
DEC10 intermediate results can be 36 bit quantities and so 100
becomes a suitable radix.
The code actually unpacks the number terms into their separate
bits for all the low level operations. At this stage the following
additional number forms are appropriately converted
- (Prolog integers)
infinity - represented as +1/0
neginfinity - represented as -1/0
undefined - represented as 0/0
The treatment of these strange things is not supposed to be
mathematically beautiful, but sensible things happen using
this representation. They are strictly an extension to the
rationals and could be removed (with eval failing should 0
denominator numbers ever get produced) if desired.
Results from eval are normalised before being returned.
This operation reverses the above transformation except that
only integers within the range -99999 to +99999 are turned
back into Prolog integers.
*/
%% TOP LEVEL PREDICATES %%
% Number recognition predicate
ratnum(N) :- integer(N), !.
ratnum(ratnum(_,_,_)) :- !.
ratnum(infinity) :- !.
ratnum(neginfinity) :- !.
ratnum(undefined) :- !.
% Simple eval interpreter with various features.
eval(Var) :- var(Var), !, long_error(eval, Var).
eval(B is Y) :- !, eval(Y, B).
eval(write(Y)) :- !, eval(Y, B), print(B).
eval(even(X)) :- !, eva2(X, A), !, mod2(A, 0).
eval(odd(X)) :- !, eva2(X, A), !, mod2(A, 1).
eval(compare(R,A,B)) :-
!, eva2(_, A), eva2(_, B), comq(A, B, 100, S), !, R=S.
eval(Term) :- Term =.. [F,X,Y], relational_op(F, R, Flag),
!, eva2(X, A), eva2(Y, B), comq(A, B, 100, S), !,
combine_ops(R, S, Flag, true).
mod2(ratnum(_,_,[]), _) :- !, fail.
mod2(ratnum(_,[],_), 0).
mod2(ratnum(_,[L|_],_), M) :- M is L mod 2.
% General evaluation of rational expressions
eval(Exp, Ans) :- % Hope for the best
eva2(Exp, N),
standardise(N, A), !,
Ans = A.
eval(Exp, Exp). % Cannot evaluate so leave alone
% ttynl, display('[Couldn''t evaluate: '),
% print(Exp), ttyput("]"), ttynl, ttynl.
eva2(Var, _) :- var(Var), !, long_error(eva2, Var).
eva2(X+Y, C) :- !, eva2(X, A), eva2(Y, B), addq(A, B, 100, C).
eva2(X-Y, C) :- !, eva2(X, A), eva2(Y, B), subq(A, B, 100, C).
eva2( -Y, C) :- !, eva2(Y, B), negq( B, 100, C).
eva2(X*Y, C) :- !, eva2(X, A), eva2(Y, B), mulq(A, B, 100, C).
eva2(X/Y, C) :- !, eva2(X, A), eva2(Y, B), divq(A, B, 100, C).
eva2(X div Y, C) :- !, eva2(X, A), eva2(Y, B), divo(A, B, 100, C, _).
eva2(X mod Y, C) :- !, eva2(X, A), eva2(Y, B), divo(A, B, 100, _, C).
eva2(X++Y, C) :- !, eva2((X+Y) mod 360, C).
eva2(X--Y, C) :- !, eva2((X-Y) mod 360, C).
eva2(X^Y, C) :- !, eva2(X, A), eva2(Y, B), powq(A, B, 100, C).
eva2(sqrt(Y), C) :- !, eva2(Y, B), nthq(2, B, 100, C).
eva2(pi,ratnum(+,[355],[113])) :- !.
eva2(log(X,Y),C) :- !, eva2(X, A), eva2(Y, B), logq(A, B, 100, C).
eva2(gcd(X,Y),C) :- !, eva2(X, A), eva2(Y, B), gcdq(A, B, 100, C).
eva2(fix(X), C) :- !, eva2(X, A), intq(A, 100, C).
eva2(sin(X), C) :- !, eva2(X, A), sineval(A, C).
eva2(cos(X), C) :- !, eva2(X, A), coseval(A, C).
eva2(tan(X), C) :- !, eva2(X, A), taneval(A, C).
eva2(arcsin(X),C):- !, eva2(X, A), arcsineval(A, C).
eva2(arccos(X),C):- !, eva2(X, A), arccoseval(A, C).
eva2(arctan(X),C):- !, eva2(X, A), arctaneval(A, C).
eva2(abs(X), ratnum(+,N, D )) :- !, eva2(X, A), A = ratnum(_,N,D).
eva2(numer(X), ratnum(+,N,[1])) :- !, eva2(X, A), A = ratnum(_,N,_).
eva2(denom(X), ratnum(+,D,[1])) :- !, eva2(X, A), A = ratnum(_,_,D).
eva2(sign(X), ratnum(S,B,[1])) :- !, eva2(X, A), A = ratnum(S,N,_),
(N=[], B=[]; B=[1]), !.
eva2(xwd(X,Y),C) :- !, U is abs(Y)>>9, V is abs(Y)/\511,
eva2((X*512+U)*512+V, C).
eva2(X, C) :- ratnum(X, 100, S, N, D), !, C = ratnum(S, N, D).
eva2(Term, C) :- Term =.. [F,X,Y], relational_op(F,R, Flag),
!, eva2(X, A), eva2(Y, B), comq(A, B, 100, S), !,
combine_ops(R, S, Flag, C).
relational_op( =, =, true).
relational_op( \=, =, false).
relational_op( <, <, true).
relational_op( >=, <, false).
relational_op( >, >, true).
relational_op( =<, >, false).
relational_op(=:=, =, true).
relational_op(=\=, =, false).
combine_ops(Sign, Sign, Flag, Ans) :- !, Ans = Flag.
combine_ops(_, _, true,false) :- !.
combine_ops(_, _, false,true) :- !.
% Pretty-Print a number.
% This now always forces parentheses. When a
% proper general portray handler is written
% this could be made cleverer (as it once was).
% The magic numbers are 40 = "(", 41 = ")",
% 45 = "-", 47 = "/", 48 = "0" {ASCII codes}.
portray_ratnum(A) :-
ratnum(A, 100, S, N, D), !,
portray_ratnum(S, N, D).
portray_ratnum(_,[], []) :- !, % 0/0 = undefined
write(undefined).
portray_ratnum(+, _, []) :- !, % +N/0 = +infinity
write(infinity).
portray_ratnum(-, _, []) :- !, % -N/0 = -infinity
write(neginfinity).
portray_ratnum(+, N, [1]) :- !, % +N/1 = a +ve integer
putn(N).
portray_ratnum(-, N, [1]) :- !, % -N/1 = a -ve integer
put(45), putn(N).
portray_ratnum(+, N, D ) :- !, % +N/D = a +ve rational
put(40), putn(N), put(47), putn(D), put(41).
portray_ratnum(-, N, D ) :- !, % -N/D = a -ve rational
put(40), put(45), putn(N), put(47), putn(D), put(41).
putn([] ) :- !, put(48).
putn([D] ) :- !, write(D).
putn([D|T]) :-
putn(T),
/* If the Radix changes from 100 you *must*
add or delete lines here
*/
D1 is (D div 10) mod 10 +48, put(D1), % D1*10^1 +
D0 is (D) mod 10 +48, put(D0). % D0*10^0 = D.
%% INTERFACE CONVERSIONS %%
% Conversion of a number, of any form, to its
% essential bits.
ratnum(infinity, _, +,[1], []) :- !.
ratnum(neginfinity, _, -,[1], []) :- !.
ratnum(undefined, _, +, [], []) :- !.
ratnum(ratnum(S, N, D), _, S, N, D) :- !.
ratnum(N, R, +, L, [1]) :- integer(N), N >= 0, !, binrad(N, R, L).
ratnum(N, R, -, L, [1]) :- integer(N), N < 0, !, M is -N, binrad(M, R, L).
binrad(0, _, []) :- !.
binrad(N, R, [M|T]) :- K is N div R, M is N mod R, !, binrad(K, R, T).
% Normalise a number
standardise(ratnum(S,[N],[1]), Ans) :- !,
( S = '+', Ans = N
; S = '-', Ans is -N
), !.
standardise(ratnum(_, [],[1]), 0 ) :- !.
standardise(ratnum(S, N, []), Ans) :- !,
( N = [], Ans = undefined
; S = '+', Ans = infinity
; S = '-', Ans = neginfinity
), !.
standardise(Number, Number).
%% LOW LEVEL INTERFACE %%
% These routines provide a low level interface
% for procedures which want to operate directly
% on pairs of numbers.
% Only currently used by TIDY (27/2/81),
% so only those necessary are provided.
add(A, B, C) :- % eval(C is A+B).
addq(A, B, 100, X),
standardise(X, C).
multiply(A, B, C) :- % eval(C is A*B).
mulq(A, B, 100, X),
standardise(X, C).
power(A, N, C) :- % eval(C is A^B).
powq(A, N, 100, X),
standardise(X, C).
%% BASIC ARITHMETIC OVER RATIONALS %%
% Integer part of a rational
intq(A, R, ratnum(S, Q, [1])) :-
ratnum(A, R, S, N, D),
divn(N, D, R, Q, _).
% The greatest common divisor of two numbers is
% defined for all pairs of non-zero rationals.
% gcd(X,Y) = Z iff Z > 0 and there are integers
% M,N relatively prime for which X=MZ & Y=NZ.
gcdq(A, B, R, ratnum(+,Nd,Dd)) :-
ratnum(A, R, _, Na, Da),
ratnum(B, R, _, Nb, Db),
gcdn(Da, Db, R, _, Ga, Gb),
muln(Gb, Na, R, Ma),
muln(Ga, Nb, R, Mb),
gcdn(Ma, Mb, R, Nd),
muln(Gb, Da, R, Dd).
/* The above seems to be right, but I'm not sure. This IS right.
gcdq(A, B, R, ratnum(+,Nd,Dd)) :-
ratnum(A, R, _, Na, Da), % |A| = Na/Da
ratnum(B, R, _, Nb, Db), % |B| = Nb/Db
muln(Na, Db, R, N1), % N1 = Na.Db
muln(Nb, Da, R, N2), % N2 = Nb.Da
gcdn(N1, N2, R, Nc), % Nc = gcd(Na.Db, Nb.Da)
muln(Da, Db, R, Dc), % Dc = Da.Db
gcdn(Nc, Dc, R, _, Nd, Dd). % Nd/Dd = Nc/Dc in standard form
*/
/* % Take the inverse of a rational
invq(A, R, ratnum(S, D, N)) :-
ratnum(A, R, S, N, D).
*/
% Multiplication of two rationals
mulq(A, B, R, ratnum(Sc, Nc, Dc)) :-
ratnum(A, R, Sa, Na, Da),
ratnum(B, R, Sb, Nb, Db),
sign(Sa, Sb, Sc),
gcdn(Na, Db, R, _, Na1, Db1),
gcdn(Da, Nb, R, _, Da1, Nb1),
muln(Na1, Nb1, R, Nc),
muln(Da1, Db1, R, Dc).
% Division of two rationals
divq(A, B, R, ratnum(Sc, Nc, Dc)) :-
ratnum(A, R, Sa, Na, Da),
ratnum(B, R, Sb, Nb, Db),
sign(Sa, Sb, Sc),
gcdn(Na, Nb, R, _, Na1, Nb1),
gcdn(Da, Db, R, _, Da1, Db1),
muln(Na1, Db1, R, Nc),
muln(Da1, Nb1, R, Dc).
% Quotient and remainder of two rationals
divo(A, B, R, ratnum(Sq,Nq,[1]), ratnum(Sx,Nx,Dx)) :-
ratnum(A, R, Sa, Na, Da), % A = Sa.Na/Da
ratnum(B, R, Sb, Nb, Db), % B = Sb.Nb/Db
muln(Na, Db, R, N1), % A/B = (Sa.Na.Db)/(Sb.Nb.Da)
muln(Nb, Da, R, D1), % = (Sa.N1)/(Sb.D1)
divz(Sa,N1, Sb,D1, R, Sq,Nq, Sx,Ny),
muln(Da, Db, R, Dy), % A/B = Q + (Sx.Ny)/(Sb.Nb.Da)
gcdn(Ny, Dy, R, _, Nx, Dx). % A = Q.B + (Sx.Ny)/Dy
% Exponentiation of rationals
% This is always defined for (positive or
% negative) integer powers, however there
% is a current implementation restiction that
% the power be between -99999 and +99999 (ie
% within the current Radix).
% This may be defined for some rational powers
% but since there are results from this which are
% not representable as rationals it will fail
% in such cases. The code for rational powers
% relies on numerator and denominator being
% relatively prime, which is standard.
powq(A, B, R, C) :-
ratnum(B, R, S, N, [1]), !,
powq(S, N, A, R, C).
powq(A, B, R, C) :-
ratnum(B, R, S, N, [D]),
nthq(D, A, R, X), !,
powq(S, N, X, R, C).
powq(_, [], _, _, ratnum(+,[1],[1])) :- !.
powq(+,[N], A, R, ratnum(Sc, Nc, Dc)) :- !,
ratnum(A, R, Sa, Na, Da),
powz(Sa, Na, N, R, Sc, Nc),
pown(N, Da,[1],R, Dc).
powq(-,[N], A, R, ratnum(Sc, Nc, Dc)) :- !,
ratnum(A, R, Sa, Na, Da),
powz(Sa, Da, N, R, Sc, Nc),
pown(N, Na,[1],R, Dc).
% Negate a rational
negq(A, R, ratnum(Sc, Nc, Dc)) :-
ratnum(A, R, Sa, Nc, Dc),
( Nc = [], Dc = [], Sc = + % -undefined=undefined
; sign(Sa, -, Sc) % -0 = -(0) now.
), !.
% Addition of two rationals
addq(A, B, R, ratnum(Sc, Nc, Dc)) :-
ratnum(A, R, Sa, Na, Da),
ratnum(B, R, Sb, Nb, Db),
muln(Na, Db, R, Xa),
muln(Nb, Da, R, Xb),
addz(Sa,Xa, Sb,Xb, R, Sc,Xc),
gcdn(Xc, Da, R, _, Nx, Ya),
gcdn(Nx, Db, R, _, Nc, Yb),
muln(Ya, Yb, R, Dc), /*Q'*/ Nc/Dc\==[]/[], !.
addq(A, B, R, ratnum(Sc, Nc, [])) :- /*Q'*/
ratnum(A, R, Sa, Na, _),
ratnum(B, R, Sb, Nb, _),
( Na\==[], Nb\==[], Sa==Sb, Sc=Sa, Nc=[1]
; Sc= +, Nc=[]
), !.
% Subtraction of two rationals
subq(A, B, R, ratnum(Sc, Nc, Dc)) :-
ratnum(A, R, Sa, Na, Da),
ratnum(B, R, Sb, Nb, Db),
muln(Na, Db, R, Xa),
muln(Nb, Da, R, Xb),
subz(Sa,Xa, Sb,Xb, R, Sc,Xc),
gcdn(Xc, Da, R, _, Nx, Ya),
gcdn(Nx, Db, R, _, Nc, Yb),
muln(Ya, Yb, R, Dc), /*Q'*/ Nc/Dc\==[]/[], !.
subq(A, B, R, ratnum(Sc, Nc, [])) :- /*Q'*/
ratnum(A, R, Sa, Na, _),
ratnum(B, R, Sb, Nb, _),
( Na\==[], Nb\==[], Sa\==Sb, Sc=Sa, Nc=[1]
; Sc= +, Nc=[]
), !.
% Comparison of two rationals
comq(A, B, R, S) :-
ratnum(A, R, Sa, Na, Da), /*Q'*/ Na/Da \== []/[],
ratnum(B, R, Sb, Nb, Db), /*Q'*/ Nb/Db \== []/[],
muln(Na, Db, R, Xa),
muln(Nb, Da, R, Xb), !,
comz(Sa, Xa, Sb, Xb, S).
% Try to find Nth root
% This will fail in cases where the solution is
% not representable as a rational
nthq(N, A, R, ratnum(+, Nr, Dr)) :-
ratnum(A, R, +, Na, Da), !,
nthn(N, Na, R, Nr),
nthn(N, Da, R, Dr).
nthq(N, A, R, ratnum(-, Nr, Dr)) :-
ratnum(A, R, -, Na, Da), !,
1 is N mod 2,
nthn(N, Na, R, Nr),
nthn(N, Da, R, Dr).
nthn(_, [], _, []) :- !.
nthn(_, [1], _, [1]) :- !.
nthn(N, A, R, S) :-
newton(N, A, A, R, S), !,
pown(N, S, [1], R, B), !, B=A. % check that S^N=A !
newton(N, A, E, R, S) :-
M is N-1,
pown(M, E, [1], R, E1), % E1=E^(N-1)
mul1(E1,N, R, D2), % D2=N.E^(N-1)
muln(E, E1,R, E2), % E2=E^N
mul1(E2,M, R, N1), % N1=(N-1).E^N
addn(N1,A, 0, R, N2), % N2=(N-1).E^N+A
divn(N2,D2,R, F, _), % F = {(N-1).E^N+A}/{N.E^(N-1)}
comn(F, E, =, Z), !, % F Z E
newton(Z, N, A, F, R, S).
newton(<, N, A, F, R, S) :- !, newton(N, A, F, R, S).
newton(=, _, _, F, _, F) :- !.
% Take the logarithm of a rational to a rational base.
% This can be expected to fail for almost every pair
% of rational numbers. To keep the search space within
% logq(B, X, R, L) is true iff
% B, X, and L are rationals such that B^L = X.
% This does its best for strange mixtures, like log(-3,-27) = 3.
logq(B, X, R, L) :-
absq(B, S, C), % B S 0 & |B| = C
absq(X, T, Y), % X T 0 & |X| = Y
logq(S, T, C, Y, R, L).
% absq(A, R, S, B) is true iff
% A and B are rationals, |A| = B, and
% S = {+,-,0,*} as A {<,=,>} 0 or is undefined.
absq(ratnum(_,[],[]), *, ratnum(+,[],[])) :- !.
absq(ratnum(_,[],_), 0, ratnum(+,[],[1])) :- !.
absq(ratnum(Sa,Na,Da), Sa, ratnum(+,Na,Da)).
% logq(S, T, ...) is just a case analysis of logq.
logq(+, +, B, X, R, L) :- !,
ratlog(B, X, R, L).
logq(-, +, B, X, R, L) :- !,
ratlog(B, X, R, L),
mod2(L, 0). % L must be "even"
logq(-, -, B, X, R, L) :- !,
ratlog(B, X, R, L), !,
mod2(L, 1). % L must be "odd"
logq(+, -, _, _, _, ratnum(+,[],[])) :- !.
logq(*, _, _, _, _, ratnum(+,[],[])) :- !.
logq(_, *, _, _, _, ratnum(+,[],[])) :- !.
logq(0, _, _, _, _, ratnum(+,[],[])) :- !.
logq(_, 0, B, _, _, ratnum(Z, N,[])) :- !,
oneq(B, S, _),
logq(S, Z,N).
logq(+, -,[1]) :- !. % log(B,0) = -inf for 1
logq(-, +,[1]) :- !. % log(B,0) = +inf for 0
logq(_, +,[]). % log(B,0) = ???? otherwise
% oneq(A, S, B) is true when A and B are positive
% defined rationals, |log A| = log B, and S = sign(log A).
oneq(ratnum(_, _,[]), *, ratnum(+,[1],[])) :- !.
oneq(ratnum(_,Na,Na), 0, ratnum(+,Na,Na)) :- !.
oneq(ratnum(_,Na,Da), +, ratnum(+,Na,Da)) :-
comn(Na, Da, =, >), !.
oneq(ratnum(_,Na,Da), -, ratnum(+,Da,Na)).
% ratlog(B, X, R, L) is true iff
% B, X > 0 and B^L = X.
ratlog(B, X, R, L) :-
oneq(B, S, C), % B S 1 & |log B| = log C
oneq(X, T, Y),!,% X T 1 & |log X| = log Y
ratlog(S, T, C, Y, R, L).
% ratlog(S,T, ...) is just a case analysis
ratlog(+, +, B, X, R, ratnum(+,N,D)) :- !,
lograt(B, X, R, N, D).
ratlog(+, -, B, X, R, ratnum(-,N,D)) :- !,
lograt(B, X, R, N, D).
ratlog(-, +, B, X, R, ratnum(-,N,D)) :- !,
lograt(B, X, R, N, D).
ratlog(-, -, B, X, R, ratnum(+,N,D)) :- !,
lograt(B, X, R, N, D).
ratlog(0, _, _, _, _, ratnum(+,[], [])) :- !.
ratlog(_, 0, _, _, _, ratnum(+,[],[1])) :- !.
ratlog(+, *, _, _, _, ratnum(+,[1],[])) :- !.
ratlog(-, *, _, _, _, ratnum(-,[1],[])) :- !.
ratlog(_, *, _, _, _, ratnum(+,[], [])) :- !.
ratlog(*, _, _, _, _, ratnum(+,[], [])) :- !.
% lograt(B, X, R, N, D) is true iff
% B > 1, X > 1 are rationals, B^N = X^D, and gcd(N,D) = 1.
lograt(ratnum(+,Nb,Db), ratnum(+,Nx,Dx), R, [N], [D] ) :-
gcdn(Db, Nx, R, U), !, U = [1], % Db co-prime Nx
gcdn(Nb, Dx, R, V), !, V = [1], % Nb co-prime Dx
loop(Nb, Nx, R, G), !,
logn(G, 1, G, Nb, R, D), !, % D=log(G,Nb)
logn(G, 1, G, Nx, R, N), !, % N=log(G,Nx)
pown(N, Db, [1], R, K1),
pown(D, Dx, [1], R, K2), !,
K1 = K2. % Db^N = Dx^D
loop(A, B, R, G) :-
comn(A, B, =, S), !,
loop(S, A, B, R, G).
loop(=, A, _, _, A) :- !.
loop(<, A, B, R, G) :-
divn(B, A, R, Q, X), X = [], !,
loop(A, Q, R, G).
loop(>, A, B, R, G) :-
divn(A, B, R, Q, X), X = [], !,
loop(Q, B, R, G).
% logn(B, N, P, X, R, L) is true iff
% X >= B > 1, P = B^N, and X = B^L.
logn(_, N, X, X, _, N) :- !.
logn(B, N, P, X, R, L) :-
comn(P, X, =, <),
muln(B, P, R, Q),
M is N+1, !,
logn(B, M, Q, X, R, L).
%% BASIC ARITHMETIC OVER LONG INTEGERS %%
% Addition of two long integers
addz(+,A, +,B, R, +,C) :- !, addn(A, B, 0, R, C).
addz(+,A, -,B, R, S,C) :- !, subn(A, B, R, S, C).
addz(-,A, +,B, R, S,C) :- !, subn(B, A, R, S, C).
addz(-,A, -,B, R, -,C) :- !, addn(B, A, 0, R, C).
addn([D1|T1], [D2|T2], Cin, R, [D3|T3]) :-
Sum is D1+D2+Cin,
AbsSum is abs(Sum),
( AbsSum >= R, Cout = 1, D3 is AbsSum-R
; AbsSum < R, Cout = 0, D3 = Sum
), !,
addn(T1, T2, Cout, R, T3).
addn([], L, 0, _, L) :- !.
addn([], L, 1, R, M) :- !, add1(L, R, M).
addn(L, [], 0, _, L) :- !.
addn(L, [], 1, R, M) :- !, add1(L, R, M).
add1([M|T], R, [N|T]) :- N is M+1, N < R, !.
add1([M|T], R, [0|S]) :- R is M+1, !, add1(T, R, S).
add1([], _, [1]).
% Comparison of two long integers
comz(_,[],_,[],S) :- !, S = '='. % -0 = 0 now, alas.
comz(+,A, +,B, S) :- !, comn(A, B, =, S).
comz(+,_, -,_, >).
comz(-,_, +,_, <).
comz(-,A, -,B, S) :- !, comn(B, A, =, S).
comn([D1|T1], [D2|T2], D, S) :-
com1(D1, D2, D, N), !,
comn(T1, T2, N, S).
comn([], [], D, S) :- !, S = D.
comn([], _, _, <) :- !.
comn(_, [], _, >) :- !.
com1(X, X, D, D) :- !.
com1(X, Y, _, <) :- X < Y, !.
com1(X, Y, _, >) :- X > Y, !.
% Subtraction of two long integers
subz(+,A, +,B, R, S,C) :- !, subn(A, B, R, S, C).
subz(+,A, -,B, R, +,C) :- !, addn(A, B, 0, R, C).
subz(-,A, +,B, R, -,C) :- !, addn(B, A, 0, R, C).
subz(-,A, -,B, R, S,C) :- !, subn(B, A, R, S, C).
subn(A, B, R, S, C) :-
comn(A, B, =, O), !, % Oh for Ordering
subn(O, A, B, R, S, C).
subn(<, A, B, R, -, C) :- !, subp(B, A, 0, R, D), prune(D, C).
subn(>, A, B, R, +, C) :- !, subp(A, B, 0, R, D), prune(D, C).
subn(=, _, _, _, +,[]) :- !.
prune([0|L], M ) :- !,
prune(L, T),
(T = [], M = []; M = [0|T]).
prune([D|L], [D|M]) :- !,
prune(L, M).
prune([], []) :- !.
subp([D1|T1], [D2|T2], Bin, R, [D3|T3]) :-
S is D1-D2-Bin,
( S >= 0, Bout = 0, D3 = S
; S < 0, Bout = 1, D3 is S+R
), !,
subp(T1, T2, Bout, R, T3).
subp(L, [], 0, _, L) :- !.
subp(L, [], 1, R, M) :- !, sub1(L, R, M).
sub1([0|T], R, [K|S]) :- !, K is R-1, sub1(T, R, S).
sub1([N|T], _, [M|T]) :- M is N-1.
% Multiplication of Signs
sign(S, S, +) :- !.
sign(_, _, -) :- !.
% Multiplication of two long integers
/*
mulz(S,A, T,B, R, U,C) :-
sign(S, T, U), !,
muln(A, B, R, C).
*/
muln([], _, _, []) :- !.
muln(_, [], _, []) :- !.
muln(A, B, R, C) :- !, muln(A, B, [], R, C).
muln([D1|T1], N2, Ac, R, [D3|Pr]) :-
mul1(N2, D1, R, P2),
addn(Ac, P2, 0, R, Sm),
conn(D3, An, Sm), !,
muln(T1, N2, An,R, Pr).
muln([], _, Ac, _, Ac) :- !.
mul1(_, 0, _, []) :- !.
mul1(A, M, R, Pr) :- !,
mul1(A, M, 0, R, Pr).
mul1([], _, 0, _, []) :- !.
mul1([], _, C, _, [C]) :- !.
mul1([D1|T1], M, C, R, [D2|T2]) :-
D2 is (D1*M+C) mod R,
Co is (D1*M+C) div R,
mul1(T1, M, Co, R, T2).
% Exponentiation of a long integer to a short
% (Prolog) integer. Note that this means the
% power must be less than 100 (current radix).
% This code should always be called with positive
% powers.
powz(-,A, N, R, -,C) :-
N mod 2 =:= 1, !,
pown(N, A, [1], R, C).
powz(_,A, N, R, +,C) :- !,
pown(N, A, [1], R, C).
pown(0, _, M, _, M) :- !.
pown(1, A, M, R, P) :- !,
muln(A, M, R, P).
pown(N, A, M, R, P) :-
N1 is N div 2,
( N mod 2 =:= 0, M1 = M
; N mod 2 =:= 1, muln(A, M, R, M1)
),
muln(A, A, R, A1), !,
pown(N1, A1, M1, R, P).
% Division of two long integers
divz(S,A, T,B, R, U,Q, S,X) :-
sign(S, T, U), !,
divn(A, B, R, Q, X).
divn(_, [], _, _, _) :- !, fail. % division by 0 is undefined
divn(A,[1], _, A,[]) :- !. % a very common special case
divn(A,[B], R, Q, X) :- !, % nearly as common a case
div1(A, B, R, Q, Y),
conn(Y, [], X).
divn(A, B, _, Q, X) :-
comn(A, B, =, S),
( S = '<', Q = [], X = A
; S = '=', Q = [1], X = []
), !.
divn(A, B, R, Q, X) :- !,
divm(A, B, R, Q, X).
conn(0, [], []) :- !.
conn(D, T, [D|T]).
div1([D1|T1], B1, R, Q1, X1) :- !,
div1(T1, B1, R, Q2, X2),
D2 is (X2*R+D1) div B1,
X1 is (X2*R+D1) mod B1,
conn(D2, Q2, Q1).
div1([], _, _, [], 0).
% divm(A, B, R, Q, X) is called with A > B > R
divm([D1|T1], B, R, Q1, X1) :- !,
divm(T1, B, R, Q2, X2),
conn(D1, X2, T2),
div2(T2, B, R, D2, X1),
conn(D2, Q2, Q1).
divm([], _, _, [], []).
div2(A, B, R, Q, X) :-
estd(A, B, R, E), !,
chkd(A, B, R, E, 0, Q, P), !,
subn(A, P, R, _, X). % S=+
div2(A, B, _, _, _) :-
long_error(divq, A/B).
estd([_,A1,A2], [_,B1], R, E) :-
B1 >= R/2, !,
E is (A2*R+A1) div B1.
estd([A0,A1,A2], [B0,B1], R, E) :- !,
L is (A2*R+A1) div (B1+1),
mul1([B0,B1], L, R, P),
subn([A0,A1,A2], P, R, _, N), !, %_=+
estd(N, [B0,B1], R, M), !,
E is L+M.
estd([A0,A1], [B0,B1], R, E) :- !,
E is (A1*R+A0+1) div (B1*R+B0).
estd([_], _, _, 0) :- !.
estd([_|Ar], [_|Br], R, E) :- !,
estd(Ar, Br, R, E).
estd([], _, _, 0) :- !.
chkd(A, B, _, _, 3, _, _) :- !,
long_error(divq, A/B).
chkd(A, B, R, E, _, E, P) :-
mul1(B, E, R, P),
comn(P, A, <, <), !.
chkd(A, B, R, E, K, Q, P) :-
L is K+1, F is E-1, !,
chkd(A, B, R, F, L, Q, P).
% GCD of two long integers
/*
gcdz(S,A, T,B, R, D, S,M, T,N) :- !,
gcdn(A, B, R, D, M, N).
*/
gcdn([], [], _, [1], [], []) :- !.
gcdn([], B, _, B, [], [1]) :- !.
gcdn( A, [], _, A, [1], []) :- !.
gcdn([1], B, _, [1], [1], B) :- !. % common case
gcdn( A,[1], _, [1], A, [1]) :- !. % common case
gcdn( A, B, R, D, M, N) :- % A, B > 1
gcdn(A, B, R, D),
divn(A, D, R, M, _),
divn(B, D, R, N, _).
gcdn(A, B, R, D) :- % A, B >= 1 !!
comn(A, B, =, S), !,
gcdn(S, A, B, R, D).
gcdn(<,[], B, _, B) :- !.
gcdn(<, A, B, R, D) :-
estg(B, A, R, E),
muln(E, A, R, P),
subn(B, P, R, _, M), !,
gcdn(A, M, R, D).
gcdn(>, A,[], _, A) :- !.
gcdn(>, A, B, R, D) :-
estg(A, B, R, E),
muln(E, B, R, P),
subn(A, P, R, _, M), !,
gcdn(M, B, R, D).
gcdn(=, A, _, _, A).
estg( A, [B], R, E) :- !,
div1(A, B, R, Q, X),
( X*2 =< B, E = Q
; add1(Q, R, E)
), !.
estg([_|A], [_|B], R, E) :- !,
estg(A, B, R, E).
%% TRIGONOMETRIC EVALUATION %%
% This stuff needs some work done on it, and the mode
% declarations haven't been written yet. Taihoa.
% To do:
% Since at this stage all the argumentss are known to be
% numbers we shouldn't waste time using the general eval.
% Approximations should be used so that the routines work
% for ANY argument. Care is needed, since little is known
% about rational approximations, lest the numbers explode.
sineval(X, S) :-
eval(X < 0), !,
eva2(-X, Y),
sineval(Y, T),
eva2(-T, S).
sineval(X, S) :-
eval(X > 90), !,
eva2(180-X, Y),
sineval(Y, S).
sineval(X, S) :- % 0 <= X <= 90
sineval1(X, S).
sineval1(ratnum(+,[],[1]), ratnum(+,[],[1])).
sineval1(ratnum(+,[30],[1]), ratnum(+,[1],[2])).
sineval1(ratnum(+,[45],[1]), ratnum(+,[99],[140])).
sineval1(ratnum(+,[60],[1]), ratnum(+,[45],[52])).
sineval1(ratnum(+,[90],[1]), ratnum(+,[1],[1])).
coseval(X, C) :-
eva2(90-X, Y), !,
sineval(Y, C).
taneval(X, T) :-
sineval(X, S),
coseval(X, C), !,
eva2(S/C, T).
arcsineval(S, X) :-
eval(S >= 0), !,
sineval1(X, S).
arcsineval(S, X) :-
eval(S < 0),
eva2(-S, T), !,
sineval1(Y, T),
eva2(-Y, X).
arccoseval(C, X) :-
arcsineval(C, Y), !,
eva2(90-Y, X).
arctaneval(ratnum(S,N,D), ratnum(S,M,C)) :-
arctaneval(N, D, M, C).
arctaneval([],X, [], X) :- !. % arctan(0) = 0, arctan(undef) = undef
arctaneval(X, X, [45], [1]) :- !. % arctan(1) = 45`
arctaneval(_,[], [90], [1]) :- !. % arctan(inf) = 90`
arctaneval(N, D, M, C) :-
R = 100, % the common radix
muln(N, N, R, Nsq),
muln(D, D, R, Dsq),
addn(Nsq, Dsq, 0, R, Sq), !,
nthn(2, Sq, R, Den), !,
arcsineval(ratnum(+,N,Den), S),
S = ratnum(+,M,C).
%% ERROR HANDLING %%
long_error(Culprit, Expression) :-
long_error_message(Culprit, Message),
display('** '), display(Message), display(': '),
print(Expression), ttynl,
break, fail.
long_error_message(eval, 'EVAL given a variable').
long_error_message(eva2, 'EVAL given an expression containing a variable').
long_error_message(divq, 'Unexpected rational division problem').
ttynl:-nl.
break.