% Pure, declarative, and constructive binary arithmetics % % aka: Addition, Multiplication, Division with the remainder, % Discrete Logarithm as sound and complete, pure and declarative % relations that can be used in any mode whatsoever. % The relations define arithmetics over base-2 non-negative numerals % of *arbitrary* size. % % In particular, we define the predicate div/4 such that div(N,M,Q,R) % succeeds if and only if N = M*Q + R and 0<=R0! % % Incidentally, the relation 'div' encompasses all four operations % of arithmetics: we can use div(X,Y,Z,0) to multiply and divide and % div(X,Y,1,R) to add and subtract. % % Our log/4 predicate has the property that log(N,B,Q,R) succeeds if % and only if N = B^Q + R where 0 <= R and Q is the largest such integer. % We can use log/4 to exponentiate, to find discrete logarithms, or % to find the logarithm base. % % % At first glance, our problem seems trivial. We can use the mu-calculus % definition of, say, multiplication to write % mul(X,Y,Z) = mu(X1,Y1,Z1). X1*Y1=Z1 & X1=X & Y1=Y & Z1=Z % and then replace mu(X1,Y1,Z1) with a relation that generates all % possible triples of numerals. The problem comes when we want to evaluate % mul(two,Y,one). Obviously, there is no such Y that makes the relation % hold. The naive mu-calculus implementation will keep generating % numbers for Y forever, as none of the generated numbers will % make 2*Y=1 true. % % The naive mu-calculus approach may fail even if the relation has a solution. % For example, if mu(X1,Y1,Z1) is implemented as % gen(X1), gen(Y1), gen(Z1) % then evaluating mul(X,Y,2) will generate X1=0 first and then % will stuck in generating successive numbers for Y trying to % find one that makes 0*Y=2 hold. % % The present approach places the correct upper bounds on the % generated numbers to make sure the search process terminates. % Therefore, our arithmetic relations are not only sound % (e.g., if mul(X,Y,Z) holds then it is indeed X*Y=Z) but also % complete (if X*Y=Z is true then mul(X,Y,Z) holds) and % nearly refutationally complete (if X*Y=Z is false then mul(X,Y,Z) fails, % in finite time -- provided that X, Y and Z share no common % logic variables. % % We give two implementations of `mul', both of which have the % properties of soundness and nearly refutational completeness. % The first implementation is six times faster, but it does not % always recursively enumerate its domain if that domain is infinite. % This is the case when mul(X,Y,Z) is invoked when all three X, Y, and % Z being uninstantiated variables. The predicate in that case has % the infinite number of solutions, as expected. Alas, those solutions % look as follows: % X = 2, Y = 3, Z = 6 % X = 4, Y = 3, Z = 12 % X = 8, Y = 3, Z = 24 % X = 16, Y = 3, Z = 48 % That is, mul(X,Y,Z) keeps generating solutions where X is a power of % two. Therefore, when the solution set of the predicate `mul' is infinite, it % truly produces an infinite set of solutions -- but only a subset of % all possible solutions. In other words, `mul' does not recursively % enumerate the set of all numbers such that X*Y=Z if that set is infinite. % % Therefore, % X = [l,l], Y = [l,l], mul(X,Y,Z) % mul(X,Y,Z), X = [l,l], Y = [l,l] % work differently. The former terminates and binds Z to the representation % of 9 (the product of 3 and 3). The latter fails to terminate. % This is not generally surprising as 'commas' in Prolog are not truly % conjunctions: they are not commutative. Still, we would like our `mul' % to have the algebraic properties expected of multiplication. % % The second version of `mul', albeit slower, almost completely fixes % the problem. It recursively enumerates the domain of numbers X, Y, Z % such as X*Y=Z and Y <= X. A test at the end of this file shows that % if I and J are some definite numbers, then % X = I, Y = J, mul(X,Y,Z) % mul(X,Y,Z), (X = I, Y = J; X = J, Y = I) % are equivalent. This is true no matter which numbers I and J are. % The predicate mul(X,Y,Z) will, in finite number of steps, produce any % solution in its domain. We achieve that property _without_ giving % up refutational completeness: if X and Y and Z are so instantiated % that X*Y /= Z, then mul(X,Y,Z) fails -- in finite time! Incidentally, % this is the case when either X or Y are uninstantiated, as in % X=[o,l], Z=[l,l], mul(X,Y,Z). % because there can't be any number Y such that 2*Y=3. Also, % if Z is instantiated, then mul(X,Y,Z) delivers _all_ factorizations % of Z. % % We should note that conjunctions of arithmetic predicates can't be always % terminating. Ditto if arguments of predicates share common logic % variables. Indeed, mul([l,l],[l|X],[l,o|X]) will fail to terminate. % This is because mul([l,l],[l|X],[l,o|X]) is equivalent to % the conjunction of primitive goals % mul([l,l],N,P), add([l],X2,N), add([l],X4,P), add(X,X,X2), add(X2,X2,X4) % The mul predicate will produce the infinite stream of numbers % with the property P = 3N, but none of these numbers will satisfy % the constraint P = 4X+1, N = 2X+1. % Well, in the case above the impossibility is easy to see. However, that % case is just a particular case of a Diofantine equation. In general, % any Diofantine equation can be represented as a conjunction of % arithmetic functions. Refutational completeness of such a conjunction % is therefore equivalent to Hilbert's 10th problem, which is undecidable. % % Arithmetic predicates are easy to implement in an impure % system such as Prolog, with the help of the predicate 'var'. The latter % can tell if its argument is an uninstantiated variable. However, 'var' % is impure. The present file shows the implementation of arithmetic % relations in a _pure_ logic system -- no cuts, no 'var' predicate, % no negations. % % We should point out that in the Kanren system % http://kanren.sourceforge.net/ % it is possible to implement the recursively enumerating `mul' % predicate without the restriction Y <= X. This is because Kanren % offers a `fair-scheduling' disjunction: any-interleave. That implementation % of arithmetics is still pure and constructive. % % % The pure and constructive arithmetics (addition/subtraction) % on _decimal_ numbers is described in Section 6 of % http://pobox.com/~oleg/ftp/Haskell/number-parameterized-types.ps.gz % That paper implements the arithmetics in the type system. % Therefore, we can represent numerical equality and inequality constraints % in the type system of Haskell. We can statically type the function vappend % (that appends two length-number-parameterized vectors and makes sure % the size of the result is the sum of the sizes of the arguments -- % statically). We can statically make sure that the function tail is applied % to a non-empty vector. The implementation in the paper is assuredly pure % as the type system of Haskell has no 'cut', 'var', etc. predicates. % % % % The numerals are represented in the binary little-endian % (least-significant bit first) notation. The higher-order bit must be 1. % Zero bit is denoted by 'o' (low-case oh) and the set bit is denoted by % 'l' (low-case el). % [] represents 0 % [l] represents 1 % [o,l] represents 2 % [l,l] represents 3 % [o,o,l] represents 4 % etc. % Please see the tests at the very end. % % This file was tested with SWI-Prolog Version 5.4.5 % % $Id: pure-bin-arithm.prl,v 1.12 2007/06/11 01:12:34 oleg Exp oleg $ % Auxiliary predicate to build and show binary numerals. % The following predicate is not the part of the system, strictly speaking. % It is also `impure'. None of the other predicates need +/- annotations % because they can be used in any mode whatsoever, and none of the % other predicates use cut. % The auxiliary predicate is "input-output", to input and output numerals % from/into more convenient representation. I must say that % the binary notation is not too bad per se. ?- op(700,xfx,isb). % +PrologNum isb -OurNum % -PrologNum isb +OurNum 0 isb []. N isb L :- number(N), N > 0, !, N1 is N // 2, P is N mod 2, (P = 0, L = [o|L1], N1 isb L1; P = 1, L = [l|L1], N1 isb L1), !. N isb [o|L1] :- !, N1 isb L1, N is 2*N1. N isb [l|L1] :- !, N1 isb L1, N is 2*N1 + 1. %---------------------------------------------------------------------- % Addition and subtraction % zero: N holds if N is zero zero([]). % not a zero pos([_|_]). % at least two gtl([_,_|_]). % Generator of binary numerals, aka the type predicate genb([]). % 0 is a number genb([l|X]) :- genb(X). % if X is a number, so is 2X+1 genb([o|X]) :- pos(X), genb(X). % if X>0 is a number, so is 2X % One-bit full-adder: carry-in a b r carry-out % The relation holds if % carry-in + a + b = r + 2*carry-out % where carry-in a b r carry-out are all either o or l. full1_adder(o,o,o,o,o). full1_adder(o,l,o,l,o). full1_adder(o,o,l,l,o). full1_adder(o,l,l,o,l). full1_adder(l,o,o,l,o). full1_adder(l,l,o,o,l). full1_adder(l,o,l,o,l). full1_adder(l,l,l,l,l). % Multiple-bit full-adder (ripple-carry adder): carry-in a b r % holds if carry-in + a + b = r % where a, b, and r are binary numerals and carry-in is either o or l fulln_adder(o,A,[],A). % 0 + a + 0 = a fulln_adder(o,[],B,B) :- pos(B). % 0 + 0 + b = b fulln_adder(l,A,[],R) :- fulln_adder(o,A,[l],R). % 1 + a + 0 = 0 + a + 1 fulln_adder(l,[],B,R) :- pos(B), fulln_adder(o,[l],B,R). % 1 +0 + b = 0 + 1 + b % The following three clauses are needed to make all numbers % well-formed by construction, that is, to make sure the % higher-order bit is one. fulln_adder(Cin,[l],[l],R) :- % c + 1 + 1 >= 2, always two-bit R = [R1,R2], full1_adder(Cin,l,l,R1,R2). % c + 1 + (2*br + bb) = (2*rr + rb) where br > 0, and so is rr > 0 fulln_adder(Cin,[l],[BB|BR],[RB|RR]) :- pos(BR), pos(RR), full1_adder(Cin,l,BB,RB,Cout), fulln_adder(Cout,[],BR,RR). fulln_adder(Cin,A,[l],R) :- % symmetric case for the above gtl(A), % a > 1, the case of a=1 is handled already gtl(R), % and so is r > 1 fulln_adder(Cin,[l],A,R). % carry-in + (2*ar + ab) + (2*br + bb) % = (carry-in + ab + bb) (mod 2) % + 2*(ar + br + (carry-in + ab + bb)/2) % The cases of ar= 0 or br = 0 have already been handled. % So, now we require ar >0 and br>0. That implies that rr>0. fulln_adder(Cin,[AB|AR],[BB|BR],[RB|RR]) :- pos(AR), pos(BR), pos(RR), full1_adder(Cin,AB,BB,RB,Cout), fulln_adder(Cout,AR,BR,RR). % a + b = c add(A,B,C) :- fulln_adder(o,A,B,C). % a - b = c sub(A,B,C) :- add(B,C,A). % less: N M % holds if N < M, or % if exists X>0, N + X = M less(N,M) :- pos(X), add(N,X,M). % The above addition and subtraction are not recursively enumerating, % but can be made to, using the technique described below. %---------------------------------------------------------------------- % Multiplication and division % lessl a b % holds if a=0 and b>0, or if (floor (log2 a)) < (floor (log2 b)) % That is, we compare the length (logarithms) of two numerals % For a positive numeral, its bitlength = (floor (log2 n)) + 1 lessl([],[_|_]). lessl([_|X],[_|Y]) :- lessl(X,Y). % samel a b % holds if both a and b are zero % or if (floor (log2 a)) = (floor (log2 b)) samel([],[]). samel([_|X],[_|Y]) :- samel(X,Y). % lessl3 p1 p n m % holds iff % p1 = 0 and p > 0 or % length(p1) < min(length(p), length(n) + length(m) + 1) lessl3([],[_|_],_,_). lessl3([_|P1R],[_|PR],[],[_|MR]) :- lessl3(P1R,PR,[],MR). lessl3([_|P1R],[_|PR],[_|NR],M) :- lessl3(P1R,PR,NR,M). % n * m = p %% mul([],M,[]). % 0 * m = 0 %% mul(N,[],[]) :- pos(N). % n * 0 = 0, n > 0 %% mul([l],M,M) :- pos(M). % 1 * m = m, m > 0 % (2*nr) * m = 2*(nr*m), m>0 (the case of m=0 is taken care of already) % nr > 0, otherwise the number is ill-formed %% mul([o|NR],M,[o|PR]) :- %% pos(M), pos(NR), pos(PR), %% mul(NR,M,PR). % (2*nr+1) * m = 2*(n*m) + m % m > 0; also nr>0 for well-formedness % the result is certainly greater than 1. % we note that m > 0 and so 2*(nr*m) < 2*(nr*m) + m % and (floor (log2 (nr*m))) < (floor (log2 (2*(nr*m) + m))) % Furthermore, (floor (log2 (* n m))) + 1 <= % (floor (log2 n)) + 1 + (floor (log2 m)) + 1 %% mul([l|NR],M,P) :- %% pos(M), pos(NR), gtl(P), %% lessl3(P1,P,[l|NR],M), %% mul(NR,M,P1), %% add([o|P1],M,P). % Now, the last clause needs some explanation. % Relations lessl3 and mul are infinitary: their solution set can be infinite. % (or, in other words, they can be backtracked infinitely many times). % If 'add' on the last line keeps failing, we may become stuck. % First we note that if P is ground, then lessl3(P1,P,_,_) is % a finitary relation -- its solution set is finite. Therefore, if % mul(NR,M,P1) and add keep failing, we eventually exhaust all solutions to % lessl3(P1,P,_,_) and fail the 'mul' clause. % Let us make the claim more precise: % Proposition 1. If P is definite number, mul(N,M,P) has a finite % number of solutions. % Proof: by induction. Base case: P = 0. % Inductive case: proposition true for all P10 and prove by induction on N. The only non-trivial % case is N odd. Relation mul(NR,M,P1) is finite (semi-deterministic, % actually). Relation lessl3(P1,P,[l|NR],M) has the finite number % of solutions if both NR and M are definite numbers. Thus the last % clause of mul is terminating and is semi-deterministic % under conditions of Proposition 2. % % If both N and P (or M and P) are unbounded, the relation mul(N,M,P) % will have the infinite number of solutions. However, given an arbitrary % bound on the number of solutions, that number of solutions will be found % in finite time. The reason is that lessl3(P1,P,[l|NR],M) (where P and % either NR or M are unbound) generates (P1=k, P=k+1+n where n is unbound). % There is always a term in that sequence (a solution of lessl3) % that makes both mul(NR,M,P1) and add([o|P1],M,P) hold % for infinitely many M (or NR). % % Please note that the above proposition said: "given an arbitrary _bound_ % on the _number_ of solutions, that _number_ of solutions will be found % in finite time." The proposition did not say: given an arbitrary % _solution_. This is because the above `mul' relation is not recursively % enumerating, as was explained in the title comments. Given below % is a recursively enumerating version (with a restriction that doesn't % seem significant). Alas, it is also about six times slower. % Recursively-enumerating predicate mul(X,Y,Z). % It recursively enumerates its domain X*Y = Z, Y <= X. % Note that if X and Y are instantiated, then Y may still be greater % than X. So, the following mul() does behave like the mul above, % but in addition it enumerates all solutions such that X*Y = Z, Y <= X. % We still maintain the nearly refutational completeness. % % Please note the analogy with space-filling curves and the proof % that rational numbers are countable. % The key is to write such a predicate pred(N,M,P) that enumerates % all numerals such that M <= N. % The first, more complex attempt: % lessl31([],[_|_],_,_). % lessl31([_|P1R],[_|PR],N,[_|MR]) :- lessl31(P1R,PR,N,MR). % lessl31([_|P1R],[_|PR],[_|NR],M) :- lessl31(P1R,PR,NR,M). % % Triangular: % % lessl31(C,P,N,M), lessl(M,[_|C]), lessl(N,[_|C]). % lessl2(P,N,M) holds if length(P) <= length(N) + length(M) lessl2([],_,_). lessl2([_|PR],[],[_|MR]) :- lessl2(PR,[],MR). lessl2([_|PR],[_|NR],M) :- lessl2(PR,NR,M). % A better solution follows. % Let us consider the properties of the predicate lessl(N,P). % It generates solutions such as % N = [], P = [_G293|_G294] ; % N = [_G293], P = [_G296, _G299|_G300] ; % N = [_G293, _G299], P = [_G296, _G302, _G305|_G306] % Let us introduce definitions: % A variable X is called L-instantiated if it is instantiated % to a list of finite length, possibly with un-instantiated % elements % A variable is called L-uninstantiated if it is instantiated % to a list with an uninstantiated tail. % For example, [_,_] is L-instantiated and [_,_|_] is L-uninstantiated. % % Proposition lessl: all the solutions of the predicate lessl(N,P) % have N L-instantiated. OTH, P is at most is L-uninstantiated, with % progressively longer L-instantiated prefixes. % Proof is by trivial induction on the definition of lessl. % Corollary lessl: if P is L-instantiated, lessl(N,P) has % the finite number of solutions. % Corollary lessl2: if N is L-instantiated, lessl(N,P) has % only one solution. % Let us define enum(N,M,P) :- lessl(N,[_|P]), lessl(M,[_|N]), lessl2(P,N,M). % Here lessl(M,[_|N]) is a shorter way to write (lessl(M,N); samel(M,N)) % The relation enum(N,M,P) has the following properties. % % Proposition enum1: if P is L-instantiated, enum(N,M,P) has the % finite number of solutions. % Proof: if P is L-instantiated, lessl(N,[_|P]) has the finite number % of solutions, by Corollary lessl. Each of those solutions % have N L-instantiated. Therefore, lessl(M,[_|N]) % has the finite number of solutions. So does lessl2(P,N,M), % as is clear by inspection. % % Proposition enum2: if N and M are both L-instantiated, enum(N,M,P) % has the finite number of solutions. % If N is L-instantiated, lessl(N,[_|P]) has only one solution. P may still % be L-uninstantiated. However, lessl2(P,N,M) will have only finite % number of solutions if both N and M are L-instantiated. % % Proposition enum3: If N is L-instantiated, enum(N,M,P) has the % finite number of solutions. % Proof: if N is L-instantiated, lessl(N,[_|P]) has only one solution, % with P at most L-uninstantiated. lessl(M,[_|N]) will % have exactly N+1 solutions. In all these solutions, M will be % L-instantiated. Therefore, lessl2(P,N,M) will have the finite % number of solutions in that case as well. % % Proposition enum4: If M is L-instantiated, enum(N,M,P) recursively % enumerates its domain, of L-instantiated N, M and P so that % length(M) <= length(N) and length(P) <= length(N) + length(M). % Proof: If both N and P are uninstantiated, lessl(N,[_|P]) will generate % the infinite set, in which N is L-instantiated % and length(N) <= length(P). % The conjunction "lessl(N,[_|P]), lessl(M,[_|N])" will have only those % L-instantiated N such that length(M) <= length(N). % The predicate lessl2(P,N,M) then will produce all those L-instantiated % P such that length(P) <= length(N) + length(M). % % Proposition enum5: enum(N,M,P) recursively enumerates its domain. % Proof: similar to the above, noting that after lessl(N,[_|P]), % N will be L-instantiated, and so lessl(M,[_|N]) will produce all such % M so that length(M) <= length(N) -- and there is the finite number of % those, for each L-instantiated N. Similarly, for each L-instantiated % N and M, lessl2(P,N,M) produces all those L-instantiated P % so that length(P) <= length(N) + length(M) -- and again, there % is the finite number of those, for each L-instantiated N and M. mul([],_,[]). % 0 * m = 0 mul(N,[],[]) :- pos(N). % n * 0 = 0, n > 0 mul(N,[l],N) :- pos(N). % n * 1 = n, n > 0 % Now, for all the clauses below, we assert M > 1. We also know that N>0 % This implies that width(N) < width(P) mul([l],M,M) :- gtl(M). % 1 * m = m, m > 1 % The first clause below uses a version of the enum() predicate. % We know that N,M,P must have the width at least 1, and % width(N) < width(P). Note that we did not write lessl2() explicitly: % it is implicit in the predicate mul1() (see, e.g., lessl3 below). mul(N,M,P) :- gtl(M), gtl(N), lessl(N,P), lessl(M,[_|N]), mul1(N,M,P). % This clause below is to handle the case of mul(N,M,P) when % both N and M are instantiated and length(N) < length(M) -- % or when P is instantiated. % When either (N and P) or (M and P) or all three are uninstantiated, % the clause above starves the clause below. mul(M,N,P) :- gtl(M), gtl(N), lessl(N,P), lessl(M,N), mul1(N,M,P). % The following is just like the `mul' in the previous section. mul1([l],M,M) :- gtl(M). % 1 * m = m, m > 1 mul1([o|NR],M,[o|PR]) :- gtl(M), pos(NR), pos(PR), mul1(NR,M,PR). % (2*nr+1) * m = 2*(n*m) + m % m > 1; also nr>0 for well-formedness % the result is certainly greater than 1. % we note that m > 1 and so 2*(nr*m) < 2*(nr*m) + m % and (floor (log2 (nr*m))) < (floor (log2 (2*(nr*m) + m))) % Furthermore, (floor (log2 (* n m))) + 1 <= % (floor (log2 n)) + 1 + (floor (log2 m)) + 1 mul1([l|NR],M,P) :- gtl(M), pos(NR), gtl(P), lessl3(P1,P,[l|NR],M), mul1(NR,M,P1), add([o|P1],M,P). %------------------------------------------------------------------------ % Division % div n m q r holds iff % n = q*m + r % where 0<=r 2*m, then lessl(n,m) fails quickly, and we don't need to do % the expensive test less(n,m). % Obviously, 0<=r=M div1(N,M,[l],R) :- % N >= M samel(N,M), % N = 1*M + R, R = N - M add(R,M,N), less(R,M). div1(N,M,[],N) :- % N < M, R=N, Q=0 samel(N,M), % R < M, obviously. less(N,M). % Now, N has more digits than M, so Q >0. Therefore, Q*M <= N % and Q*M < 2*N and so (floor (log2 (Q*M))) < (floor (log2 (2*N))) div1(N,M,Q,R) :- lessl(M,N), less(R,M), lessl(P,[o|N]), add(P,R,N), mul(Q,M,P). % First note that div is non-recursive. % If N is a definite number and the first div clause does not apply then % div has the finite number of solutions. % In particular, the last clause has the finite number of solutions % because for a definite N, lessl(M,N), less(R,M), lessl(P,[o|N]), % and add(P,R,N) are all finite. % If N is not a definite number, then the last clause has the infinite % number of solutions -- but all of them are enumerable. % One may think that if we remove the test "lessl(M,N)" from the last % clause, then we don't need the three previous ones. If we did that % however, the following simple question (finding if dividing 5 can ever % give us 7) % 5 isb A5, 7 isb A7, div(A5,X,A7,R). % would diverge (rather than fail in finite time). This is because % when Q is a definite number greater than 1, M must be less than N. % If that condition is not asserted, relation less(R,M) will keep generating % bigger and bigger numbers. In other words, when N is definite, we % can no longer guarantee that div() will have the finite number of solutions. % % One may also think that we can get a simpler div() by a more rearrangement % of clauses: div3(N,M,Q,R) :- lessl(P,[o|N]), mul(Q,M,P), less(R,M), add(P,R,N). % This div3 too has the finite number of solutions whenever N is instantiated. % However, div3 does not finitely fails (rather it loops forever) % when we try divide by zero: div3(N,[],Q,[]). % The original div() doesn't have this drawback. % Here's another example: % N=[], div(N,M,Q,R) % it gives one solution, with the original, `complex' implementation of % div -- but not for any simple, one-clause div, such as div3. % However, all is not well % 7 isb Q, 5 isb M, div(N,M,Q,[]) % diverges with all divs so far, if we ask for more than 1 solution. % We wish that multiplication of two definite numbers, 7 and 5, gave one, % and only one result. So, we need a more complex div. % A faster and a more refutationally complete div algorithm % Again, div(N,M,Q,R) holds iff N = M*Q + R % Let l be the bit-length of r (if r=0, l=0). % Let n = 2^(l+1) * n1 + n2 % q = 2^(l+1) * q1 + q2 % Note that n1 or q1 may be zero. % We obtain that % n = m*q + r % is equivalent to the conjunction of the following two relations % q2*m + r - n2 is divisible by 2^(l+1) % n1 = q1*m + (q2*m + r - n2)/2^(l+1) % We note that by construction (see the mentioning of lessl(M,N) below) % all numbers in 'q2*m + r - n2' are length-limited. Therefore, we can % obtain the success or failure of 'q2*m + r - n2' in finite time. % This fact let us fail the 'div' relation, in finite time, % when it obviously does not hold, as in % div([0|X],[o,l],Q,[l]) % (because no even number can give the remainder 1 upon the % division by two). % % We should note that if n1=0, we obtain that q1 must be zero % (because m>0) and q2*m + r = n2. The latter can be solved in finite % time. % We also note that (q2*m + r - n2)/2^(l+1) < m % because r - n2 < (2^(l+1) - q2)* m % because 2^(l+1) - q2 >=1 and m > r by construction. Therefore, to % solve the relation n1 = q1*m + (q2*m + r - n2)/2^(l+1) we use % div itself: div(n1,m,q1,(q2*m + r - n2)/2^(l+1)) % Thus our division algorithm is recursive. On each stage we determine at % least one bit of the quotient (if r=0, l=0 and q2 is either 0 or 1), % in finite time. % Chung-chieh Shan has pointed out that the present algorithm is % akin to the elementary school long-division algorithm, only % done right-to-left. % the divisor M is bigger than the dividend N: Q=0,N=R div(N,M,[],R) :- N = R, less(N,M). % N is at least M and certainly has the same number of digits as M div(N,M,[l],R) :- samel(N,M), add(R,M,N), less(R,M). % N has more digits than the divisor M. div(N,M,Q,R) :- lessl(M,N), % By the property of lessl/2, M has become L-instantiated % (if it wasn't already) less(R,M), % By the property of less/2, R is now L-instantiated pos(Q), % Q must be positive then split(N,R,N1,N2), split(Q,R,Q1,Q2), ( N1 = [], Q1 = [], % q2*m + r = n2 sub(N2,R,Q2M), mul(Q2,M,Q2M) % provably terminates ; pos(N1), mul(Q2, M, Q2M), add(Q2M,R, Q2MR), sub(Q2MR, N2, RR), % rr = q2*m + r - n2 split(RR, R, R1, []), % r1 = rr/2^(l+1), evenly div(N1, M, Q1, R1) ). % split N R N1 N2 % holds if n = 2^(l+1)*n1 + n2 where l = bitlength(r) % This relation makes sense to use only when 'r' is L-instantiated. % In that case, the relation has only the finite number of solutions, in % all of which n2 is L-instantiated. % We take trouble to assure that we produce only well-formed numbers: % the major bit must be one. split([], _, [], []). split([o,B|N], [], [B|N], []). split([l|N], [], N, [l]). split([o,B|N], [_|R], N1, []) :- split([B|N], R, N1, []). split([l|N], [_|R], N1, [l]) :- split(N, R, N1, []). split([B|N], [_|R], N1, [B|N2]) :- pos(N2), split(N, R, N1, N2). %------------------------------------------------------------------------ % Exponentiation and discrete logarithm % n = b^q + r < b^(q+1) % % From the above condition we obtain the upper bound on r: % n >= b^q, n < b^(q+1) = b^q * b = (n-r)* b % r*b < n*(b-1) % % We can also obtain the bounds on q: % if |b| is the bitwidth of b and |n| is the bitwidth of n, % we have, by the definition of the bitwidth: % (1) 2^(|b|-1) <= b < 2^|b| % (2) 2^(|n|-1) <= n < 2^|n| % % Raising (1) to the power of q: % 2^((|b|-1)*q) <= b^q % % OTH, b^q <= n, and n < 2^|n|. So we obtain % (3) (|b|-1)*q < |n| % which defines the upper bound on |q|. % % OTH, raising (1) to the power of (q+1): % b^(q+1) < 2^(|b|*(q+1)) % But n < b^(q+1) by definition of exponentiation, and keeping in mind (1) % (4) |n|-1 < |b|*(q+1) % which is the lower bound on q. % When b = 2, exponentiation and discrete logarithm are easier to obtain % n = 2^q + r, 0<= 2*r < n % Here, we just relate n and q. % exp2 n b q % holds if: n = (|b|+1)^q + r, q is the largest such number, and % (|b|+1) is a power of two. % % Side condition: (|b|+1) is a power of two and b is L-instantiated. % To obtain the binary exp/log relation, invoke the relation as % (exp2 n '() q) % Properties: if n is L-instantiated, one solution, q is fully instantiated. % If q is fully instantiated: one solution, n is L-instantiated. % In any event, q is always fully instantiated in any solution % and n is L-instantiated. % We depend on the properties of split. exp2([l],_,[]). % 1 = b^0 exp2(N,B,[l]) :- gtl(N), split(N,B,[l],_). % n = (|b|+1)^1 + r exp2(N,B,[o|Q1]) :- % n = (2^k)^(2*q) + r pos(Q1), % = (2^(2*k))^q + r lessl(B,N), append(B,[l|B],B2), exp2(N,B2,Q1). exp2(N,B,[l|Q1]) :- % n = (2^k)^(2*q+1) + r pos(Q1), % n/(2^k) = (2^(2*k))^q + r' pos(N1), split(N,B,N1,_), append(B,[l|B],B2), exp2(N1,B2,Q1). % nq = n^q where n is L-instantiated and q is fully instantiated repeated_mul([_|_],[],[l]). % n^0 = 1 where n > 0 repeated_mul(N,[l],N). % n^1 = 1 repeated_mul(N,Q,NQ) :- gtl(Q), add(Q1,[l],Q), repeated_mul(N,Q1,NQ1), mul(NQ1,N,NQ). % We call this predicate log rather than exp due to its close similarity % to div. As the tests at the end show, log can be used for determining % the exact discrete logarithm, logarithm with a residual, exponentiation, % and the n-th root. log([l],B,[],[]) :- pos(B). % 1 = b^0 + 0, b >0 log(N,B,[],R) :- less(N,B), pos(R), add(R,[l],N). % n = b^0 + (n-1) log(N,B,[l],R) :- % n = b + r, n and b the same sz gtl(B), samel(N,B), add(R,B,N). log(N,[l],Q,R) :- pos(Q), add(R,[l],N). % n = 1^q + (n-1), q>0 log(N,[],Q,N) :- pos(Q). % n = 0^q + n, q>0 % in the rest, n is longer than b log(N,[o,l],Q,R) :- % b = 2 pos(N1), N = [_,_|N1], % n is at least 4 exp2(N,[],Q), % that will L-instantiate n and n1 split(N,N1,_,R). log(N,B,Q,R) :- % the general case (B = [l,l]; B = [_,_,_|_]), % B >= 3 lessl(B,N), % b becomes L-instantiated % If b was L-instantiated, the previous % goal had only *one* solution exp2(B,[],Bw1), add(Bw1,[l],Bw), lessl(Q,N), % A _very_ lose bound, but makes % sure q will be L-instantiated % Now, we can use b and q to bound n % |n|-1 < |b|*(q+1) add(Q,[l],Q1), mul(Bw,Q1,Bwq1), % |b|*(q+1) less(Nw1, Bwq1), exp2(N,[],Nw1), % n becomes L-instantiated % Now we have only finite number of ans add(Nw1,[l],Nw), div(Nw, Bw, Ql1, _), % low boundary on q: add(Ql,[l],Ql1), % |n| = |b|(ql+1) + c (samel(Q,Ql); lessl(Ql,Q)), % Tighten the estimate for q repeated_mul(B,Ql,Bql), % bql = b^ql div(Nw,Bw1,Qh,_), % upper boundary on q-1 add(Ql,Qdh,Qh), add(Ql,Qd,Q), (Qd = Qdh; less(Qd,Qdh)), % qd is bounded repeated_mul(B,Qd,Bqd), % b^qd mul(Bql,Bqd,Bq), % b^q mul(B, Bq, Bq1), % b^(q+1) add(Bq,R,N), less(N, Bq1). % check the r condition %---------------------------------------------------------------------- % Tests % A simple unit testing framework shln(N) :- N isb R, print(R), nl. sh1(L,LR) :- maplist((isb),LR,L). sh2(L,LR) :- maplist(sh1,L,LR). okn(B,N) :- N1 isb B, N1 =:= N *-> (print('as expected: '),print(N),nl); throw(mismatch(B,N)). okt(T1,T2) :- T1 = T2 *-> (print('as expected: '),print(T1),nl); throw(mismatch(T1,T2)). should_fail(G) :- G, !, throw(should_have_failed(G)). should_fail(_). ?- nl, print('Addition'), nl. ?- 29 isb A29, 3 isb A3, add(A29,A3,X), okn(X,32). ?- 29 isb A29, 3 isb A3, add(A3,X,A29), okn(X,26). ?- 29 isb A29, 3 isb A3, add(X,A3,A29), okn(X,26). ?- 29 isb A29, 3 isb A3, should_fail(add(_,A29,A3)). ?- print('all numbers that sum to 6'), nl. ?- 6 isb A6, findall([X,Y],add(X,Y,A6),R), sh2(R,LR), okt(LR,[[6, 0], [0, 6], [1, 5], [5, 1], [2, 4], [4, 2], [3, 3]]). ?- print('a few numbers such as X + 1 = Y'), nl. ?- call_with_depth_limit((add(X,[l],Y),print([X,Y]), nl, fail),10,_). % Especially note the third line: 2*x and 2*x +1 are successors, for all x>0! % [[], [l]] % [[l], [o, l]] % [[o, _G523|_G524], [l, _G523|_G524]] % [[l, l], [o, o, l]] % [[l, o, _G538|_G539], [o, l, _G538|_G539]] % [[l, l, l], [o, o, o, l]] % [[l, l, o, _G547|_G548], [o, o, l, _G547|_G548]] ?- nl, print('Subtraction'), nl. ?- 29 isb A29, 3 isb A3, sub(A29,A3,X), okn(X,26). ?- 29 isb A29, 3 isb A3, sub(A29,X,A3), okn(X,26). ?- 26 isb A26, 3 isb A3, sub(X,A3,A26), okn(X,29). ?- 29 isb A29, sub(A29,A29,X), okn(X,0). ?- 29 isb A29, 30 isb A30, should_fail(sub(A29,A30,_)). ?- print('a few numbers such as Y - Z = 4'), nl. ?- call_with_depth_limit((sub(Y,Z,[o,o,l]),print([Y,Z]), nl, fail),10,_). ?- print('a few numbers such as X - Y = Z'), nl. %?- call_with_depth_limit((sub(X,Y,Z),print([X,Y,Z]), nl, fail),10,_). ?- call_with_depth_limit(findall([X,Y,Z],sub(X,Y,Z),R),8,_), print(R), nl. ?- nl, print('Comparisons'), nl. ?- 4 isb A4, less(X,A4), okn(X,0). ?- 40 isb A40, 42 isb A42, less(A40,A42). ?- 40 isb A40, 42 isb A42, should_fail(less(A42,A40)). ?- 6 isb A6, findall(Z,(less(Z,A6)),R), sh1(R,LR), okt(LR,[0, 1, 5, 2, 4, 3]). ?- print('a few numbers that are greater than 4'), nl. ?- call_with_depth_limit((less([o,o,l],Y),print(Y), nl, fail),7,_). % [l, o, l] % [o, o, o, l] % [o, l, o, l] % [l, o, o, l] % [l, l, o, l] ?- nl, print('Multiplication'), nl. ?- 2 isb A2, 3 isb A3, mul(A2,A3,X), okn(X,6). ?- 12 isb A12, 3 isb A3, mul(A3,X,A12), okn(X,4). ?- 12 isb A12, 3 isb A3, mul(X,A3,A12), okn(X,4). ?- 12 isb A12, 5 isb A5, should_fail(mul(_,A5,A12)). ?- print('All factors of 24'), nl. ?- 24 isb A24, findall([X,Y],mul(X,Y,A24),R), sh2(R,LR), okt(LR,[[24, 1], [1, 24], [4, 6], [6, 4], [8, 3], [12, 2], [3, 8], [2, 12]]). ?- print('2 * 3 must produce only one solution'), nl. ?- 2 isb A2, 3 isb A3, findall(Z,mul(A2,A3,Z),R), sh1(R,LR), okt(LR,[6]), nl. ?- nl, print('a few numbers such as 3*X = Y'), nl. ?- call_with_depth_limit((mul([l,l],X,Y),print([X,Y]), nl, fail),9,_). ?- nl, print('a few numbers such as X*3 = Y'), nl. ?- call_with_depth_limit((mul(X,[l,l],Y),print([X,Y]), nl, fail),7,_). ?- nl, print('a few numbers such as X*Y = Z'), nl. ?- call_with_depth_limit((mul(X,Y,Z),print([X,Y,Z]), nl, fail),5,_). ?- nl, print('Recursive enumerability of multiplication'), nl. ?- mul(X,[l,l],Y),X=[l,l], okn(Y,9). ?- mul(X,[l,l],Y),Y=[l,o,o,l], okn(X,3). ?- X=[l,l], mul(X,[l,l],Y), okn(Y,9). ?- mul(N,M,P), (N=[l, o, l], M=[o, l, o, o, l, o, l]; M=[l, o, l], N=[o, l, o, o, l, o, l]), okn(P,410). ?- print('Test to see that mul(X,Y,Z) contains all numbers '), print('such as X*Y=Z, Y<=X'), nl. et(I,J,P) :- mul(X,Y,Z), (X = I, Y = J; X = J, Y = I), Z = P, sh1([X,Y,Z],R), print(R), nl, !. et(L) :- less(I,L), less(J,L), et(I,J,_), fail. et(_). ?- 8 isb L, et(L). ?- nl, print('All Factorizations of 16565'), nl. %?- 16565 isb A,mul(Y,Z,A), sh1([Y,Z],L), print(L), nl, fail. ?- nl, print('Finding all non-primitive factors of 16567'), nl. %?- time((16567 isb A,mul(Y,Z,A), gtl(Y), gtl(Z), sh1([Y,Z],L))). % 63,953,351 inferences in 32.09 seconds (1992704 Lips) % The recursively enumerating mul is slower: %?- time((16567 isb A,mul(Y,Z,A), gtl(Y), gtl(Z), sh1([Y,Z],L))). % SWI-Prolog, Version 5.0: % 453,383,545 inferences in 183.55 seconds (2470124 Lips) % SWI-Prolog, Version 5.4.5: % 453,383,551 inferences, 163.80 CPU in 164.60 seconds (100% CPU, 2767830 Lips) ?- nl, print('Split'), nl. ?- 4 isb A4, findall(Z,(Z=((N1,N2)),split(A4,[],N1,N2)),R), okt(R, [([o,l],[])]). ?- 4 isb A4, findall(Z,(Z=((N1,N2)),split(A4,[l],N1,N2)),R), okt(R, [([l],[])]). ?- 4 isb A4, findall(Z,(Z=((N1,N2)),split(A4,[l,l],N1,N2)),R), okt(R, [([],[o,o,l])]). ?- 4 isb A4, findall(Z,(Z=((N1,N2)),split(A4,[l,l,l],N1,N2)),R), okt(R, [([],[o,o,l])]). ?- 5 isb A5, findall(Z,(Z=((N1,N2)),split(A5,[l],N1,N2)),R), okt(R, [([l],[l])]). ?- 5 isb A5, findall(N,split(N,A5,[],[l]),R), okt(R, [ [l] ]). ?- nl, print('Division'), nl. ?- 2 isb A2, 4 isb A4, div(A4,A2,X,R), sh1([X,R],L), okt(L,[2,0]). ?- print('division by zero must fail.'), nl. ?- 4 isb A4, should_fail(div(A4,[],_,_)). ?- 4 isb A4, div(A4,A4,X,R), sh1([X,R],L), okt(L,[1,0]). ?- 3 isb A3, 4 isb A4, div(A4,A3,X,R), sh1([X,R],L), okt(L,[1,1]). ?- 5 isb A5, 4 isb A4, div(A4,A5,X,R), sh1([X,R],L), okt(L,[0,4]). ?- 2 isb A2, 5 isb A5, div(A5,A2,X,R), sh1([X,R],L), okt(L,[2,1]). ?- 3 isb A3, 33 isb A33, div(A33,A3,X,R), sh1([X,R],L), okt(L,[11,0]). ?- 11 isb A11, 33 isb A33, div(A33,X,A11,R), sh1([X,R],L), okt(L,[3,0]). ?- 11 isb A11, 3 isb A3, div(X,A3,A11,R), sh1([X,R],L), okt(L,[33,0]). ?- 5 isb A5, 33 isb A33, div(A33,A5,X,R), sh1([X,R],L), okt(L,[6,3]). ?- 5 isb A5, 7 isb A7, should_fail(div(A5,_,A7,_)). % Cannot divide anything by 5 and get 5 in the remainder. ?- 5 isb A5, should_fail(div(_,A5,_,A5)). ?- nl, print('all numbers such as 5/Z = 1'), nl. ?- 1 isb A1, 5 isb A5, findall(Z,div(A5,Z,A1,_),R), sh1(R,LR), okt(LR,[5, 4, 3]). ?- nl, print('all inexact factorizations of 24'), nl. ?- 24 isb A24, findall([M,Q,Res],(lessl(M,A24),div(A24,M,Q,Res)),R), sh2(R,LR), okt(LR,[[1, 24, 0], [3, 8, 0], [2, 12, 0], [6, 4, 0], [4, 6, 0], [5, 4, 4], [7, 3, 3], [12, 2, 0], [8, 3, 0], [14, 1, 10], [10, 2, 4], [13, 1, 11], [15, 1, 9], [11, 2, 2], [9, 2, 6]]). ?- nl, print('Even divide'), nl. ?- call_with_depth_limit((div([l|Y],[o,l],Q,R),print([Y,Q,R]), nl, fail),7,_). ?- nl, print('a few numbers such as N = M*Q+R'), nl. ?- call_with_depth_limit((div(N,M,Q,R),print([N,M,Q,R]), nl, fail),6,_). ?- nl, print('N=[], div(N,M,Q,R) gives only one solution!'), nl. ?- findall([M,Q,Res],(N=[], div(N,M,Q,Res)),R), okt(R,[[[_G386|_G387], [], []]]). % The following diverges for simpler division algorithms, such as % div1/4 and div3/4 ?- nl, print('N = 5*7+0 should definitely have only one solution.'), nl. ?- 7 isb Q, 5 isb M, R = [], findall(N,div(N,M,Q,R),NS), sh1(NS,RR), okt(RR,[35]). ?- nl, print('even number divided by 2 cannot give 1 in the remainder'), nl. ?- should_fail(div([o|_],[o,l],_,[l])). ?- nl, print('odd number divided by 2 cannot give 0 in the remainder'), nl. ?- should_fail(div([l,o|_],[o,l],_,[])). ?- nl, print('Binary exponentiation and logarithm'), nl. ?- findall(Q,exp2([l,l,l,l],[],Q),QS), okt(QS,[[l,l]]). ?- findall(Q,exp2([l,o,l,l,l],[],Q),QS), okt(QS,[[o,o,l]]). ?- nl, print('_all_ numbers n such that n = 2^5 + r'), nl. ?- 5 isb A5, findall(N,exp2(N,[],A5),NS), okt(NS,[[o, o, o, o, o, l], [o, l, o, o, o, l], [o, _G924, l, o, o, l], [o, _G945, _G948, l, o, l], [o, _G966, _G969, _G972, l, l], [l, o, o, o, o, l], [l, l, o, o, o, l], [l, _G1029, l, o, o, l], [l, _G1050, _G1053, l, o, l], [l, _G1071, _G1074, _G1077, l, l]]). ?- nl, print('a few numbers such as N = 2^Q+R'), nl. ?- call_with_depth_limit((exp2(N,[],R),sh1([N,R],P), print(P),nl, fail),9,_). ?- nl, print('Exponentiation, logarithm, and base determination'), nl. ?- print('Solving 15 = B^Q + R for different B'), nl. ?- 15 isb A15, 2 isb B, findall([Q,R],log(A15,B,Q,R),RS), sh2(RS,P),okt(P,[[3,7]]). ?- 15 isb A15, 3 isb B, findall([Q,R],log(A15,B,Q,R),RS), sh2(RS,P),okt(P,[[2,6]]). ?- 15 isb A15, 4 isb B, findall([Q,R],log(A15,B,Q,R),RS), sh2(RS,P),okt(P,[[1,11]]). ?- 15 isb A15, 5 isb B, findall([Q,R],log(A15,B,Q,R),RS), sh2(RS,P),okt(P,[[1,10]]). ?- 15 isb A15, 15 isb B, findall([Q,R],log(A15,B,Q,R),RS), sh2(RS,P),okt(P,[[1,0]]). ?- 15 isb A15, 16 isb B, findall([Q,R],log(A15,B,Q,R),RS), sh2(RS,P),okt(P,[[0,14]]). ?- print('Solving 32 = B^Q + R for different B'), nl. ?- 32 isb A32, 2 isb B, findall([Q,R],log(A32,B,Q,R),RS), sh2(RS,P),okt(P,[[5,0]]). ?- 32 isb A32, 3 isb B, findall([Q,R],log(A32,B,Q,R),RS), sh2(RS,P),okt(P,[[3,5]]). ?- 32 isb A32, 4 isb B, findall([Q,R],log(A32,B,Q,R),RS), sh2(RS,P),okt(P,[[2,16]]). ?- 32 isb A32, 5 isb B, findall([Q,R],log(A32,B,Q,R),RS), sh2(RS,P),okt(P,[[2,7]]). ?- 32 isb A32, 6 isb B, findall([Q,R],log(A32,B,Q,R),RS), sh2(RS,P),okt(P,[[1,26]]). ?- print('Finding all B such that 15 = B^3 + R'), nl. ?- 15 isb A15, 3 isb Q, findall([B,R],log(A15,B,Q,R),RS), sh2(RS,P), okt(P,[[1, 14], [0, 15], [2, 7]]). ?- print('Finding all B such that 32 = B^4 + R'), nl. ?- 32 isb N, 4 isb Q, findall([B,R],log(N,B,Q,R),RS), sh2(RS,P), okt(P,[[1, 31], [0, 32]]). ?- print('Finding all B such that 33 = B^5 + R'), nl. ?- 33 isb N, 5 isb Q, findall([B,R],log(N,B,Q,R),RS), sh2(RS,P), okt(P,[[1, 32], [0, 33], [2, 1]]). ?- nl, print('Given specific Q, B, R find all N = B^Q + R'), nl. ?- 1 isb B, 2 isb Q, 0 isb R, findall(N,log(N,B,Q,R),RS), sh1(RS,P),okt(P,[1]). ?- 2 isb B, 5 isb Q, 1 isb R, findall(N,log(N,B,Q,R),RS), sh1(RS,P),okt(P,[33]). ?- 3 isb B, 2 isb Q, 1 isb R, findall(N,log(N,B,Q,R),RS), sh1(RS,P),okt(P,[10]). ?- 3 isb B, 3 isb Q, 1 isb R, findall(N,log(N,B,Q,R),RS), sh1(RS,P),okt(P,[28]). ?- 3 isb B, 3 isb Q, 5 isb R, findall(N,log(N,B,Q,R),RS), sh1(RS,P),okt(P,[32]). ?- nl, print('a few numbers such as N = 3^Q+R'), nl. ?- 3 isb B, call_with_depth_limit((log(N,B,Q,R),sh1([N,B,Q,R],P), print(P),nl, fail),11,_). ?- nl, print('a few numbers such as N = B^3+R'), nl. ?- 3 isb Q, call_with_depth_limit((log(N,B,Q,R),sh1([N,B,Q,R],P), print(P),nl, fail),9,_). ?- nl, print('All done'), nl.