/* 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), !. 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').