%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% %% %% http://..................../~......./....../altdivlats.txt %% %% %% %% A set of Prolog-definitions that illustrate how the first terms %% %% of A0xxxxx are produced. %% %% %% %% Written by Antti Karttunen, 2006, http://www.iki.fi/kartturi/ %% %% %% %% Last edited April 25, 2006. %% %% %% %% This works with GNU prolog: %% %% http://www.gnu.org/software/gprolog/gprolog.html %% %% %% %% Load as: %% %% consult('/karttu/prolog/altdivlats.txt'). %% %% %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% msb(X,M) :- M is 1 << floor(log(X)/log(2)). divides(X,Y) :- 0 is Y rem X. divisors(_,1,Z,[1|Z]) :- !. divisors(N,I,Z,ZZ) :- divides(I,N), !, J is I-1, divisors(N,J,[I|Z],ZZ). divisors(N,I,Z,ZZ) :- J is I-1, divisors(N,J,Z,ZZ). divisors(N,Z,ZZ) :- divisors(N,N,Z,ZZ). divisors(N,Z) :- divisors(N,N,[],Z). %% Like above, but 1 is not included in the list of divisors: nonunitdivisors(_,1,Z,Z) :- !. nonunitdivisors(N,I,Z,ZZ) :- divides(I,N), !, J is I-1, nonunitdivisors(N,J,[I|Z],ZZ). nonunitdivisors(N,I,Z,ZZ) :- J is I-1, nonunitdivisors(N,J,Z,ZZ). nonunitdivisors(N,Z,ZZ) :- nonunitdivisors(N,N,Z,ZZ). nonunitdivisors(N,Z) :- nonunitdivisors(N,N,[],Z). %% Like above, but includes neither 1 nor N in the list of divisors: propdivisors(_,0,Z,Z) :- !. propdivisors(_,1,Z,Z) :- !. propdivisors(N,I,Z,ZZ) :- divides(I,N), !, J is I-1, propdivisors(N,J,[I|Z],ZZ). propdivisors(N,I,Z,ZZ) :- J is I-1, propdivisors(N,J,Z,ZZ). propdivisors(N,Z,ZZ) :- propdivisors(N,N,Z,ZZ). propdivisors(N,Z) :- NN is N-1, propdivisors(N,NN,[],Z). %% divisors(60,Z). %% produces Z = [1,2,3,4,5,6,10,12,15,20,30,60] %% factorize(360,Z). %% --> Z = [2,2,2,3,3,5] factorize(1,_,Z,ZZ) :- reverse(Z,ZZ), !. factorize(N,I,Z,ZZ) :- divides(I,N), !, NN is N//I, factorize(NN,I,[I|Z],ZZ). factorize(N,I,Z,ZZ) :- J is I+1, factorize(N,J,Z,ZZ). factorize(N,Z) :- factorize(N,2,[],Z). selectbymask(0,_,L2,L3) :- reverse(L2,L3), !. selectbymask(M,[E1|L1],L2,L3) :- 1 is M /\ 1, %% That is, 1 is M rem 2. !, MM is M // 2, selectbymask(MM,L1,[E1|L2],L3). selectbymask(M,[_|L1],L2,L3) :- 0 is M /\ 1, %% That is, 0 is M rem 2. !, MM is M // 2, selectbymask(MM,L1,L2,L3). selectbymask(M,L1,L3) :- selectbymask(M,L1,[],L3). buildindices([],_,_,NINS,RES) :- !, reverse(NINS,RES). buildindices([A|L1],[A|L2],I,OINS,RES) :- !, J is I+1, buildindices(L1,L2,J,[I|OINS],RES). buildindices(L1,[_|L2],I,OINS,RES) :- !, J is I+1, buildindices(L1,L2,J,OINS,RES). buildindices(L1,L2,RES) :- buildindices(L1,L2,0,[],RES). buildmask([],_,_,NMASK,NMASK) :- !. buildmask([A|L1],[A|L2],I,OMASK,RES) :- !, J is I << 1, NMASK is OMASK + I, buildmask(L1,L2,J,NMASK,RES). buildmask(L1,[_|L2],I,OMASK,RES) :- !, J is I << 1, buildmask(L1,L2,J,OMASK,RES). buildmask(L1,L2,RES) :- buildmask(L1,L2,1,0,RES). %% divisors(12,Ds). %% Ds = [1,2,3,4,6,12] %% propdivisors(60,Ds). %% Ds = [2,3,4,5,6,10,12,15,20,30] %% findall(DDs,(divisors(12,Ds),member(D,Ds),divisors(D,DDs)),ALL). %% ALL = [[1],[1,2],[1,3],[1,2,4],[1,2,3,6],[1,2,3,4,6,12]] %% findall(DDs,(propdivisors(12,Ds),member(D,Ds),propdivisors(D,DDs)),ALL). %% ALL = [[],[],[2],[2,3]] %% findall(DDs,(propdivisors(60,Ds),member(D,Ds),propdivisors(D,DDs)),ALL). %% ALL = [[],[],[2],[],[2,3],[2,5],[2,3,4,6],[3,5],[2,4,5,10],[2,3,5,6,10,15]] %% findall(M,(propdivisors(60,Ds),member(D,Ds),propdivisors(D,DDs),buildmask(DDs,Ds,M)),ALL). %% ALL = [0,0,1,0,3,9,23,10,45,187] %% findall(M,(propdivisors(60,Ds),member(D,Ds),buildmask([D],Ds,M)),ALL). %% ALL = [1,2,4,8,16,32,64,128,256,512] %% divisorbitmasks(60,M). %% M = [1,2,5,8,19,41,87,138,301,699] divisorbitmasks(N,MASKS) :- findall(M,(propdivisors(N,Ds),member(D,Ds),nonunitdivisors(D,DDs),buildmask(DDs,Ds,M)),MASKS). divisorflags(N,BITS) :- findall(M,(propdivisors(N,Ds),member(D,Ds),buildmask([D],Ds,M)),BITS). %% For each divisor, form a bit-mask of its own divisors %% (excluding it itself, and possibly also 1). %% This is then anded with the "traversed divisors" -mask %% and if the result is equal, then the search is %% allowed to continue to this divisor. %% linext([],_,VISITED,VISITED) :- !. linext(UNVIS,VISMASK,VISITED,ZVISITED) :- member(U,UNVIS), msb(U,M), UR is U - M, %% UR is the mask for the prerequisite divisors. UR is UR /\ VISMASK, delete(UNVIS,U,NUNVIS), NVISMASK is VISMASK \/ M, linext(NUNVIS,NVISMASK,[M|VISITED],ZVISITED). linext(N,VISITED) :- divisorbitmasks(N,UNVISITED), linext(UNVISITED,0,[],VISITED). chklinext([],_) :- !. chklinext(UNVIS,VISMASK) :- member(U,UNVIS), msb(U,M), UR is U - M, %% UR is the mask for the prerequisite divisors. UR is UR /\ VISMASK, delete(UNVIS,U,NUNVIS), NVISMASK is VISMASK \/ M, chklinext(NUNVIS,NVISMASK). chklinext(N) :- divisorbitmasks(N,UNVISITED), chklinext(UNVISITED,0). not(P) :- call(P), !, fail. not(_). %% for_all(member(X,[1,2,3,4,5]),X < 100). %% yes %% From Sterling, Shapiro, The Art of Prolog, page 302: for_all(G,C) :- not((G,not(C))). %% From Sterling, Shapiro, The Art of Prolog, page 200: disjoint(Xs,Ys) :- not((member(Z,Xs), member(Z,Ys))). :- dynamic(cnt/1). cnt(_) :- fail. addcount(_) :- cnt(X), !, Y is X+1, retract(cnt(X)), asserta(cnt(Y)). addcount(_) :- asserta(cnt(1)). countsolutions(P,N) :- for_all(P,addcount(1)), cnt(N), retract(cnt(N)). a114717(N,V) :- countsolutions(chklinext(N),V). seq_a114717(UPTON,VALS) :- findall(V,(for(N,1,UPTON),a114717(N,V)),VALS). %% seq_a114717(179,S), findmaxpoints(S,Z). %% %% S = [1,1,1,1,1,2,1,1,1,2,1,5,1,2,2,1,1,5,1,5,2,2,1,14,1,2,1,5,1,48,1,1,2,2,2,42,1,2,2,14,1,48,1,5,5,2,1,42,1,5,2,5,1,14,2,14,2,2,1,2452,1,2,5,1,2,48,1,5,2,48,1,462,1,2,5,5,2,48,1,42,1,2,1,2452,2,2,2,14,1,2452,2,5,2,2,2,132,1,5,5,42,1,48,1,14,48,2,1,462,1,48,2,42,1,48,2,5,5,2,2,183958,1,2,2,5,1,2452,1,1,2,48,1,2452,2,2,14,14,1,48,1,2452,2,2,2,6006,2,2,5,5,1,2452,1,14,5,48,2,2452,1,2,2,132,2,42,1,5,48,2,1,183958,1,48,5,5,1,48,5,42,2,2,1] %% Z = [1,6,12,24,30,60,120] %% %% C.f. A018894 Sigma(n)/phi(n) sets a new record: %% 2, 4, 6, 12, 24, 30, 60, 120, 180, 210, 360, 420, 840, 1260, 1680, 2520, 4620, 9240, ... findmaxpoints([],_,_,MPS,ZMPS) :- !, reverse(MPS,ZMPS). findmaxpoints([N|Ns],I,MAX,MPS,ZMPS) :- N > MAX, !, J is I+1, findmaxpoints(Ns,J,N,[I|MPS],ZMPS). findmaxpoints([N|Ns],I,MAX,MPS,ZMPS) :- N =< MAX, !, J is I+1, findmaxpoints(Ns,J,MAX,MPS,ZMPS). findmaxpoints(L,MPS) :- findmaxpoints(L,1,-1,[],MPS). divdivisors(N,DDs) :- divisors(N,Ds), member(D,Ds), divisors(D,DDs). %% findall(L,(divdivisors(24,DDs),length(DDs,L)),Z). %% --> Z = [8,6,4,4,3,2,2,1] %% To find all the prime-divisors: %% findall(D,(divisors(60,Ds),member(D,Ds),divisors(D,DDs),length(DDs,2)),ALL). %% --> ALL = [5,3,2] noZdividesN(_,[]) :- !. noZdividesN(N,[Z|Zs]) :- divides(Z,N) -> fail; noZdividesN(N,Zs). isprimedivisor(I,N,Z) :- divides(I,N), noZdividesN(I,Z). %% This is stupid, use your head! (Divide the primes out of N as you recurse!) %% (Now done in factorize, above.) %% This ending condition still needed for squareful integers: primedivisors(N,M,_,Z,ZZ) :- M > N, reverse(Z,ZZ), !. %% This cuts the searching time dramatically for square-free integers: primedivisors(N,_,N,Z,ZZ) :- reverse(Z,ZZ), !. primedivisors(N,I,P,Z,ZZ) :- isprimedivisor(I,N,Z), !, J is I+1, NP is I*P, primedivisors(N,J,NP,[I|Z],ZZ). primedivisors(N,I,P,Z,ZZ) :- J is I+1, primedivisors(N,J,P,Z,ZZ). primedivisors(N,Z) :- primedivisors(N,2,1,[],Z). %% primedivisors(900,Z). %% --> Z = [5,3,2] %% Need: mprimedivisors (with multiplicity) (i.e. factorize) !