%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% %% %% 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 May 25, 2006. %% %% %% %% This works with SWI-Prolog, version 5.6.12 %% %% by Jan Wielemaker (jan@swi-prolog.org) %% %% (Copyright (c) 1990-2006 University of Amsterdam.) %% %% See: http://www.swi-prolog.org/ %% %% %% %% Load as: %% %% consult('c:/karttu/prolog/altdivlats.txt'). %% %% %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% :- set_prolog_flag(toplevel_print_options,[quoted(true), portray(true), attributes(portray), max_depth(1000)]). msb(X,M) :- M is 1 << floor(log(X)/log(2)). %% In this case zero is also a power of 2. is_power_of_2(0) :- !. is_power_of_2(X) :- msb(X,X). divides(X,Y) :- 0 is Y mod 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). %% %% even_and_odd_divisors(12,EVENS,ODDS). %% EVENS = [4, 6] %% ODDS = [2, 3] %% %% even_and_odd_divisors(24,EVENS,ODDS). %% EVENS = [4, 6] %% ODDS = [2, 3, 8, 12] %% %% even_and_odd_divisors(36,EVENS,ODDS). %% EVENS = [4, 6, 9] %% ODDS = [2, 3, 12, 18] %% even_and_odd_divisors([],EDS,ODS,EDS,ODS) :- !. even_and_odd_divisors([D|Ds],EDS,ODS,ZEDS,ZODS) :- a001222(D,K), 0 is K /\ 1, %% That is, 0 is K rem 2. even_and_odd_divisors(Ds,[D|EDS],ODS,ZEDS,ZODS). even_and_odd_divisors([D|Ds],EDS,ODS,ZEDS,ZODS) :- a001222(D,K), 1 is K /\ 1, %% That is, 1 is K rem 2. even_and_odd_divisors(Ds,EDS,[D|ODS],ZEDS,ZODS). even_and_odd_divisors(N,ZEDS,ZODS) :- propdivisors(N,DIVS), reverse(DIVS,RDIVS), even_and_odd_divisors(RDIVS,[],[],ZEDS,ZODS). %% 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_old(60,M). %% M = [1,2,5,8,19,41,87,138,301,699] divisorbitmasks_old(N,MASKS) :- findall(M,(propdivisors(N,Ds),member(D,Ds),nonunitdivisors(D,DDs),buildmask(DDs,Ds,M)),MASKS). %% We don't hit the 2^26 tagged integer limit as soon as before if we %% build the masks from the prime factorization of each divisor, instead of its divisors: %% divisorbitmasks(60,M). %% M = [1, 4, 3, 8, 5, 9, 7, 12, 11, 13] divisorbitmasks(N,MASKS) :- findall(M,(factorize(N,Fs),propdivisors(N,Ds),member(D,Ds),factorize(D,DFs),buildmask(DFs,Fs,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,ZVISITED) :- !, reverse(VISITED,ZVISITED). linext(UNVIS,VISMASK,VISITED,ZVISITED) :- member(U,UNVIS), %% This drives the backtracking UR is ((VISMASK /\ U) xor U), %% UR is set to the mask of missing prime factor(s). is_power_of_2(UR), %% We continue only if there is max. one prime factor hop from here to there. delete(UNVIS,U,NUNVIS), NVISMASK is VISMASK \/ U, linext(NUNVIS,NVISMASK,[U|VISITED],ZVISITED). linext(N,VISITED) :- divisorbitmasks(N,UNVISITED), linext(UNVISITED,0,[],VISITED). chklinext_old([],_) :- !. chklinext_old(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_old(NUNVIS,NVISMASK). chklinext_old(N) :- divisorbitmasks_old(N,UNVISITED), chklinext_old(UNVISITED,0). chklinext([],_) :- !. chklinext(UNVIS,VISMASK) :- member(U,UNVIS), UR is ((VISMASK /\ U) xor U), %% UR is set to the mask of missing prime factor(s). is_power_of_2(UR), %% We continue only if there is max. one prime factor hop from here to there. delete(UNVIS,U,NUNVIS), NVISMASK is VISMASK \/ U, chklinext(NUNVIS,NVISMASK). chklinext(N) :- divisorbitmasks(N,UNVISITED), chklinext(UNVISITED,0). %% %% "Alternating Linear Extensions" %% altlinext(12,VISITED). %% VISITED = [1, 4, 2, 8] altlinext([],[],_,VISITED,ZVISITED) :- reverse(VISITED,ZVISITED), !. altlinext(UNVIS1,UNVIS2,VISMASK,VISITED,ZVISITED) :- member(U,UNVIS1), msb(U,M), UR is U - M, %% UR is the mask for the prerequisite divisors. UR is UR /\ VISMASK, delete(UNVIS1,U,NUNVIS1), NVISMASK is VISMASK \/ M, altlinext(UNVIS2,NUNVIS1,NVISMASK,[M|VISITED],ZVISITED). %% if N is 36, then EVENMASKS = [5, 11, 18] and ODDMASKS = [1, 2, 47, 91]. altlinext(N,VISITED) :- propdivisors(N,ALLDIVS), even_and_odd_divisors(N,EDIVS,ODIVS), findall(EVENM,(member(E,EDIVS),nonunitdivisors(E,EVENDDS),buildmask(EVENDDS,ALLDIVS,EVENM)),EVENMASKS), findall(ODDM,(member(O,ODIVS),nonunitdivisors(O,ODDDDS),buildmask(ODDDDS,ALLDIVS,ODDM)),ODDMASKS), altlinext(ODDMASKS,EVENMASKS,0,[],VISITED). chkaltlinext([],[],_) :- !. chkaltlinext(UNVIS1,UNVIS2,VISMASK) :- member(U,UNVIS1), msb(U,M), UR is U - M, %% UR is the mask for the prerequisite divisors. UR is UR /\ VISMASK, delete(UNVIS1,U,NUNVIS1), NVISMASK is VISMASK \/ M, chkaltlinext(UNVIS2,NUNVIS1,NVISMASK). chkaltlinext(N) :- propdivisors(N,ALLDIVS), even_and_odd_divisors(N,EDIVS,ODIVS), findall(EVENM,(member(E,EDIVS),nonunitdivisors(E,EVENDDS),buildmask(EVENDDS,ALLDIVS,EVENM)),EVENMASKS), findall(ODDM,(member(O,ODIVS),nonunitdivisors(O,ODDDDS),buildmask(ODDDDS,ALLDIVS,ODDM)),ODDMASKS), chkaltlinext(ODDMASKS,EVENMASKS,0). %% This is already in SWI-Prolog: %% 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)). :- module_transparent succeeds_n_times/2. succeeds_n_times(Goal, Times) :- ( flag(succeeds_n_times, Old, 0), Goal, flag(succeeds_n_times, N, N+1), fail ; flag(succeeds_n_times, Times, 0) %% Was: Old, doesn't work that way. ). %% a114717(N,V) :- %% findall(SOL,linext(N,SOL),SOLS), %% length(SOLS,V). a114717(N,V) :- succeeds_n_times(chklinext(N),V). seq_a114717(UPTON,VALS) :- findall(V,(between(1,UPTON,N),a114717(N,V)),VALS). a001222(N,V) :- factorize(N,Fs), length(Fs,V). seq_a001222(UPTON,VALS) :- findall(V,(between(1,UPTON,N),a001222(N,V)),VALS). a119842(N,V) :- succeeds_n_times(chkaltlinext(N),V). seq_a119842(UPTON,VALS) :- findall(V,(between(1,UPTON,N),a119842(N,V)),VALS). precomp_seq(a025487,[1,2,4,6,8,12,16,24,30,32,36,48,60,64,72,96,120, 128,144,180,192,210,216,240,256,288,360,384,420, 432,480,512,576,720,768,840,864,900,960,1024,1080, 1152,1260,1296,1440,1536,1680,1728,1800,1920,2048, 2160,2304,2310]). %% 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] %% %% seq_a114717(288,S), findmaxpoints(S,M), distpoints(S,DP), distvals(S,DV). %% %% 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, 4877756, 1, 48, 2, 14, 2, 48, 2, 5, 14, 48, 1, 429, 1, 2, 48, 42, 1, 2452, 1, 462, 2, 2, 2, 2452, 2, 2, 5, 42, 2, 1680384, 1, 5, 2, 2, 2, 24024, 2, 2, 2, 2452, 2, 48, 1, 132, 42, 2, 1, 2452, 1, 48, 48, 14, 1, 2452, 2, 5, 2, 48, 1, 17454844, 1, 5, 1, 5, 5, 48, 2, 14, 2, 14, 1, 4877756, 2, 2, 48, 1, 1, 48, 2, 2452, 5, 2, 1, 183958, 2, 48, 2, 5, 1, 183958, 1, 42, 48, 2, 5, 2452, 1, 2, 5, 183958, 1, 48, 1, 5, 48, 48, 2, 87516] %% MP = [1, 6, 12, 24, 30, 60, 120, 180, 240] %% MV = [1, 2, 5, 14, 48, 2452, 183958, 4877756, 17454844] %% DP = [1, 6, 12, 24, 30, 36, 60, 72, 96, 120, 144, 180, 192, 210, 216, 240, 288] (subset of A025487...) %% DV = [1, 2, 5, 14, 48, 42, 2452, 462, 132, 183958, 6006, 4877756, 429, 1680384, 24024, 17454844, 87516] %% %% A119838 --- A119849. %% 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, ... %% seq_a119842(1024,A119842), findmaxpoints(A119842,A119843), findmaxvals(A119842,A119844), distpoints(A119842,A119845), distvals(A119842,A119846), findeqpoints(A119842,0,A119847), findgtpoints(A119842,0,A119848), findgtvals(A119842,0,A119849), findgtpoints(A119842,1,A119850), findgtvals(A119842,1,A119851). 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). %% Find Maximum Values. findmaxvals([],_,MVS,ZMVS) :- !, reverse(MVS,ZMVS). findmaxvals([N|Ns],MAX,MVS,ZMVS) :- N > MAX, !, findmaxvals(Ns,N,[N|MVS],ZMVS). findmaxvals([N|Ns],MAX,MVS,ZMVS) :- N =< MAX, !, findmaxvals(Ns,MAX,MVS,ZMVS). findmaxvals(L,MVS) :- findmaxvals(L,-1,[],MVS). %% find points where the value is greater than GT. findgtpoints([],_,_,MPS,ZMPS) :- !, reverse(MPS,ZMPS). findgtpoints([N|Ns],I,GT,MPS,ZMPS) :- N > GT, !, J is I+1, findgtpoints(Ns,J,GT,[I|MPS],ZMPS). findgtpoints([N|Ns],I,GT,MPS,ZMPS) :- N =< GT, !, J is I+1, findgtpoints(Ns,J,GT,MPS,ZMPS). findgtpoints(L,GT,MPS) :- findgtpoints(L,1,GT,[],MPS). %% %% find values which are greater than GT. findgtvals([],_,_,MVS,ZMVS) :- !, reverse(MVS,ZMVS). findgtvals([N|Ns],I,GT,MVS,ZMVS) :- N > GT, !, J is I+1, findgtvals(Ns,J,GT,[N|MVS],ZMVS). findgtvals([N|Ns],I,GT,MVS,ZMVS) :- N =< GT, !, J is I+1, findgtvals(Ns,J,GT,MVS,ZMVS). findgtvals(L,GT,MVS) :- findgtvals(L,1,GT,[],MVS). %% find points where the value is equal to EQ. findeqpoints([],_,_,MPS,ZMPS) :- !, reverse(MPS,ZMPS). findeqpoints([EQ|Ns],I,EQ,MPS,ZMPS) :- !, J is I+1, findeqpoints(Ns,J,EQ,[I|MPS],ZMPS). findeqpoints([_|Ns],I,EQ,MPS,ZMPS) :- J is I+1, findeqpoints(Ns,J,EQ,MPS,ZMPS). findeqpoints(L,EQ,MPS) :- findeqpoints(L,1,EQ,[],MPS). %% distpoints -- Find Points with Distinct values. distpoints([],_,_,DPS,ZDPS) :- !, reverse(DPS,ZDPS). distpoints([N|Ns],I,DVS,DPS,ZDPS) :- not(memberchk(N,DVS)), !, J is I+1, distpoints(Ns,J,[N|DVS],[I|DPS],ZDPS). distpoints([_|Ns],I,DVS,DPS,ZDPS) :- J is I+1, distpoints(Ns,J,DVS,DPS,ZDPS). distpoints(L,DPS) :- distpoints(L,1,[],[],DPS). %% distvals -- Like above, but returns the distinct values themselves. distvals([],_,DVS,_,ZDVS) :- !, reverse(DVS,ZDVS). distvals([N|Ns],I,DVS,DPS,ZDVS) :- not(memberchk(N,DVS)), !, J is I+1, distvals(Ns,J,[N|DVS],[I|DPS],ZDVS). distvals([_|Ns],I,DVS,DPS,ZDVS) :- J is I+1, distvals(Ns,J,DVS,DPS,ZDVS). distvals(L,DVS) :- distvals(L,1,[],[],DVS). 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) !