######################################################################
##TOWERS: Save this file as  TOWERS                                  #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read TOWERS<Enter>                                 #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#zeilberg at math dot rutgers dot edu                                #
######################################################################
 
#Created: Oct.-Dec., 2012
 
print(`Created:  Oct.-Dec., 2012`):
print(` This is TOWERS `):
print(`It is a  companion of the article `):
print(`Automated Counting of Towers (A La Bordelaise)`):
print(`[Or: Footnote to p. 81 of the Flajolet-Sedgewick Chef-d'oeuvre]`):
print(`by Shalosh B. Ekhad and Doron Zeilberger`):
print(``):
print(`Pour Philippe Flajolet, in Memoriam `):

print(`and also available from Zeilberger's website`):
print(``):
print(`Please report bugs to zeilberg at math dot rutgers dot edu`):
print(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://www.math.rutgers.edu/~zeilberg/  .`):

 print(`For a list of the supporting procedures type ezra1();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`--------------------------------------------------------------`):

print(`--------------------------------------------------------------`):
 print(`For a list of the MAIN procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`--------------------------------------------------------------`):


print(`--------------------------------------------------------------`):
 print(`For a list of the procedures for specific problems`):
 print(` type ezraS();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`--------------------------------------------------------------`):

print(`--------------------------------------------------------------`):
 print(`For a list of the Ad hoc procedures`):
 print(` type ezraAd();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`--------------------------------------------------------------`):


print(`--------------------------------------------------------------`):
 print(`For a list of the  bijection procedures`):
 print(` type ezraB();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`--------------------------------------------------------------`):

with(combinat):

ezraB:=proc()

if args=NULL then
 print(` The bijection procedures are: `):
 print(` CanPass, Land, Land1 `):
 print(` SetH12nz `):

else
ezra(args):
fi:

end:

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: CheckBricks, Cont`):
 print(`ContNF, ContSF, Empir `):
 print(` FindBrick, Findrec, IsComp, Meda,`):
 print(` Mishkal, PfromH, PfromHL, PfromHLt, ShiftR,`):
 print(` Size1, TouchingBricks, Wt, Zinn  `):
 print(``):

else
ezra(args):
fi:

end:


ezraAd:=proc()

if args=NULL then
 print(` The Ad hoc procedures are:`):
 print(`A12f, H12Bab, H12nzBab, P12, M12 `):
 print(` `):
 print(` H12 `):
 print(` H12nz, Zset `):
 print(``):

else
ezra(args):
fi:

end:

ezraS:=proc()

if args=NULL then
 print(` The procedures about specific countings are:`):
 print(` ITP `):
 print(` KamaOneBH, KamaOneBP, KamaOneBM, `):
 print(` KamaOneBHfe, KamaOneBPfe, KamaOneBMfe `):
 print(` KamaOneBHnz, KamaOneBPnz, KamaOneBMnz, `):
 print(` KamaOneBHnzFe, KamaOneBPnzFe, KamaOneBMnzFe, `):
 print(` KamaManyBH, KamaManyBP, KamaManyBM, `):
  print(` KamaManyBHnz, KamaManyBPnz, KamaManyBMnz `):
 print(` ManyBH, ManyBP, ManyBM, `):
 print(` ManyBHnz, ManyBPnz, ManyBMnz, `):
 print(` WtManyBH, WtManyBP, WtManyBM, `):
 print(` WtManyBHnz, WtManyBPnz, WtManyBMnz `):
 print(``):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are: KamaH, KamaM, KamaP, `):
 print(`Hk,HL, Pk, PL, Mk, ML `):
 print(`HPyramidot, Migdalim,`):
 print(` Mamar1,  Mamar2, Mamar3 `):
 print(` Mishpat1,  Mishpat2, Mishpat3 `):
 print(` PolH, Pyramidot `):
 print(`SeqH,  SeqM, SeqP `):
 print(` `):



elif nops([args])=1 and op(1,[args])=A12f then
print(`A12f(t): the explicit formulas for the `):
print(`generating functions in t, for the sequences`):
print(`number of (Half-Pyramids, Pyramids,and TOWERS)`):
print(`with bricks {1,2} with`):
print(`all interfaces allowed`):
print(`Try A12f(t);`):


elif nops([args])=1 and op(1,[args])=CanPass then
print(`CanPass(R,p): can the piece p pass through the row of bricks R?`):
print(`For example, try:`):
print(`CanPass([[0,2],[2,4]],[1,3]);`):

elif nops([args])=1 and op(1,[args])=CheckBricks then
print(`CheckBricks(L,M): given a list L, of positive integers`):
print(`and a list of length nops(L) of lists of of length nops(L)`):
print(`checks whether it is a valid input for buliding bricks`):
print(`the condition is that M[i][j] should be`):
print(`a subset of {-({L[j]-1)...(L[i]-1)} indicating the allowed`):
print(`interfaces of placing brick L[j] on TOP of brick L[i]`):
print(`For example, try:`):
print(`CheckBricks([2],[[{-1,1}]]);`):
print(`CheckBricks([1, 2], [[{0}, {-1}], [{0, 1}, {-1, 1}]]);`):
print(``):

elif nops([args])=1 and op(1,[args])=Cont then
print(`Cont(L,M,T,N): given a list of bricks (set of pos. integers) L`):
print(`and a list-of-lists of sets M, describing the construction`):
print(`rules for building towers, and a tower T, and a positive integer`):
print(` N,outputs the set of towers obtained by adding ONE brick`):
print(`to the current top floor of T, with size still <=N `):
print(` For example, try:`):
print(`Cont([2],[[{-1,1}]],[[[0,2],[2,4]],[[-1,1]]],20);`);
print(`Cont([1,2],[ [{0},{-1,0}],[{0,1},{-1,0,1}]], [ [[0,2],[2,3],[3,5]], [[1,3]]],20 );`):

elif nops([args])=1 and op(1,[args])=ContNF then
print(`ContNF(L,M,T,N): Given a list of bricks (set of pos. integers) L`):
print(`and a list-of-lists of sets M, describing the construction`):
print(`rules for building towers, and a tower T and a positive integer N`):
print(`outputs the set of towers obtained by creating a new floor`):
print(`by adding ONE brick`):
print(`For example, try:`):
print(`ContNF([1,2],[ [{0},{-1,0}],[{0,1},{-1,0,1}]], [ [[0,2],[2,3],[3,5]], [[1,3]]],20 );`):

elif nops([args])=1 and op(1,[args])=ContSF then
print(`ContSF(L,M,T,N): given a list of bricks (set of pos. integers) L`):
print(`and a list-of-lists of sets M, describing the construction`):
print(`rules for building towers, and a tower T, and a positive integer`):
print(` N,outputs the set of towers obtained by adding ONE brick`):
print(`to the current top floor of T, with size still <=N `):
print(` For example, try:`):
print(`ContSF([2],[[{-1,1}]],[[[0,2],[2,4]],[[-1,1]]],20);`);
print(`ContSF([1,2],[ [{0},{-1,0}],[{0,1},{-1,0,1}]], [ [[0,2],[2,3],[3,5]], [[1,3]]],20 );`):


elif nops([args])=1 and op(1,[args])=Empir then
print(`Empir(gu,x,P): given a sequence of numbers gu, `):
print(`guesses an algebraic`):
print(`equation satisfied by P=1+add(gu[i]*x^i,i=1..infinity).`):
print(`Try:`):
print(`Empir([seq(binomial(2*i,i)/(i+1),i=1..40)],x,P);`):

elif nops([args])=1 and op(1,[args])=FindBrick then
print(`FindBrick(Koma,x): Given a floor Koma, and an integer x`):
print(`finds which brick of Koma lies under it (if it is above)`):
print(`a brick, or if it is between two bricks it`):
print(`returns an integer+1/2. For example, try:`):
print(`FindBrick([[-2,3],[5,7]], -1);`):


elif nargs=1 and args[1]=Findrec then
print(`Findrec(f,n,N,MaxC): Given a list f tries to find a linear recurrence equation with`):
print(`poly coffs. ope(n,N), where n is the discrete variable, and N is the shift operator `):
print(`of maximum DEGREE+ORDER<=MaxC`):
print(`e.g. try Findrec([1,1,2,3,5,8,13,21,34,55,89],n,N,2);`):

elif nops([args])=1 and op(1,[args])=H12Bab then
print(`H12Bab(z,N): the first N terms of the sequence of`):
print(`polynomials "weight-enumerator of half-pyramids`):
print(`(with all interfaces allowed) with bricks of size`):
print(`1 and 2.  Using Ad-hoc quadratic equation solving `):
print(`known to the Babylonians. `):
print(`Try:`):
print(`H12Bab(z,10);`):

elif nops([args])=1 and op(1,[args])=H12ShortCut then
print(`H12ShortCut(z,N): the first N terms of the sequence of`):
print(`polynomials "weight-enumerator of half-pyramids`):
print(`(with all interfaces allowed) with bricks of size`):
print(`1 and 2.  Using Iterating the Functional Equation`):
print(`Try:`):
print(`H12ShortCut(z,10);`):


elif nops([args])=1 and op(1,[args])=H12 then
print(`H12(z,N): the first N terms of the sequence of`):
print(`polynomials "weight-enumerator of half-pyramids`):
print(`(with all interfaces allowed) with bricks of size`):
print(`1 and 2.  Using Iterating the Functional Equation`):
print(`For the system of size 2, i.e. without using the shortcut`):
print(`Try:`):
print(`H12(z,10);`):

elif nops([args])=1 and op(1,[args])=H12nz then
print(`H12nz(z,N): the first N terms of the sequence of`):
print(`polynomials "weight-enumerator of half-pyramids`):
print(`(with all interfaces allowed EXCEPT 0) with bricks of size`):
print(`1 and 2.  Using Iterating the Functional Equation`):
print(`For the system of size 2, i.e. without using the shortcut`):
print(`Try:`):
print(`H12nz(z,10);`):



elif nops([args])=1 and op(1,[args])=Hk then
print(`Hk(k,z,N): the first N terms of the sequence of`):
print(`polynomials "weight-enumerator of half-pyramids`):
print(`(with all interfaces allowed) with bricks of size`):
print(`1, 2,..k. Try:`):
print(`Hk(3,z,10);`):

elif nops([args])=1 and op(1,[args])=HL then
print(`HL(L,z,N): the first N terms of the sequence of`):
print(`polynomials "weight-enumerator of half-pyramids`):
print(`(with all interfaces allowed) with bricks of sizes`):
print(`in the list L. Try:`):
print(`HL([1,3,4],z,10);`):


elif nops([args])=1 and op(1,[args])=Land1 then
print(`Land1(T,p): landing a piece p from heaven onto`):
print(`a tower T, getting a larger tower`):
print(`try:`):
print(`Land1([[[0,2]]],[0,1]);`):

elif nops([args])=1 and op(1,[args])=Land then
print(`Land(M1,M2): landing a tower M2 from heaven onto`):
print(`a tower M1, getting a larger tower`):
print(`try:`):
print(`Land([[[0,2]]],[[[0,1]]]);`):


elif nops([args])=1 and op(1,[args])=Mamar1 then
print(`Mamar1(L1,t,H,P,M,n,MaxC,K,Godel1,Godel2):`):
print(`outputs a paper about Mishpat1(L,t,H,P,M,n,MaxC,K,Godel1,Godel2):`):
print(`for all subsets of {1, ..., L1}  of size>=2. Try`):
print(`Mamar1(3,t,H,P,M,n,15,1000,30,100):`):

elif nops([args])=1 and op(1,[args])=Mamar2 then
print(`Mamar2(L1,t,H,P,M,n,MaxC,K,Godel1,Godel2):`):
print(`outputs a paper about Mishpat2(L,t,H,P,M,n,MaxC,K,Godel1,Godel2):`):
print(`for L from 2  to L1, try:`):
print(`Mamar2(3,t,H,P,M,n,15,1000,30,100):`):

elif nops([args])=1 and op(1,[args])=Mamar3 then
print(`Mamar3(L1,t,H,P,M,n,MaxC,K,Godel1,Godel2):`):
print(`outputs a paper about Mishpat3(L,t,H,P,M,n,MaxC,K,Godel1,Godel2):`):
print(`for L from 2  to L1, try:`):
print(`Mamar3(3,t,H,P,M,n,15,1000,30,100):`):

elif nops([args])=1 and op(1,[args])=Meda then
print(`Meda(Li,x,P,n,N,MaxC,K): Given a sequence Li, outputs, if possible`):
print(`(i) a (conjectured, but often proved [depending on the context])`):
print(` algebraic equation satisfied by its ordinary generating function`):
print(`denoted by P(x)=1+add(L[i]*x^i,i=1..infinity) (but of course only`):
print(`using nops(Li) terms) (it gives FAIL if it fails)`):
print(``):
print(`(ii) a (conjectured, but often routinely provable`):
print(` [depending on the context]) linear recurrence operator `):
print(` in n and the shift-operator N annihilating the sequence `):
print(`of complexity <=MaxC`):
print(`(iii) Its K-th term (from the operator, if found) otherwise FAIL`):
print(`(iv) empirically derived (but once again provable) asymptotics`):
print(`For example, try:`):
print(`Meda([seq(binomial(2*i,i)/(i+1),i=1..50)],x,P,n,N,10,1000);`):

elif nops([args])=1 and op(1,[args])=M12 then
print(`M12(z,N): the first N terms of the sequence of`):
print(`polynomials "weight-enumerator of half-pyramids`):
print(`(with all interfaces allowed) with bricks of size`):
print(`1 and 2, using Ad-hoc quadratic equation solving.`):
print(` Try:`):
print(`M12(z,10);`):

elif nops([args])=1 and op(1,[args])=Mishpat1 then
print(`Mishpat1(L,t,H,P,M,n,MaxC,K,Godel1,Godel2,co):`):
print(`inputs:`):
print(`(i) a list, L, of increasing positive integers`):
print(`(ii) a symbol t`):
print(`(iii, iv, v) symbols H,P,M (for ordinary generating functions)`):
print(`(vi) symbol n`):
print(`(vii) a positive integer MaxC (for the maximal complexity of the recurrence)`):
print(`(viii) a large positive integer (for the asympotics and for displaying`):
print(`a specific value)`):
print(`(ix,x) positive integers Godel1 and Godel2 (Godel1<Godel2)`):
print(`(xi): a pos. integer co, for labelling the Theorem`):
print(`Outputs: a theorem about the sequences enumerating`):
print(`half-pyramids, pyramids, and towers consisting of k-mers`):
print(`whose sizes,k, are drawn from the list L.`):
print(`For example, try:`):
print(`Mishpat1([1,2],t,H,P,M,n,10,1000,30,100,1):`):



elif nops([args])=1 and op(1,[args])=Mishpat2 then
print(`Mishpat2(L,t,H,P,M,n,MaxC,K,Godel1,Godel2,co):`):
print(`inputs:`):
print(`(i) a pos. integer L, `):
print(`(ii) a symbol t`):
print(`(iii, iv, v) symbols H,P,M (for ordinary generating functions)`):
print(`(vi) symbol n`):
print(`(vii) a positive integer MaxC (for the maximal complexity of the recurrence)`):
print(`(viii) a large positive integer (for the asympotics and for displaying`):
print(`a specific value)`):
print(`(ix,x) positive integers Godel1 and Godel2 (Godel1<Godel2)`):
print(`(xi): a pos. integer co, for labelling the Theorem`):
print(`Outputs: a theorem about the sequences enumerating`):
print(`half-pyramids, pyramids, and towers consisting of L-mers`):
print(`all of the same size, where ALL interfaces are allowed.`):
print(`For example, try:`):
print(`Mishpat2(2,t,H,P,M,n,10,1000,30,100,1):`):

elif nops([args])=1 and op(1,[args])=Mishpat3 then
print(`Mishpat3(L,t,H,P,M,n,MaxC,K,Godel1,Godel2,co):`):
print(`inputs:`):
print(`(i) a pos. integer L, `):
print(`(ii) a symbol t`):
print(`(iii, iv, v) symbols H,P,M (for ordinary generating functions)`):
print(`(vi) symbol n`):
print(`(vii) a positive integer MaxC (for the maximal complexity of the recurrence)`):
print(`(viii) a large positive integer (for the asympotics and for displaying`):
print(`a specific value)`):
print(`(ix,x) positive integers Godel1 and Godel2 (Godel1<Godel2)`):
print(`(xi): a pos. integer co, for labelling the Theorem`):
print(`Outputs: a theorem about the sequences enumerating`):
print(`half-pyramids, pyramids, and towers consisting of L-mers`):
print(`all of the same size, where ALL interfaces are allowed.`):
print(`except EXACTLY lined-up`):
print(`For example, try:`):
print(`Mishpat3(2,t,H,P,M,n,10,1000,30,100,1):`):

elif nops([args])=1 and op(1,[args])=Mk then
print(`Mk(k,z,N): the first N terms of the sequence of`):
print(`polynomials "weight-enumerator of towers`):
print(`(with all interfaces allowed) with bricks of size`):
print(`1, 2,..k. Try:`):
print(`Mk(3,z,5);`):

elif nops([args])=1 and op(1,[args])=ML then
print(`ML(L,z,N): the first N terms of the sequence of`):
print(`polynomials "weight-enumerator of towers`):
print(`(with all interfaces allowed) with bricks of sizes`):
print(`given by the list L. Try:`):
print(`ML([1,2,3],z,5);`):


elif nops([args])=1 and op(1,[args])=P12 then
print(`P12(z,N): the first N terms of the sequence of`):
print(`polynomials "weight-enumerator of pyramids`):
print(`(with all interfaces allowed) with bricks of size`):
print(`1 and 2, using Ad-hoc quadratic equation solving.`):
print(` Try:`):
print(`P12(z,10);`):


elif nops([args])=1 and op(1,[args])=PfromH then
print(`PfromH(k,z,t,H): expressing the generating function`):
print(`of pyramids with bricks of size {1, ....,k}`):
print(`in terms of half-pyramids where H is the symbol for`):
print(`the generating function of half-pyramids`):
print(`Try:`):
print(`PfromH(2,z,t,H);`):


elif nops([args])=1 and op(1,[args])=PfromHL then
print(`PfromHL(L,z,t,H): expressing the generating function`):
print(`of pyramids with bricks of sizes given by the list L`):
print(`in terms of half-pyramids where H is the symbol for`):
print(`the generating function of half-pyramids`):
print(`with weight given by z[i]'s`):
print(`Try:`):
print(`PfromHL([1,2],z,t,H);`):

elif nops([args])=1 and op(1,[args])=PfromHLt then
print(`PfromHLt(L,t,H): expressing the generating function`):
print(`of pyramids with bricks of sizes given by the list L`):
print(`in terms of half-pyramids where H is the symbol for`):
print(`the generating function of half-pyramids`):
print(`with all weight 1 (just enumeratin)`):
print(`Try:`):
print(`PfromHLt([1,2],t,H);`):

elif nops([args])=1 and op(1,[args])=Pk then
print(`Pk(k,z,N): the first N terms of the sequence of`):
print(`polynomials "weight-enumerator of pyramids`):
print(`(with all interfaces allowed) with bricks of size`):
print(`1, 2,..k. Try:`):
print(`Pk(3,z,5);`):

elif nops([args])=1 and op(1,[args])=PL then
print(`PL(L,z,N): the first N terms of the sequence of`):
print(`polynomials "weight-enumerator of pyramids`):
print(`(with all interfaces allowed) with bricks of sizes`):
print(`given in the list L. Try:`):
print(`PL([1,2,3],z,5);`):

elif nops([args])=1 and op(1,[args])=PolH then
print(`PolH(L,H,t): the polynomial, P, in H and t such that`):
print(`P(H(t),t)=0, where H(t) is the ordinary generating function`):
print(`for enumerating towers of k-mers of sizes drawn from L`):
print(`For example, try:`):
print(`PolH([1,3],H,t);`):

elif nops([args])=1 and op(1,[args])=HPyramidot then
print(`HPyramidot(L,M,N): Given a list of bricks (set of pos. integers) L`):
print(`and a list-of-lists of sets M, outputs a list of length N`):
print(`whose i-th item is the set of HALF-PYRAMIDS of size i`):
print(`built from the bricks L using the interface rules M`):
print(`A pyramid is a tower whose first floor only has ONE brick`):
print(`For example, try:`):
print(`HPyramidot([2],[[{-1,1}]],5);`):
print(`HPyramidot([1,2],[ [{0},{-1,0}],[{0,1},{-1,0,1}]],5 );`):

elif nops([args])=1 and op(1,[args])=IsComp then
print(`IsComp(L,M,B1,B2): Is it legal to place brick B2 on top`):
print(`of brick B1, with the list of bricks L and `):
print(`rules given by the list-of-lists M? Try:`):
print(`IsComp([1,2],[ [{0},{-1,0}],[{0,1},{-1,0,1}]],[4,6],[5,7]);`):

elif nops([args])=1 and op(1,[args])=ITP then
print(`ITP(L,F,t,N): Given a list of polynomials L`):
print(`(let's n=nops(L))`):
print(`in the variable t and the indexed variable`):
print(`F[1], ..., F[N], describing a polynomial transformation`):
print(`of formal power series. Iterates it N times. For`):
print(`example try`):
print(`ITP([z[1]*t+z[1]*t*H[1]+z[1]*t*H[2],z[2]*t^2+2*z[2]*t^2*(H[1]+H[2])+z[2]*t^2*(H[1]+H[2])^2],H,t,5);`):

elif nops([args])=1 and op(1,[args])=KamaH then
print(`KamaH(L,M,N): Given a list of bricks (set of pos. integers) L`):
print(`and a list-of-lists of sets M, outputs a list of length N`):
print(`whose i-th item is the NUMBER of Half-Pyramids of size i`):
print(`built from the bricks L using the interface rules M.`):
print(`For example, try:`):
print(`KamaH([2],[[{-1,1}]],10);`);
print(`KamaH([1,2],[ [{0},{-1,0}],[{0,1},{-1,0,1}]],10 );`):


elif nops([args])=1 and op(1,[args])=KamaOneBH then
print(`KamaOneBH(k,N): The first N terms of the sequence:`):
print(`Number of Half-Pyramids with N bricks, where only one kind of brick, `):
print(` of dimensions 1 by k participate and `):
print(`all interfaces are allowed, try:`):
print(`KamaOneBH(3,5);`):

elif nops([args])=1 and op(1,[args])=KamaOneBHfe then
print(`KamaOneBHfe(k,N): The first N terms of the sequence:`):
print(`Number of Half-Pyramids with N bricks, where only one kind of brick, `):
print(` of dimensions 1 by k participate and `):
print(`all interfaces are allowed. Using the Functional Equation. Try:`):
print(`KamaOneBHfe(3,5);`):

elif nops([args])=1 and op(1,[args])=KamaOneBPfe then
print(`KamaOneBPfe(k,N): The first N terms of the sequence:`):
print(`Number of Pyramids with N bricks, where only one kind of brick, `):
print(` of dimensions 1 by k participate and `):
print(`all interfaces are allowed. Using the Functional Equation. Try:`):
print(`KamaOneBPfe(3,5);`):

elif nops([args])=1 and op(1,[args])=KamaOneBMfe then
print(`KamaOneBMfe(k,N): The first N terms of the sequence:`):
print(`Number of Towers with N bricks, where only one kind of brick, `):
print(` of dimensions 1 by k participate and `):
print(`all interfaces are allowed. Using the Functional Equation. Try:`):
print(`KamaOneBMfe(3,5);`):

elif nops([args])=1 and op(1,[args])=KamaOneBM then
print(`KamaOneBM(k,N): The first N terms of the sequence:`):
print(`Number of towers with N bricks, where only one kind of brick, `):
print(` of dimensions 1 by k participate and `):
print(`all interfaces are allowed, try:`):
print(`KamaOneBM(3,5);`):

elif nops([args])=1 and op(1,[args])=KamaOneBP then
print(`KamaOneBP(k,N): The first N terms of the sequence:`):
print(`Number of Pyramids with N bricks, where only one kind of brick, `):
print(` of dimensions 1 by k participate and `):
print(`all interfaces are allowed, try:`):
print(`KamaOneBP(3,5);`):




elif nops([args])=1 and op(1,[args])=KamaOneBHnz then
print(`KamaOneBHnz(k,N): The first N terms of the sequence:`):
print(`Number of Half-Pyramids with N bricks, where only one kind of brick, `):
print(` of dimensions 1 by k participate and `):
print(`all interfaces are allowed, except the 0 one. Try:`):
print(`KamaOneBHnz(3,5);`):

elif nops([args])=1 and op(1,[args])=KamaOneBHnzFe then
print(`KamaOneBHnzFe(k,N): The first N terms of the sequence:`):
print(`Number of Half-Pyramids with N bricks, where only one kind of brick, `):
print(` of dimensions 1 by k participate and `):
print(`all interfaces are allowed, except the 0 one. `):
print(`Using the Functional Equation.`):
print(`H=z*((1+H)^k-H)`):
print(` Try: `):
print(`KamaOneBHnzFe(3,5);`):


elif nops([args])=1 and op(1,[args])=KamaOneBPnzFe then
print(`KamaOneBPnzFe(k,N): The first N terms of the sequence:`):
print(`Number of Pyramids with N bricks, where only one kind of brick, `):
print(` of dimensions 1 by k participate and `):
print(`all interfaces are allowed, except the 0 one. `):
print(`Using the Functional Equation.`):
print(`H=z*((1+H)^k-k), P=H/(1-(k-1)*H)`):
print(` Try: `):
print(`KamaOneBPnzFe(3,5);`):


elif nops([args])=1 and op(1,[args])=KamaOneBMnzFe then
print(`KamaOneBMnzFe(k,N): The first N terms of the sequence:`):
print(`Number of Pyramids with N bricks, where only one kind of brick, `):
print(` of dimensions 1 by k participate and `):
print(`all interfaces are allowed, except the 0 one. `):
print(`Using the Functional Equation.`):
print(`H=z*((1+H)^k-k), M=H/(1-(k-1)*H)/(1-H)`):
print(` Try: `):
print(`KamaOneBMnzFe(3,5);`):

elif nops([args])=1 and op(1,[args])=KamaOneBMnz then
print(`KamaOneBMnz(k,N): The first N terms of the sequence:`):
print(`Number of towers with N bricks, where only one kind of brick, `):
print(` of dimensions 1 by k participate and `):
print(`all interfaces are allowed, except the 0 one. Try:`):
print(`KamaOneBMnz(3,5);`):

elif nops([args])=1 and op(1,[args])=KamaOneBPnz then
print(`KamaOneBPnz(k,N): The first N terms of the sequence:`):
print(`Number of Pyramids with N bricks, where only one kind of brick, `):
print(` of dimensions 1 by k participate and `):
print(`all interfaces are allowed, except the 0 one. Try:`):
print(`KamaOneBPnz(3,5);`):


elif nops([args])=1 and op(1,[args])=KamaManyBH then
print(`KamaManyBH(L,N): The first N terms of the sequence:`):
print(`Number of Half-Pyramids with N bricks, where the list of `):
print(` bricks is L and`):
print(`all interfaces are allowed, try:`):
print(`KamaManyBH([1,2],5);`):


elif nops([args])=1 and op(1,[args])=KamaManyBHnz then
print(`KamaManyBHnz(L,N): The first N terms of the sequence:`):
print(`Number of Half-Pyramids with N bricks, where the list of `):
print(` bricks is L and`):
print(`all interfaces are allowed, except 0. Try:`):
print(`KamaManyBHnz([1,2],5);`):




elif nops([args])=1 and op(1,[args])=KamaManyBP then
print(`KamaManyBP(L,N): The first N terms of the sequence:`):
print(`Number of Pyramids with N bricks, where the list of `):
print(` bricks is L and`):
print(`all interfaces are allowed, try:`):
print(`KamaManyBP([1,2],5);`):


elif nops([args])=1 and op(1,[args])=KamaManyBPnz then
print(`KamaManyBPnz(L,N): The first N terms of the sequence:`):
print(`Number of Pyramids with N bricks, where the list of `):
print(` bricks is L and`):
print(`all interfaces are allowed, except 0. Try:`):
print(`KamaManyBPnz([1,2],5);`):


elif nops([args])=1 and op(1,[args])=KamaManyBM then
print(`KamaManyBM(L,N): The first N terms of the sequence:`):
print(`Number of towers with N bricks, where the list of `):
print(` bricks is L and`):
print(`all interfaces are allowed, try:`):
print(`KamaManyBM([1,2],5);`):


elif nops([args])=1 and op(1,[args])=KamaManyBMnz then
print(`KamaManyBMnz(L,N): The first N terms of the sequence:`):
print(`Number of towers with N bricks, where the list of `):
print(` bricks is L and`):
print(`all interfaces are allowed, except 0. Try:`):
print(`KamaManyBMnz([1,2],5);`):

elif nops([args])=1 and op(1,[args])=KamaP then
print(`KamaP(L,M,N): Given a list of bricks (set of pos. integers) L`):
print(`and a list-of-lists of sets M, outputs a list of length N`):
print(`whose i-th item is the NUMBER of Pyramids of size i`):
print(`built from the bricks L using the interface rules M.`):
print(`For example, try:`):
print(`KamaP([2],[[{-1,1}]],10);`);
print(`KamaP([1,2],[ [{0},{-1,0}],[{0,1},{-1,0,1}]],10 );`):

elif nops([args])=1 and op(1,[args])=KamaM then
print(`KamaM(L,M,N): Given a list of bricks (set of pos. integers) L`):
print(`and a list-of-lists of sets M, outputs a list of length N`):
print(`whose i-th item is the NUMBER of towers of size i`):
print(`built from the bricks L using the interface rules M.`):
print(`For example, try:`):
print(`KamaM([2],[[{-1,1}]],10);`);
print(`KamaM([1,2],[ [{0},{-1,0}],[{0,1},{-1,0,1}]],10 );`):

elif nops([args])=1 and op(1,[args])=ManyBH then
print(`ManyBH(L,N): The first N terms of the sequence:`):
print(`SET of  Half-Pyramids with N bricks, where the list of `):
print(` bricks is L and`):
print(`all interfaces are allowed, try:`):
print(`ManyBH([1,2],5);`):

elif nops([args])=1 and op(1,[args])=ManyBHnz then
print(`ManyBHnz(L,N): The first N terms of the sequence:`):
print(`SET of  Half-Pyramids with N bricks, where the list of `):
print(` bricks is L and`):
print(`all interfaces are allowed EXCEPT 0. Try:`):
print(`ManyBHnz([1,2],5);`):

elif nops([args])=1 and op(1,[args])=ManyBP then
print(`ManyBP(L,N): The first N terms of the sequence:`):
print(`SET of  Pyramids with N bricks, where the list of `):
print(` bricks is L and`):
print(`all interfaces are allowed, try:`):
print(`ManyBP([1,2],5);`):

elif nops([args])=1 and op(1,[args])=ManyBPnz then
print(`ManyBPnz(L,N): The first N terms of the sequence:`):
print(`SET of  Pyramids with N bricks, where the list of `):
print(` bricks is L and`):
print(`all interfaces are allowed, EXCEPT 0. try:`):
print(`ManyBPnz([1,2],5);`):

elif nops([args])=1 and op(1,[args])=ManyBM then
print(`ManyBM(L,N): The first N terms of the sequence:`):
print(`SET of  Towers with N bricks, where the list of `):
print(` bricks is L and`):
print(`all interfaces are allowed, try:`):
print(`ManyBM([1,2],5);`):

elif nops([args])=1 and op(1,[args])=ManyBMnz then
print(`ManyBMnz(L,N): The first N terms of the sequence:`):
print(`SET of  Towers with N bricks, where the list of `):
print(` bricks is L and`):
print(`all interfaces are allowed, EXCEPT 0. try:`):
print(`ManyBMnz([1,2],5);`):


elif nops([args])=1 and op(1,[args])=Migdalim then
print(`Migdalim(L,M,N): Given a list of bricks (set of pos. integers) L`):
print(`and a list-of-lists of sets M, outputs a list of length N`):
print(`whose i-th item is the set of towers of size i`):
print(`built from the bricks L using the interface rules M`):
print(`For example, try:`):
print(`Migdalim([2],[[{-1,1}]],5);`):
print(`Migdalim([1,2],[ [{0},{-1,0}],[{0,1},{-1,0,1}]],5 );`):

elif nops([args])=1 and op(1,[args])=Mishkal then
print(` Mishkal(S,z): the sum of the weights of the towers of bricks  `):
print(`in the set S`):


elif nops([args])=1 and op(1,[args])=Pyramidot then
print(`Pyramidot(L,M,N): Given a list of bricks (set of pos. integers) L`):
print(`and a list-of-lists of sets M, outputs a list of length N`):
print(`whose i-th item is the set of PYRAMIDS of size i`):
print(`built from the bricks L using the interface rules M`):
print(`A pyramid is a tower whose first floor only has ONE brick`):
print(`For example, try:`):
print(`Pyramidot([2],[[{-1,1}]],5);`):
print(`Pyramidot([1,2],[ [{0},{-1,0}],[{0,1},{-1,0,1}]],5 );`):


elif nops([args])=1 and op(1,[args])=SeqH then
print(`SeqH(L,N): the first N terms of the sequence of`):
print(`whose n-th term is`):
print(` "number of half-pyramids with area n"`):
print(`(with all interfaces allowed) with bricks of sizes`):
print(`in L. Try:`):
print(`SeqH([1,4,5],10);`):

elif nops([args])=1 and op(1,[args])=SeqM then
print(`SeqM(L,N): the first N terms of the sequence of`):
print(`whose n-th term is`):
print(` "number of towers with area n"`):
print(`(with all interfaces allowed) with bricks of sizes`):
print(`in L. Try:`):
print(`SeqM([1,4,5],10);`):

elif nops([args])=1 and op(1,[args])=SeqP then
print(`SeqP(L,N): the first N terms of the sequence of`):
print(`whose n-th term is`):
print(` "number of pyramids with area n"`):
print(`(with all interfaces allowed) with bricks of sizes`):
print(`in L. Try:`):
print(`SeqP([1,4,5],10);`):

elif nops([args])=1 and op(1,[args])=SetH12nz then
print(`SetH12nz(N): the set of half-pyramids with tiles of length {1,2}`):
print(`without zero-interface of size n. Should be the same`):
print(`ManyBHnz([1,2],N)[N]. Try:`):
print(`SetH12nz(14);`):

elif nops([args])=1 and op(1,[args])=ShiftR then
print(`ShiftR(M,a): shifting the tower M by a units to the right`):
print(` try: `):
print(` ShiftR([[[0,3]]],2); `):

elif nops([args])=1 and op(1,[args])=Size1 then
print(`Size1(T): the number of cells in the tower T`):
print(`For example, try:`):
print(`Size1([[[0,2],[2,4]],[[1,3],[3,6]]] );`):

elif nops([args])=1 and op(1,[args])=TouchingBricks then
print(`TouchingBricks(Koma,B1): inputs a floor, Koma, and`):
print(`a potential brick B1 of the form [i,j]`):
print(`returns the list of bricks in Koma that`):
print(`have a non-empty intersection with B1. Try:`):
print(`TouchingBricks([[-2,2],[2,4],[4,5]],[1,6]);`):

elif nops([args])=1 and op(1,[args])=Wt then
print(` Wt(L,z): the weight of a tower of bricks  `):
print(` (the product of the weights of all bricks there) `):


elif nops([args])=1 and op(1,[args])=WtManyBH then
print(`WtManyBH(L,N,z): The list of length N whose i-th item`):
print(`is the sum of weights of the`):
print(`members of the set of half-pyramids with brick list given by L`):
print(`where all interfaces are allowed, try:`):
print(`WtManyBH([1,2],5,z);`):

elif nops([args])=1 and op(1,[args])=WtManyBHnz then
print(`WtManyBHnz(L,N,z): The list of length N whose i-th item`):
print(`is the sum of weights of the`):
print(`members of the set of half-pyramids with brick list given by L`):
print(`where all interfaces are allowed, EXCEPT 0. Try:`):
print(`WtManyBHnz([1,2],5,z);`):


elif nops([args])=1 and op(1,[args])=WtManyBP then
print(`WtManyBP(L,N,z): The list of length N whose i-th item`):
print(`is the sum of weights of the`):
print(`members of the set of pyramids with brick list given by L`):
print(`where all interfaces are allowed, try:`):
print(`WtManyBP([1,2],5,z);`):

elif nops([args])=1 and op(1,[args])=WtManyBPnz then
print(`WtManyBPnz(L,N,z): The list of length N whose i-th item`):
print(`is the sum of weights of the`):
print(`members of the set of pyramids with brick list given by L`):
print(`where all interfaces are allowed, EXCEPT 0. Try:`):
print(`WtManyBPnz([1,2],5,z);`):

elif nops([args])=1 and op(1,[args])=WtManyBM then
print(`WtManyBM(L,N,z): The list of length N whose i-th item`):
print(`is the sum of weights of the`):
print(`members of the set of TOWERS with brick list given by L`):
print(`where all interfaces are allowed, try:`):
print(`WtManyBM([1,2],5,z);`):

elif nops([args])=1 and op(1,[args])=WtManyBMnz then
print(`WtManyBMnz(L,N,z): The list of length N whose i-th item`):
print(`is the sum of weights of the`):
print(`members of the set of TOWERS with brick list given by L`):
print(`where all interfaces are allowed, EXCEPT 0. Try:`):
print(`WtManyBMnz([1,2],5,z);`):

elif nops([args])=1 and op(1,[args])=Zinn then
print(`Zinn(L): given a list of integers L`):
print(`outputs theta,mu, such that L[n] is roughly`):
print(`mu^n*n^theta.`):
prin(`Try:`):
print(`Zinn([seq(binomial(2*i,i),i=1..200)]);`):

elif nops([args])=1 and op(1,[args])=Zset then
print(`Zset(N): the set of half-pyramids with tiles of length {1,2}`):
print(`without zero-interface of size n. Should be the same`):
print(`ManyBHnz([1,2],N)[N]. Try:`):
print(`Zset(14);`):

else
print(`There is no ezra for`,args):
fi:
 
end:


###start stuff from Findrec
sn:=proc(resh,n):
-1/log(op(n+1,resh)*op(n-1,resh)/op(n,resh)^2):
#evalf(%):
end:

Zinn:=proc(resh)
local s1,s2,n,ka:
if nops(resh)<=10 then
 RETURN(FAIL):
fi:

if member(0,{op(nops(resh)-4..nops(resh),resh)}) then
  RETURN(FAIL):
fi:

ka:=abs(resh[nops(resh)]/resh[nops(resh)-1]-
resh[nops(resh)-1]/resh[nops(resh)-2]):
if ka <1/10^7 then
RETURN(0,evalf(resh[nops(resh)]/resh[nops(resh)-1])):
fi:
n:=nops(resh)-2:
s1:=sn(resh,n):
s2:=sn(resh,n-1):
evalf(2*(s1+s2)/(s1-s2)^2),
evalf(sqrt(op(n+1,resh)/op(n-1,resh))*exp(-(s1+s2)/((s1-s2)*s1))):
end:


###Findrec 
#findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating
#the sequence f of degree DEGREE and order ORDER
#For example, try: findrec([seq(i,i=1..10)],0,2,n,N);
findrecVerbose:=proc(f,DEGREE,ORDER,n,N)
local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a:
if (1+DEGREE)*(1+ORDER)+3+ORDER>nops(f) then
ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE):
fi:
ope:=0:
var:={}:
 
for i from 0 to ORDER do
 for j from 0 to DEGREE do
  ope:=ope+a[i,j]*n^j*N^i:
  var:=var union {a[i,j]}:
 od:
od:
    
eq:={}:
 
for n0 from 1 to (1+DEGREE)*(1+ORDER)+2 do
  eq1:=0:
 
  for i from 0 to ORDER do
     eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f):
  od:
 
   eq:= eq union {eq1}:
od:
 
var1:=solve(eq,var):
 
kv:={}:
 
for i from 1 to nops(var1) do
 mu:=op(i,var1):
 
if op(1,mu)=op(2,mu) then
   kv:= kv union {op(1,mu)}:
 fi:
 
od:

ope:=subs(var1,ope):

if ope=0 then
  RETURN(FAIL):
fi:

ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}:

if nops(ope)>1 then
print(`There is some slack, there are `, nops(ope)):
print(ope):
RETURN(Yafe(ope[1],N)[2]):
elif nops(ope)=1 then
RETURN(Yafe(ope[1],N)[2]):
else
 RETURN(FAIL):
fi:

end:
 
#findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating
#the sequence f of degree DEGREE and order ORDER
#For example, try: findrec([seq(i,i=1..10)],0,2,n,N);
findrec:=proc(f,DEGREE,ORDER,n,N)
local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a:
option remember:


if not findrecEx(f,DEGREE,ORDER,ithprime(20)) then
 RETURN(FAIL):
fi:

if not findrecEx(f,DEGREE,ORDER,ithprime(40)) then
 RETURN(FAIL):
fi:

if not findrecEx(f,DEGREE,ORDER,ithprime(80)) then
 RETURN(FAIL):
fi:


if (1+DEGREE)*(1+ORDER)+5+ORDER>nops(f) then
ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE):
fi:
ope:=0:
var:={}:
 
for i from 0 to ORDER do
 for j from 0 to DEGREE do
  ope:=ope+a[i,j]*n^j*N^i:
  var:=var union {a[i,j]}:
 od:
od:
    
eq:={}:
 
for n0 from 1 to (1+DEGREE)*(1+ORDER)+4 do
  eq1:=0:
 
  for i from 0 to ORDER do
     eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f):
  od:
 
   eq:= eq union {eq1}:
od:
 
var1:=solve(eq,var):
 
kv:={}:
 
for i from 1 to nops(var1) do
 mu:=op(i,var1):
 
if op(1,mu)=op(2,mu) then
   kv:= kv union {op(1,mu)}:
 fi:
 
od:

ope:=subs(var1,ope):

if ope=0 then
  RETURN(FAIL):
fi:


if {seq(
expand(add(subs(n=n1,coeff(ope,N,i))*f[n1+i],i=0..ORDER)),
n1=1..nops(f)-ORDER)}<>{0} then
 RETURN(FAIL):
fi:

ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}:

if nops(ope)>1 then
RETURN(Yafe(ope[1],N)[2]):
elif nops(ope)=1 then
RETURN(Yafe(ope[1],N)[2]):
else
 RETURN(FAIL):
fi:

end:
 


Yafe:=proc(ope,N) local i,ope1,coe1,L: 
if ope=0 then
 RETURN(1,0):
fi:
ope1:=expand(ope):
L:=degree(ope1,N):
coe1:=coeff(ope1,N,L):
ope1:=normal(ope1/coe1):
ope1:=normal(ope1):
ope1:=
convert(
[seq(factor(coeff(ope1,N,i))*N^i,i=ldegree(ope1,N)..degree(ope1,N))],`+`):
factor(coe1),ope1:
end:


#Findrec(f,n,N,MaxC): Given a list f tries to find a linear recurrence equation with
#poly coffs.
#of maximum DEGREE+ORDER<=MaxC
#e.g. try Findrec([1,1,2,3,5,8,13,21,34,55,89],n,N,2);
Findrec:=proc(f,n,N,MaxC)
local DEGREE, ORDER,ope,L:

for L from 0 to MaxC do
for ORDER from 0 to L do
 DEGREE:=L-ORDER:
if (2+DEGREE)*(1+ORDER)+4>=nops(f) then
 print(`Insufficient data for degree`, DEGREE, `and order `,ORDER):
 RETURN(FAIL):
fi:
 ope:=findrec([op(1..(2+DEGREE)*(1+ORDER)+4,f)],DEGREE,ORDER,n,N):

if  ope<>FAIL and
{seq(
expand(add(subs(n=n1,coeff(ope,N,i))*f[n1+i],i=0..ORDER)),
n1=1..nops(f)-ORDER)}={0} then
      RETURN(ope):
 fi:
 od:
od:
FAIL:

end:





 
#SeqFromRec(ope,n,N,Ini,K): Given the first L-1
#terms of the sequence Ini=[f(1), ..., f(L-1)]
#satisfied by the recurrence ope(n,N)f(n)=0
#extends it to the first K values
SeqFromRec:=proc(ope,n,N,Ini,K)
local ope1,gu,L,n1,j1:
ope1:=Yafe(ope,N)[2]:
L:=degree(ope1,N):
if nops(Ini)<>L then
 ERROR(`Ini should be of length`, L):
fi:

ope1:=expand(subs(n=n-L,ope1)/N^L):

gu:=Ini:

for n1 from nops(Ini)+1 to K do
gu:=[op(gu), -add(gu[nops(gu)+1-j1]*subs(n=n1,coeff(ope1,N,-j1)),
j1=1..L)]:
od:

gu:

end:
 
#end Findrec


with(linalg):

#findrecEx(f,DEGREE,ORDER,m1): Explores whether thre
#is a good chance that there is a recurrence of degree DEGREE
#and order ORDER, using the prime m1
#For example, try: findrecEx([seq(i,i=1..10)],0,2,n,N,1003);
findrecEx:=proc(f,DEGREE,ORDER,m1)
local ope,var,eq,i,j,n0,eq1,a,A1,
D1,E1,Eq,Var,f1,n,N:
option remember:
f1:=f mod m1:
if (1+DEGREE)*(1+ORDER)+5+ORDER>nops(f) then
ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE):
fi:
ope:=0:
var:={}:
 
for i from 0 to ORDER do
 for j from 0 to DEGREE do
  ope:=ope+a[i,j]*n^j*N^i:
  var:=var union {a[i,j]}:
 od:
od:
    
eq:={}:
 
for n0 from 1 to (1+DEGREE)*(1+ORDER)+4 do
  eq1:=0:
 
  for i from 0 to ORDER do
     eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f1) mod m1:
  od:
 
   eq:= eq union {eq1}:
od:


Eq:= convert(eq,list):
Var:= convert(var,list):


D1:=nops(Var):
E1:=nops(Eq):
if E1<D1 then
  RETURN(true):
fi:

A1:=matrix(D1,D1):

for i from 1 to D1-1 do
for j from 1 to D1 do
  A1[i,j]:=coeff(Eq[i],Var[j]):
od:
od:

for j from 1 to nops(Var) do
  A1[D1,j]:=coeff(Eq[D1],Var[j]):
od:

if det(A1) mod m1 <>0  then
 RETURN(false):
fi:

if E1-D1>=1 then
 for j from 1 to nops(Var) do
  A1[D1,j]:=coeff(Eq[D1+1],Var[j]):
 od:

if det(A1) mod m1 <>0  then
 RETURN(false):
fi:
fi:

if E1-D1>=2 then
 for j from 1 to nops(Var) do
  A1[D1,j]:=coeff(Eq[D1+2],Var[j]):
 od:

if det(A1) mod m1 <>0  then
 RETURN(false):
fi:
fi:

true:

end:
###End stuff from Findrec

#CheckBricks(L,M): given a list L, of positive integers
#and a list of length nops(L) of lists of of length nops(L)
#checks whether it is a valid input for buliding bricks
#the condition is that M[i][j] should be
#a subset of {-({L[j]-1)...(L[i]-1)} indicating the allowed
#interfaces of placing brick L[j] on TOP of brick L[i]
#For example, try:
#CheckBricks([2],[[{-1,1}]]);
CheckBricks:=proc(L,M) local i,j,a,b,k:

if not (type(L,list) and type(M,list) and 
{seq(type(M[i],list),i=1..nops(M))}={true}) then
 RETURN(false):
fi:

if nops(L)<>nops(M) then
 RETURN(false):
fi:

if {seq(nops(M[i]),i=1..nops(M))}<>{nops(L)} then
 RETURN(false):
fi:

if {seq(type(L[i],integer),i=1..nops(L))}<>{true} then
 RETURN(false):
fi:

if {seq(evalb(L[i]>0), i=1..nops(L))}<>{true} then
 RETURN(false):
fi:


for i from 1 to nops(L) do
   a:=L[i]: 
 for j from 1 to nops(L) do
  b:=L[j]:

  if M[i][j] minus {seq(k,k=-(b-1)..(a-1))}<>{} then
   RETURN(false):
  fi:
 od:
od:

true:

end:


#Size1(T): the number of cells in the tower T
#For example, try:
#Size1([[[0,2],[2,4]],[[1,3],[3,6]]);
Size1:=proc(T) local i,j:
add(
add(T[i][j][2]-T[i][j][1],j=1..nops(T[i])),i=1..nops(T)):
end:



#FindBrick(Koma,x): Given a floor Koma, and an integer x
#finds which brick of Koma lies under it (if it is above)
#a brick, or if it is between two bricks it
#returns the gap. For example, try:
#FindBrick([[-2,3],[5,7]], -1);

FindBrick:=proc(Koma,x) local i:

if x<Koma[1][1] then
RETURN(1/2):
fi:

for i from 1 to nops(Koma)-1 do

 if x>=Koma[i][1] and x<Koma[i][2] then
   RETURN(i):
 fi:

 if x>=Koma[i][2] and x<Koma[i+1][1] then
   RETURN(i+1/2):
 fi:

od:


i:=nops(Koma):

 if x>=Koma[i][1] and x<Koma[i][2] then
   RETURN(i):
 fi:

if x>=Koma[i][2] then
 RETURN(i+1/2):
fi:



end:

#TouchingBricks(Koma,B1): inputs a floor, Koma, and
#a potential brick B1 of the form [i,j]
#returns the list of bricks in Koma that
#have a non-empty intersection with B1. Try:
#TouchingBricks([[-2,2],[2,4],[4,5]],[1,6]);
TouchingBricks:=proc(Koma,B1) local gu,i,j:

gu:=[]:

for i from 1 to nops(Koma) do
 if {seq(j,j=Koma[i][1]+1..Koma[i][2])} intersect 
     {seq(j,j=B1[1]+1..B1[2])}<>{} then
     gu:=[op(gu),Koma[i]]:
 fi:
od:

gu:

end:

#IsComp(L,M,B1,B2): Is it legal to place brick B2 on top
#of brick B1, with the list of bricks L and 
#rules given by the list-of-lists M? Try:
#IsComp([1,2],[ [{0},{-1,0}],[{0,1},{-1,0,1}]],[4,6],[5,7]);
IsComp:=proc(L,M,B1,B2) local s1,s2,i1,i2,gu:
s1:=B1[2]-B1[1]:

if not member(s1, {op(L)}) then
 RETURN(false):
fi:

s2:=B2[2]-B2[1]:

if not member(s2, {op(L)}) then
 RETURN(false):
fi:

for i1 from 1 to nops(L) while s1<>L[i1] do od:

for i2 from 1 to nops(L) while s2<>L[i2] do od:

gu:=M[i1][i2]:

if member(B2[1]-B1[1],gu) then
  RETURN(true):
else
 RETURN(false):
fi:

end:



#ContSF(L,M,T,N): Given a list of bricks (set of pos. integers) L
#and a list-of-lists of sets M, describing the construction
#rules for building towers, and a tower T and a positive integer N
#outputs the set of towers obtained by adding ONE brick
#to the current top floor of T, whose size is still <=N .
#For example, try:
#ContSF([1,2],[ [{0},{-1,0}],[{0,1},{-1,0,1}]], [ [[0,2],[2,3],[3,5]], [[1,3]]],20 );
ContSF:=proc(L,M,T,N) local i, gu,sof,gu1,hatkh,B2,B1,Koma,lu,a:

if not CheckBricks(L,M) then
 RETURN(FAIL):
fi:

sof:=T[nops(T)][nops(T[nops(T)])][2]:

gu:={}:

if nops(T)=1 then

for i from 1 to nops(L) do
 gu1:=[[op(T[1]),[sof,sof+L[i]] ] ]:
  if Size1(gu1)<=N then
    gu:=gu union { gu1 }:
  fi:
od:

RETURN(gu):

fi:


hatkh:=sof:
Koma:=T[nops(T)-1]:
sof:=Koma[nops(Koma)][2]:
gu:={}:

for a from hatkh to sof-1 do
 for i from 1 to nops(L) do
    B2:=[a,a+L[i]]:

    lu:=TouchingBricks(Koma,B2):

    
  if {seq(IsComp(L,M,B1,B2),B1 in lu)}={true} then
       gu1:=[op(1..nops(T)-1,T),[op(T[nops(T)]),B2]]:

     if Size1(gu1)<=N then
       gu:=gu union {gu1}:
     fi:
  fi:
 od:
od:

gu:  

end:





#ContNF(L,M,T,N): Given a list of bricks (set of pos. integers) L
#and a list-of-lists of sets M, describing the construction
#rules for building towers, and a tower T and a positive integer N
#outputs the set of towers obtained by creating a new floor
#by adding ONE brick
#For example, try:
#ContNF([1,2],[ [{0},{-1,0}],[{0,1},{-1,0,1}]], [ [[0,2],[2,3],[3,5]], [[1,3]]],20 );
ContNF:=proc(L,M,T,N) local sof,hatkh,i,j,smol,B1,x,Koma,mu,mu1,gu:


if not CheckBricks(L,M) then
 RETURN(FAIL):
fi:

smol:=min(seq(seq(op(M[i][j]),j=1..nops(M[i])),i=1..nops(M))):


hatkh:=T[nops(T)][1][1]+smol:

sof:=T[nops(T)][nops(T[nops(T)])][2]:

gu:={}:

Koma:=T[nops(T)]:

for x from hatkh to sof do

 for i from 1 to nops(L) do

 B1:=[x,x+L[i]]:

 mu:=TouchingBricks(Koma,B1):


  if mu<>{} and {seq(IsComp(L,M,mu1,B1),mu1 in mu)}={true} 
    and Size1(T)+L[i]<=N then
       gu:=gu union {[op(T),[B1]]}:
  fi:
 od:
od:

gu:

end:


#Cont(L,M,T,N): Given a list of bricks (set of pos. integers) L
#and a list-of-lists of sets M, describing the construction
#rules for building towers, and a tower T and a positive integer N
#outputs the set of towers obtained by adding one extra brick
#either to the current top floor, or creating a new floor
#For example, try:
#Cont([1,2],[ [{0},{-1,0}],[{0,1},{-1,0,1}]], [ [[0,2],[2,3],[3,5]], [[1,3]]],20 );
Cont:=proc(L,M,T,N):

if not CheckBricks(L,M) then
 RETURN(FAIL):
fi:

ContSF(L,M,T,N) union ContNF(L,M,T,N):

end:



#Migdalim(L,M,N): Given a list of bricks (set of pos. integers) L
#and a list-of-lists of sets M, outputs a list of length N
#whose i-th item is the set of towers of size i
#built from the bricks L using the interface rules M
#For example, try:
#Migdalim([2],[[{-1,1}]],10);`);
#Migdalim([1,2],[ [{0},{-1,0}],[{0,1},{-1,0,1}]],10 );
Migdalim:=proc(L,M,N) local i,TA,gu,gu1:

if not CheckBricks(L,M) then RETURN(FAIL): fi:


for i from 1 to N do
 TA[i]:={}:
od:

 gu:={seq([[[0,L[i]]]],i=1..nops(L))}:

while gu<>{} do

for gu1 in gu do
 TA[Size1(gu1)]:=TA[Size1(gu1)] union {gu1}:
od:


gu:={seq(op(Cont(L,M,gu1,N)),gu1 in gu)}:

od:

[seq(TA[i],i=1..N)]:

end:



#KamaM(L,M,N): Given a list of bricks (set of pos. integers) L
#and a list-of-lists of sets M, outputs a list of length N
#whose i-th item is the NUMBER of towers of size i
#built from the bricks L using the interface rules M
#For example, try:
#KamaM([2],[[{-1,1}]],10);`);
#KamaM([1,2],[ [{0},{-1,0}],[{0,1},{-1,0,1}]],10 );
KamaM:=proc(L,M,N) local gu,i:

if not CheckBricks(L,M) then RETURN(FAIL): fi:

gu:=Migdalim(L,M,N):

[seq(nops(gu[i]),i=1..N)]:

end:


#Pyramidot(L,M,N): Given a list of bricks (set of pos. integers) L
#and a list-of-lists of sets M, outputs a list of length N
#whose i-th item is the set of pyramids of size i
#built from the bricks L using the interface rules M
#For example, try:
#Pyramidot([2],[[{-1,1}]],10);`);
#Pyramidot([1,2],[ [{0},{-1,0}],[{0,1},{-1,0,1}]],10 );
Pyramidot:=proc(L,M,N) local i,TA,gu,gu1:

if not CheckBricks(L,M) then RETURN(FAIL): fi:

for i from 1 to N do
 TA[i]:={}:
od:

 gu:={seq([[[0,L[i]]]],i=1..nops(L))}:

for gu1 in gu do
 TA[Size1(gu1)]:=TA[Size1(gu1)] union {gu1}:
od:

gu:={seq(op(ContNF(L,M,gu1,N)),gu1 in gu)}:

while gu<>{} do

for gu1 in gu do
 TA[Size1(gu1)]:=TA[Size1(gu1)] union {gu1}:
od:


gu:={seq(op(Cont(L,M,gu1,N)),gu1 in gu)}:

od:

[seq(TA[i],i=1..N)]:

end:

#KamaP(L,M,N): Given a list of bricks (set of pos. integers) L
#and a list-of-lists of sets M, outputs a list of length N
#whose i-th item is the NUMBER of pyramids of size i
#built from the bricks L using the interface rules M
#For example, try:
#KamaP([2],[[{-1,1}]],10);`);
#KamaP([1,2],[ [{0},{-1,0}],[{0,1},{-1,0,1}]],10 );
KamaP:=proc(L,M,N) local gu,i:

if not CheckBricks(L,M) then RETURN(FAIL): fi:

gu:=Pyramidot(L,M,N):

[seq(nops(gu[i]),i=1..N)]:

end:



#ContNFpos(L,M,T,N): Given a list of bricks (set of pos. integers) L
#and a list-of-lists of sets M, describing the construction
#rules for building towers, and a tower T and a positive integer N
#outputs the set of towers obtained by creating a new floor
#by adding ONE brick to the right of the y-axis
#For example, try:
#ContNFpos([1,2],[ [{0},{-1,0}],[{0,1},{-1,0,1}]], [ [[0,2],[2,3],[3,5]], [[1,3]]],20 );
ContNFpos:=proc(L,M,T,N) local sof,i,B1,x,Koma,mu,mu1,gu:


if not CheckBricks(L,M) then
 RETURN(FAIL):
fi:

sof:=T[nops(T)][nops(T[nops(T)])][2]:

gu:={}:

Koma:=T[nops(T)]:

for x from 0 to sof do

 for i from 1 to nops(L) do

 B1:=[x,x+L[i]]:

 mu:=TouchingBricks(Koma,B1):


  if mu<>{} and {seq(IsComp(L,M,mu1,B1),mu1 in mu)}={true} 
    and Size1(T)+L[i]<=N then
       gu:=gu union {[op(T),[B1]]}:
  fi:
 od:
od:

gu:

end:



#ContPos(L,M,T,N): Given a list of bricks (set of pos. integers) L
#and a list-of-lists of sets M, describing the construction
#rules for building towers, and a tower T and a positive integer N
#outputs the set of towers obtained by adding one extra brick
#either to the current top floor, or creating a new floor
#but always staying in the region x>=0
#For example, try:
#ContPos([1,2],[ [{0},{-1,0}],[{0,1},{-1,0,1}]], [ [[0,2],[2,3],[3,5]], [[1,3]]],20 );
ContPos:=proc(L,M,T,N):

if not CheckBricks(L,M) then
 RETURN(FAIL):
fi:

ContSF(L,M,T,N) union ContNFpos(L,M,T,N):

end:



#HPyramidot(L,M,N): Given a list of bricks (set of pos. integers) L
#and a list-of-lists of sets M, outputs a list of length N
#whose i-th item is the set of half-pyramids of size i
#built from the bricks L using the interface rules M
#For example, try:
#HPyramidot([2],[[{-1,1}]],10);`);
#HPyramidot([1,2],[ [{0},{-1,0}],[{0,1},{-1,0,1}]],10 );
HPyramidot:=proc(L,M,N) local i,TA,gu,gu1:

if not CheckBricks(L,M) then RETURN(FAIL): fi:

for i from 1 to N do
 TA[i]:={}:
od:

 gu:={seq([[[0,L[i]]]],i=1..nops(L))}:

for gu1 in gu do
 TA[Size1(gu1)]:=TA[Size1(gu1)] union {gu1}:
od:

gu:={seq(op(ContNFpos(L,M,gu1,N)),gu1 in gu)}:

while gu<>{} do

for gu1 in gu do
 TA[Size1(gu1)]:=TA[Size1(gu1)] union {gu1}:
od:


gu:={seq(op(ContPos(L,M,gu1,N)),gu1 in gu)}:

od:

[seq(TA[i],i=1..N)]:

end:

#KamaH(L,M,N): Given a list of bricks (set of pos. integers) L
#and a list-of-lists of sets M, outputs a list of length N
#whose i-th item is the NUMBER of half-pyramids of size i
#built from the bricks L using the interface rules M
#For example, try:
#KamaH([2],[[{-1,1}]],10);`);
#KamaH([1,2],[ [{0},{-1,0}],[{0,1},{-1,0,1}]],10 );
KamaH:=proc(L,M,N) local gu,i:

if not CheckBricks(L,M) then RETURN(FAIL): fi:
gu:=HPyramidot(L,M,N):

[seq(nops(gu[i]),i=1..N)]:

end:


#KamaOneBH(k,N): The first N terms of the sequence:
#Number of half-pyramids with one 1 by k brick where
#all interfaces are allowed, try:
#KamaOneBH(3,5);
KamaOneBH:=proc(k,N) local i,gu:
gu:=KamaH([k],[[{seq(i,i=-(k-1)..k-1)}]],N*k):
[seq(gu[k*i],i=1..N)]:
end:


#KamaOneBP(k,N): The first N terms of the sequence:
#Number of pyramids with one 1 by k brick where
#all interfaces are allowed, try:
#KamaOneBP(3,5);
KamaOneBP:=proc(k,N) local i,gu:
gu:=KamaP([k],[[{seq(i,i=-(k-1)..k-1)}]],N*k):
[seq(gu[k*i],i=1..N)]:
end:

#KamaOneBM(k,N): The first N terms of the sequence:
#Number of towers with one 1 by k brick where
#all interfaces are allowed, try:
#KamaOneBM(3,5);
KamaOneBM:=proc(k,N) local i,gu:
gu:=KamaM([k],[[{seq(i,i=-(k-1)..k-1)}]],N*k):
[seq(gu[k*i],i=1..N)]:
end:



#KamaOneBHnz(k,N): The first N terms of the sequence:
#Number of half-pyramids with one 1 by k brick where
#all interfaces are allowed, except for being right on top (i.e. 0) try:
#KamaOneBHnz(3,5);
KamaOneBHnz:=proc(k,N) local i,gu:
gu:=KamaH([k],[[{seq(i,i=-(k-1)..k-1)} minus {0}]],N*k):
[seq(gu[k*i],i=1..N)]:
end:


#KamaOneBPnz(k,N): The first N terms of the sequence:
#Number of pyramids with one 1 by k brick where
#all interfaces are allowed, except for being right on top (i.e. 0) try:
#try:
#KamaOneBPnz(3,5);
KamaOneBPnz:=proc(k,N) local i,gu:
gu:=KamaP([k],[[{seq(i,i=-(k-1)..k-1)} minus {0}]],N*k):
[seq(gu[k*i],i=1..N)]:
end:

#KamaOneBMnz(k,N): The first N terms of the sequence:
#Number of towers with one 1 by k brick where
#all interfaces are allowed, except for being right on top (i.e. 0) try:
#KamaOneBMnz(3,5);
KamaOneBMnz:=proc(k,N) local i,gu:
gu:=KamaM([k],[[{seq(i,i=-(k-1)..k-1)} minus {0}]],N*k):
[seq(gu[k*i],i=1..N)]:
end:





#KamaManyBH(L,N): The first N terms of the sequence:
#Number of half-pyramids with brick list given by L
#where all interfaces are allowed, try:
#KamaManyBH([1,2],5);
KamaManyBH:=proc(L,N) local M,M1,i,j,a,b,i1:
M:=[]:

for i from 1 to nops(L) do
 a:=L[i]:
 M1:=[]:
 
 for j from 1 to nops(L) do
  b:=L[j]:
  M1:=[op(M1),{seq(i1,i1=-(b-1)..a-1)}]:
 od:

M:=[op(M),M1]:

od:

KamaH(L,M,N):

end:


#KamaManyBHnz(L,N): The first N terms of the sequence:
#Number of half-pyramids with brick list given by L
#where all interfaces are allowed, except for 0. Try:
#KamaManyBHnz([1,2],5);
KamaManyBHnz:=proc(L,N) local M,M1,i,j,a,b,i1:
M:=[]:

for i from 1 to nops(L) do
 a:=L[i]:
 M1:=[]:
 
 for j from 1 to nops(L) do
  b:=L[j]:
  M1:=[op(M1),{seq(i1,i1=-(b-1)..a-1)} minus {0} ]:
 od:

M:=[op(M),M1]:

od:

KamaH(L,M,N):

end:



#KamaManyBP(L,N): The first N terms of the sequence:
#Number of pyramids with brick list given by L
#where all interfaces are allowed, try:
#KamaManyBP([1,2],5);
KamaManyBP:=proc(L,N) local M,M1,i,j,a,b,i1:
M:=[]:

for i from 1 to nops(L) do
 a:=L[i]:
 M1:=[]:
 
 for j from 1 to nops(L) do
  b:=L[j]:
  M1:=[op(M1),{seq(i1,i1=-(b-1)..a-1)}]:
 od:

M:=[op(M),M1]:

od:

KamaP(L,M,N):

end:


#KamaManyBPnz(L,N): The first N terms of the sequence:
#Number of half-pyramids with brick list given by L
#where all interfaces are allowed, except for 0. Try:
#KamaManyBPnz([1,2],5);
KamaManyBPnz:=proc(L,N) local M,M1,i,j,a,b,i1:
M:=[]:

for i from 1 to nops(L) do
 a:=L[i]:
 M1:=[]:
 
 for j from 1 to nops(L) do
  b:=L[j]:
  M1:=[op(M1),{seq(i1,i1=-(b-1)..a-1)} minus {0} ]:
 od:

M:=[op(M),M1]:

od:

KamaP(L,M,N):

end:



#KamaManyBM(L,N): The first N terms of the sequence:
#Number of towers with brick list given by L
#where all interfaces are allowed, try:
#KamaManyBM([1,2],5);
KamaManyBM:=proc(L,N) local M,M1,i,j,a,b,i1:
M:=[]:

for i from 1 to nops(L) do
 a:=L[i]:
 M1:=[]:
 
 for j from 1 to nops(L) do
  b:=L[j]:
  M1:=[op(M1),{seq(i1,i1=-(b-1)..a-1)}]:
 od:

M:=[op(M),M1]:

od:

KamaM(L,M,N):

end:


#KamaManyBMnz(L,N): The first N terms of the sequence:
#Number of towers with brick list given by L
#where all interfaces are allowed, except for 0. Try:
#KamaManyBMnz([1,2],5);
KamaManyBMnz:=proc(L,N) local M,M1,i,j,a,b,i1:
M:=[]:

for i from 1 to nops(L) do
 a:=L[i]:
 M1:=[]:
 
 for j from 1 to nops(L) do
  b:=L[j]:
  M1:=[op(M1),{seq(i1,i1=-(b-1)..a-1)} minus {0} ]:
 od:

M:=[op(M),M1]:

od:

KamaM(L,M,N):

end:




#KamaOneBHfe(k,N): The first N terms of the sequence:
#Number of half-pyramids with one 1 by k brick where
#all interfaces are allowed, 
#using the Functional Equation
#H=z*(1+H)^k
#try:
#KamaOneBHfe(3,5);
KamaOneBHfe:=proc(k,N) local gu,z,i,n:
gu:=z:

for  n from 1 to N do
gu:=expand(z*(1+gu)^k):
gu:=add(coeff(gu,z,i)*z^i,i=1..N):
od:

[seq(coeff(gu,z,i),i=1..N)]:

end:


#KamaOneBHnzFe(k,N): The first N terms of the sequence:
#Number of half-pyramids with one 1 by k brick where
#all interfaces are allowed, except 0
#using the Functional Equation
#H=z*((1+H)^k-H)
#try:
#KamaOneBHnzFe(3,5);
KamaOneBHnzFe:=proc(k,N) local gu,z,i,n:
gu:=z:

for  n from 1 to N do
gu:=expand(z*((1+gu)^k-gu)):
gu:=add(coeff(gu,z,i)*z^i,i=1..N):
od:

[seq(coeff(gu,z,i),i=1..N)]:

end:


#KamaOneBPfe(k,N): The first N terms of the sequence:
#Number of pyramids with one 1 by k brick where
#all interfaces are allowed, 
#using the Functional Equation
#H=z*(1+H)^k and P=H/(1-(k-1)*H)
#try:
#KamaOneBPfe(3,5);
KamaOneBPfe:=proc(k,N) local gu,z,i,n:

gu:=KamaOneBHfe(k,N):
gu:=add(gu[i]*z^i,i=1..N):
gu:=gu/(1-(k-1)*gu):
gu:=taylor(gu,z=0,N+1):
[seq(coeff(gu,z,i),i=1..N)]:

end:


#KamaOneBPnzFe(k,N): The first N terms of the sequence:
#Number of pyramids with one 1 by k brick where
#all interfaces are allowed except 0
#using the Functional Equation
#H=z*((1+H)^k-H) and P=H/(1-(k-1)*H)
#try:
#KamaOneBPnzFe(3,5);
KamaOneBPnzFe:=proc(k,N) local gu,z,i,n:

gu:=KamaOneBHnzFe(k,N):
gu:=add(gu[i]*z^i,i=1..N):
gu:=gu/(1-(k-1)*gu):
gu:=taylor(gu,z=0,N+1):
[seq(coeff(gu,z,i),i=1..N)]:

end:



#KamaOneBMfe(k,N): The first N terms of the sequence:
#Number of pyramids with one 1 by k brick where
#all interfaces are allowed, 
#using the Functional Equation
#H=z*(1+H)^k and M=H/(1-H)/(1-(k-1)*H)
#try:
#KamaOneBMfe(3,5);
KamaOneBMfe:=proc(k,N) local gu,z,i,n:

gu:=KamaOneBHfe(k,N):
gu:=add(gu[i]*z^i,i=1..N):
gu:=gu/(1-(k-1)*gu)/(1-gu):
gu:=taylor(gu,z=0,N+1):
[seq(coeff(gu,z,i),i=1..N)]:

end:


#KamaOneBMnzFe(k,N): The first N terms of the sequence:
#Number of towers with one 1 by k brick where
#all interfaces are allowed, 
#using the Functional Equation
#H=z*((1+H)^k-H) and M=H/(1-H)/(1-(k-1)*H)
#try:
#KamaOneBMnzFe(3,5);
KamaOneBMnzFe:=proc(k,N) local gu,z,i,n:

gu:=KamaOneBHnzFe(k,N):
gu:=add(gu[i]*z^i,i=1..N):
gu:=gu/(1-(k-1)*gu)/(1-gu):
gu:=taylor(gu,z=0,N+1):
[seq(coeff(gu,z,i),i=1..N)]:

end:





#ManyBH(L,N): The list of length N whose i-th item
#is the SET of half-pyramids with brick list given by L
#where all interfaces are allowed, try:
#ManyBH([1,2],5);
ManyBH:=proc(L,N) local M,M1,i,j,a,b,i1:
M:=[]:

for i from 1 to nops(L) do
 a:=L[i]:
 M1:=[]:
 
 for j from 1 to nops(L) do
  b:=L[j]:
  M1:=[op(M1),{seq(i1,i1=-(b-1)..a-1)}]:
 od:

M:=[op(M),M1]:

od:

HPyramidot(L,M,N):

end:


#ManyBP(L,N): The list of length N whose i-th item
#is the SET of pyramids with brick list given by L
#where all interfaces are allowed, try:
#ManyBP([1,2],5);
ManyBP:=proc(L,N) local M,M1,i,j,a,b,i1:
M:=[]:

for i from 1 to nops(L) do
 a:=L[i]:
 M1:=[]:
 
 for j from 1 to nops(L) do
  b:=L[j]:
  M1:=[op(M1),{seq(i1,i1=-(b-1)..a-1)}]:
 od:

M:=[op(M),M1]:

od:

Pyramidot(L,M,N):

end:


#ManyBM(L,N): The list of length N whose i-th item
#is the SET of Towers with brick list given by L
#where all interfaces are allowed, try:
#ManyBM([1,2],5);
ManyBM:=proc(L,N) local M,M1,i,j,a,b,i1:
M:=[]:

for i from 1 to nops(L) do
 a:=L[i]:
 M1:=[]:
 
 for j from 1 to nops(L) do
  b:=L[j]:
  M1:=[op(M1),{seq(i1,i1=-(b-1)..a-1)}]:
 od:

M:=[op(M),M1]:

od:

Migdalim(L,M,N):

end:

#Wt2(B,z): the weight of single a brick B=[a,b] being z[b-a]
Wt2:=proc(B,z) z[B[2]-B[1]]: end:

#Wt1(L,z): the weight of a row of bricks 
Wt1:=proc(L,z) local i:
mul(Wt2(L[i],z),i=1..nops(L)): 
end:

#Wt(L,z): the weight of a tower of bricks 
#(the product of the weights of all bricks there)
Wt:=proc(T,z) local i:mul(Wt1(T[i],z),i=1..nops(T)): end:

#Mishkal(S,z): the weight of a set of towers S in terms
#of the indexed variable z
Mishkal:=proc(S,z) local s:
add(Wt(s,z), s in S):
end:



#WtManyBH(L,N,z): The list of length N whose i-th item
#is the sum of weights of the
#members of the set of half-pyramids with brick list given by L
#where all interfaces are allowed, try:
#WtManyBH([1,2],5,z);
WtManyBH:=proc(L,N,z) local gu,i:
gu:=ManyBH(L,N,z):
[seq(Mishkal(gu[i],z),i=1..nops(gu))]:
end:

#WtManyBHnz(L,N,z): The list of length N whose i-th item
#is the sum of weights of the
#members of the set of half-pyramids with brick list given by L
#where all interfaces are allowed, EXCEPT 0. Try:
#WtManyBHnz([1,2],5,z);
WtManyBHnz:=proc(L,N,z) local gu,i:
gu:=ManyBHnz(L,N,z):
[seq(Mishkal(gu[i],z),i=1..nops(gu))]:
end:


#WtManyBP(L,N,z): The list of length N whose i-th item
#is the sum of weights of the
#members of the set of pyramids with brick list given by L
#where all interfaces are allowed, try:
#WtManyBP([1,2],5,z);
WtManyBP:=proc(L,N,z) local gu,i:
gu:=ManyBP(L,N,z):
[seq(Mishkal(gu[i],z),i=1..nops(gu))]:
end:

#WtManyBPnz(L,N,z): The list of length N whose i-th item
#is the sum of weights of the
#members of the set of pyramids with brick list given by L
#where all interfaces are allowed, EXCEPT 0. Try:
#WtManyBPnz([1,2],5,z);
WtManyBPnz:=proc(L,N,z) local gu,i:
gu:=ManyBPnz(L,N,z):
[seq(Mishkal(gu[i],z),i=1..nops(gu))]:
end:

#WtManyBM(L,N,z): The list of length N whose i-th item
#is the sum of weights of the
#members of the set of towers with brick list given by L
#where all interfaces are allowed, try:
#WtManyBM([1,2],5,z);
WtManyBM:=proc(L,N,z) local gu,i:
gu:=ManyBM(L,N,z):
[seq(Mishkal(gu[i],z),i=1..nops(gu))]:
end:

#WtManyBMnz(L,N,z): The list of length N whose i-th item
#is the sum of weights of the
#members of the set of towers with brick list given by L
#where all interfaces are allowed, EXCEPT 0. Try:
#WtManyBMnz([1,2],5,z);
WtManyBMnz:=proc(L,N,z) local gu,i:
gu:=ManyBMnz(L,N,z):
[seq(Mishkal(gu[i],z),i=1..nops(gu))]:
end:

#H12Bab(z,N): the first N terms of the sequence of
#polynomials "weight-enumerator of half-pyramids
#(with all interfaces allowed) with bricks of size
#1 and 2. Using Babylonian knowledge. Try:
#H12Bab(z,10);
H12Bab:=proc(z,N) local t,F,gu,i:

gu:=solve(F-t*z[1]-t^2*z[2]-(z[1]*t+2*z[2]*t^2)*F-z[2]*t^2*F^2,F)[2]:

gu:=taylor(gu,t=0,N+5):
[seq(expand(coeff(gu,t,i)),i=1..N)]:
end:



Chop1:=proc(f,t,n) local i:
 add(coeff(f,t,i)*t^i,i=0..n):
end:

#ITP(L,F,t,N): Given a list of polynomials L
#(let's n=nops(L))
#in the variable t and the indexed variable
#F[1], ..., F[N], describing a polynomial transformation
#of formal power series. Iterates it N times. For
#example try
#ITP([z[1]*t+z[1]*t*H[1]+z[1]*t*H[2],z[2]*t^2+2*z[2]*t^2*(H[1]+H[2])+z[2]*t^2*(H[1]+H[2])^2],H,t,5);
ITP:=proc(L,F,t,N) local gu,n,i,j:
n:=nops(L):
gu:=[0$n]:

 for i from 1 to N do
  gu:=expand(subs({seq(F[j]=gu[j], j=1..n)},L)):
  gu:=[seq(Chop1(gu[j],t,i),j=1..nops(gu))]:
 od:

expand(gu):

end:




#P12(z,N): the first N terms of the sequence of
#polynomials "weight-enumerator of pyramids
#(with all interfaces allowed) with bricks of size
#1 and 2. Try:
#P12(z,10);
P12:=proc(z,N) local t,F,gu,i,lu:

gu:=solve(F-t*z[1]-t^2*z[2]-(z[1]*t+2*z[2]*t^2)*F-z[2]*t^2*F^2,F)[2]:

lu:=gu/(1+z[1]*t-(1-z[1]*t)*gu):

lu:=taylor(lu,t=0,N+5):
[seq(expand(coeff(lu,t,i)),i=1..N)]:
end:


#M12(z,N): the first N terms of the sequence of
#polynomials "weight-enumerator of towers
#(with all interfaces allowed) with bricks of size
#1 and 2. Try:
#M12(z,10);
M12:=proc(z,N) local t,F,gu,i,lu,ku:

gu:=solve(F-t*z[1]-t^2*z[2]-(z[1]*t+2*z[2]*t^2)*F-z[2]*t^2*F^2,F)[2]:

lu:=gu/(1+z[1]*t-(1-z[1]*t)*gu):

ku:=lu/(1-gu):

ku:=taylor(ku,t=0,N+5):
[seq(expand(coeff(ku,t,i)),i=1..N)]:
end:



#A12f(t): the explicit formulas for the 
#generating functions in t, for the sequences
#number of (Half-Pyramids, Pyramids,and TOWERS)
#with bricks {1,2} with
#all interfaces allowed
#Try A12f(t);
A12f:=proc(t) local F,gu,lu,ku:

gu:=solve(F-t-t^2-(t+2*t^2)*F-t^2*F^2,F)[2]:

lu:=normal(gu/(1+t-(1-t)*gu)):

ku:=normal(lu/(1-gu)):

[gu,lu,ku]:

end:





#H12ShortCut(z,N): the first N terms of the sequence of
#polynomials "weight-enumerator of half-pyramids
#(with all interfaces allowed) with bricks of size
#1 and 2. Try:
#H12ShortCut(z,10);
H12ShortCut:=proc(z,N) local t,H,gu,i:

gu:=ITP([z[1]*t+z[2]*t^2+(z[1]*t+2*z[2]*t^2)*H[1]+z[2]*t^2*H[1]^2],H,t,N)[1]:


[seq(coeff(gu,t,i),i=1..degree(gu,t))]:

end:


#H123i(z,N): the first N terms of the sequence of
#polynomials "weight-enumerator of half-pyramids
#(with all interfaces allowed) with bricks of size
#1, 2, and 3. Try:
#H123i(z,10);
H123i:=proc(z,N) local t,H,gu,i:

gu:=ITP(
[z[1]*t+z[2]*t^2+z[3]*t^3+
(z[1]*t+2*z[2]*t^2+3*z[3]*t^3)*H[1]+
(z[2]*t^2+3*z[3]*t^3)*H[1]^2+z[3]*t^3*H[1]^3],H,t,N)[1]:


[seq(coeff(gu,t,i),i=1..degree(gu,t))]:

end:



#Hk(k,z,N): the first N terms of the sequence of
#polynomials "weight-enumerator of half-pyramids
#(with all interfaces allowed) with bricks of size
#1, 2,..k. Try:
#Hk(3,z,10);
Hk:=proc(k,z,N) local t,H,gu,i:

gu:=ITP([1+add(z[i]*t^i*H[1]^i,i=1..k)],H,t,N+1)[1]:

[seq(coeff(gu,t,i),i=1..N)]:

end:


#HL(L,z,N): the first N terms of the sequence of
#polynomials "weight-enumerator of half-pyramids
#(with all interfaces allowed) with bricks of sizes
#in L. Try:
#HL([1,4,5],z,10);
HL:=proc(L,z,N) local t,H,gu,i:

gu:=ITP([1+add(z[L[i]]*t^L[i]*H[1]^L[i],i=1..nops(L))],H,t,N+1)[1]:

[seq(coeff(gu,t,i),i=1..N)]:

end:


#SeqH(L,N): the first N terms of the sequence of
#whose n-th term is
# "number of half-pyramids with area n"
#(with all interfaces allowed) with bricks of sizes
#in L. Try:
#SeqH([1,4,5],10);
SeqH:=proc(L,N) local t,H,gu,i:

gu:=ITP([1+add(t^L[i]*H[1]^L[i],i=1..nops(L))],H,t,N+1)[1]:

[seq(coeff(gu,t,i),i=1..N)]:

end:




#ManyBHnz(L,N): The first N terms of the sets of
#half-pyramids with brick list given by L
#where all interfaces are allowed, except for 0. Try:
#ManyBHnz([1,2],5);
ManyBHnz:=proc(L,N) local M,M1,i,j,a,b,i1:
M:=[]:

for i from 1 to nops(L) do
 a:=L[i]:
 M1:=[]:
 
 for j from 1 to nops(L) do
  b:=L[j]:
  M1:=[op(M1),{seq(i1,i1=-(b-1)..a-1)} minus {0} ]:
 od:

M:=[op(M),M1]:

od:

HPyramidot(L,M,N):
end:


#ManyBPnz(L,N): The first N terms of the sets of
#pyramids with brick list given by L
#where all interfaces are allowed, except for 0. Try:
#ManyBPnz([1,2],5);
ManyBPnz:=proc(L,N) local M,M1,i,j,a,b,i1:
M:=[]:

for i from 1 to nops(L) do
 a:=L[i]:
 M1:=[]:
 
 for j from 1 to nops(L) do
  b:=L[j]:
  M1:=[op(M1),{seq(i1,i1=-(b-1)..a-1)} minus {0} ]:
 od:

M:=[op(M),M1]:

od:

Pyramidot(L,M,N):
end:


#ManyBMnz(L,N): The first N terms of the sets of
#towers with brick list given by L
#where all interfaces are allowed, except for 0. Try:
#ManyBMnz([1,2],5);
ManyBMnz:=proc(L,N) local M,M1,i,j,a,b,i1:
M:=[]:

for i from 1 to nops(L) do
 a:=L[i]:
 M1:=[]:
 
 for j from 1 to nops(L) do
  b:=L[j]:
  M1:=[op(M1),{seq(i1,i1=-(b-1)..a-1)} minus {0} ]:
 od:

M:=[op(M),M1]:

od:

Migdalim(L,M,N):
end:



#H12(z,N): the first N terms of the sequence of
#polynomials "weight-enumerator of half-pyramids
#(with all interfaces allowed) with bricks of size
#1 and 2, without using the shortcut Try:
#H12(z,10);
H12:=proc(z,N) local t,H,gu,i:

gu:=ITP
([z[1]*t+z[1]*t*(H[1]+H[2]),z[2]*t^2+2*z[2]*t^2*(H[1]+H[2])+
z[2]*t^2*(H[1]+H[2])^2],
H,t,N):


[seq(coeff(gu[1]+gu[2],t,i),i=1..N)]:

end:



#H12nz(z,N): the first N terms of the sequence of
#polynomials "weight-enumerator of half-pyramids
#(with all interfaces allowed except 0) with bricks of size
#1 and 2, without using the shortcut Try:
#H12nz(z,10);
H12nz:=proc(z,N) local t,H,gu,i:

gu:=ITP
([z[1]*t,z[2]*t^2+z[2]*t^2*(H[1]+H[2])+
z[2]*t^2*(H[1]*H[2]+H[2]^2)],
H,t,N):


[seq(coeff(gu[1]+gu[2],t,i),i=1..N)]:

end:





#H12nzBab(z,N): the first N terms of the sequence of
#polynomials "weight-enumerator of half-pyramids
#(with all interfaces allowed except 0) with bricks of size
#1 and 2. Using Babylonian knowledge. Try:
#H12nzBab(z,10);
H12nzBab:=proc(z,N) local t,H2,gu,i:

gu:=solve(t^2*z[2]*H2^2+(t^3*z[1]*z[2]+t^2*z[2]-1)*H2
+t^2*z[2]+t^3*z[1]*z[2],H2)[2]:
gu:=gu+t*z[1]:

gu:=taylor(gu,t=0,N+5):
[seq(expand(coeff(gu,t,i)),i=1..N)]:
end:





#PfromHold(k,z,t,H): expressing the generating function
#of pyramids with bricks of size {1, ....,k}
#in terms of half-pyramids where H is the symbol for
#the generating function of half-pyramids
#Try:
#PfromHold(2,z,t,H);
PfromHold:=proc(k,z,t,H) local P,eq,var,i,j:
var:={seq(P[i],i=1..k)}:
eq:={seq(P[i]-t^i*z[i]*H^i*(1+add((j-1)*P[j],j=2..k)),i=1..k)}:
eq:=expand(eq):
var:=solve(eq,var):
1+normal(add(subs(var,P[i]),i=1..k)):
end:


#PfromH(k,z,t,H): expressing the generating function
#of pyramids with bricks of size {1, ....,k}
#in terms of half-pyramids where H is the symbol for
#the generating function of half-pyramids
#Try:
#PfromH(2,z,t,H);
PfromH:=proc(k,z,t,H) local i:
1+add(t^i*z[i]*H^i,i=1..k)/
(1-add((i-1)*t^i*z[i]*H^i,i=1..k)):
end:


#PfromHL(L,z,t,H): expressing the generating function
#of pyramids with bricks of size L
#in terms of half-pyramids where H is the symbol for
#the generating function of half-pyramids
#Try:
#PfromHL([1,2,3],z,t,H);
PfromHL:=proc(L,z,t,H) local i:
normal(
1+add(t^L[i]*z[L[i]]*H^L[i],i=1..nops(L))/
(1-add((L[i]-1)*t^L[i]*z[L[i]]*H^L[i],i=1..nops(L)))):
end:

#PfromHLt(L,t,H): expressing the generating function
#of pyramids with bricks of size L
#in terms of half-pyramids where H is the symbol for
#the generating function of half-pyramids (no z[i]'s)
#Try:
#PfromHLt([1,2,3],t,H);
PfromHLt:=proc(L,t,H) local i:
normal(
1+add(t^L[i]*H^L[i],i=1..nops(L))/
(1-add((L[i]-1)*t^L[i]*H^L[i],i=1..nops(L)))):
end:



#Pk(k,z,N): the first N terms of the sequence of
#polynomials "weight-enumerator of pyramids
#(with all interfaces allowed) with bricks of size
#1, 2,..k. Try:
#Pk(3,z,10);
Pk:=proc(k,z,N) local t,H,gu,i,mu,H1:

gu:=PfromH(k,z,t,H):


mu:=Hk(k,z,N):

H1:=1+add(mu[i]*t^i,i=1..N):

gu:=subs(H=H1,gu):

gu:=taylor(gu,t=0,N+1):

[seq(expand(coeff(gu,t,i)),i=1..N)]:

end:

#PL(L,z,N): the first N terms of the sequence of
#polynomials "weight-enumerator of pyramids
#(with all interfaces allowed) with bricks of sizes
#given by the list L. Try:
#PL([1,2,3],z,10);
PL:=proc(L,z,N) local t,H,gu,i,mu,H1:

gu:=PfromHL(L,z,t,H):


mu:=HL(L,z,N):

H1:=1+add(mu[i]*t^i,i=1..N):

gu:=subs(H=H1,gu):

gu:=taylor(gu,t=0,N+1):

[seq(expand(coeff(gu,t,i)),i=1..N)]:

end:



#Mk(k,z,N): the first N terms of the sequence of
#polynomials "weight-enumerator of towers
#(with all interfaces allowed) with bricks of size
#1, 2,..k. Try:
#Mk(3,z,10);
Mk:=proc(k,z,N) local gu,mu,P1,H1,X1,i,t:

mu:=Hk(k,z,N):

gu:=Pk(k,z,N):

H1:=add(mu[i]*t^i,i=1..N):
P1:=add(gu[i]*t^i,i=1..N):
X1:=P1/(1-H1):
X1:=taylor(X1,t=0,N+1):

[seq(expand(coeff(X1,t,i)),i=1..N)]:

end:


#ML(L,z,N): the first N terms of the sequence of
#polynomials "weight-enumerator of towers
#(with all interfaces allowed) with bricks of sizes
#given by the list L.  Try:
#ML([1,2,3],z,10);
ML:=proc(L,z,N) local gu,mu,P1,H1,X1,i,t:

mu:=HL(L,z,N):

gu:=PL(L,z,N):

H1:=add(mu[i]*t^i,i=1..N):
P1:=add(gu[i]*t^i,i=1..N):
X1:=P1/(1-H1):
X1:=taylor(X1,t=0,N+1):

[seq(expand(coeff(X1,t,i)),i=1..N)]:

end:

#CanPass(R,p): can the piece p pass through the row of bricks R?
#For example, try:
#CanPass([[0,2],[2,4]],[1,3]);
CanPass:=proc(R,p) local i,a,b:

if not (type(p,list) and nops(p)=2 and type(p[1],integer) and
       type(p[2],integer) ) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

a:=p[1]: b:=p[2]:

for i from 1 to nops(R) while R[i][2]<=a do od:
i:=i-1:
if i=nops(R) or b<=R[i+1][1] then
 true:
else
 false:
fi:

end:




#Land1(M,p): landing a piece p from heaven onto
#a tower M, getting a larger tower
#try:
#Land1([[[0,2]]],[0,1]);
Land1:=proc(M,p) local i,R,j:
for i from nops(M) by -1 to 1 while CanPass(M[i],p) do od:

if i=nops(M) then
RETURN([op(M),[p]]):
fi:

R:=M[i+1]:


for j from 1 to nops(R) while R[j][2]<=p[1] do od:

 R:=[op(1..j-1,R),p,op(j..nops(R),R)]:


[op(1..i,M),R,op(i+2..nops(M),M)]:


end:



#Land(M1,M2): landing a tower M2 from heaven onto
#a tower M1, getting a larger tower
#try:
#Land([[[0,2]]],[[[0,1]]]);
Land:=proc(M1,M2) local T,i,j:

T:=M1:

for i from 1 to nops(M2) do
 for j from 1 to nops(M2[i]) do
  T:=Land1(T,M2[i][j]):
 od:
od:

T:
end:



#ShiftR(M,a): shifting the tower M by a units to the right
#try:
#ShiftR([[[0,3]]],2);
ShiftR:=proc(M,a) local i,j:

[seq([seq([M[i][j][1]+a,M[i][j][2]+a],j=1..nops(M[i]))],i=1..nops(M))]:

end:





#Zset(N): the set of half-pyramids with tiles of length {1,2}
#without zero-interface of size n. Should be the same
#ManyBHnz([1,2],N)[N]. Try:
#Zset(14);
Zset:=proc(N) local gu,mu,mu1,lu,lu1,r:
option remember:

if N<2 then
 RETURN({}):
fi:

if N=2 then
 RETURN({[ [[0,2]] ]}):
fi:


if N=3 then
 RETURN({[ [[0,2]], [[1,2]] ]}):
fi:


gu:={}:

mu:=Zset(N-2):

for mu1 in mu do
 gu:=gu union {Land([[[0,2]]], ShiftR(mu1,1))}:
od:

mu:=Zset(N-3):

for mu1 in mu do
 gu:=gu union {Land([ [[0,2]],[[1,2]] ], mu1)}:
od:


for r from 2 to N-4 do
 mu:=Zset(r):
 lu:=Zset(N-r-2):

  for mu1 in mu do
   for lu1 in lu do
    gu:=gu union {Land(Land([ [[0,2]] ], ShiftR(mu1,1) ),lu1)  }:
   od:
 od:

od:

gu:

end:



#SeqP(L,N): the first N terms of the sequence whose n-th term is
#"the number of pyramids
#(with all interfaces allowed) with bricks of sizes
#given by the list L". Try:
#SeqP([1,2,3],10);
SeqP:=proc(L,N) local t,H,gu,i,mu,H1:

gu:=PfromHLt(L,t,H):


mu:=SeqH(L,N):

H1:=1+add(mu[i]*t^i,i=1..N):

gu:=subs(H=H1,gu):

gu:=taylor(gu,t=0,N+1):

[seq(expand(coeff(gu,t,i)),i=1..N)]:

end:

#SeqM(L,N): the first N terms of the sequence of
#polynomials "number of towers
#(with all interfaces allowed) with bricks of sizes
#given by the list L.  Try:
#SeqM([1,2,3],10);
SeqM:=proc(L,N) local gu,mu,P1,H1,X1,i,t:

mu:=SeqH(L,N):

gu:=SeqP(L,N):

H1:=add(mu[i]*t^i,i=1..N):
P1:=add(gu[i]*t^i,i=1..N):
X1:=P1/(1-H1):
X1:=taylor(X1,t=0,N+1):

[seq(expand(coeff(X1,t,i)),i=1..N)]:

end:


#empir(gu,degx,degP,x,P) 
#to "fit" an algebraic equation F(P(x),x)=0 of degree
#degP in P(x) and degx in n for P(x):=sum_i gu[i]*x^i
 
empir:=proc(gu,degx,degP,x,P)
local i1,i2,F,a,cand,lu,eq,var,mu,zak,mu1,yofee:
if (degx+1)*(degP+1)+3 > nops(gu) then
 RETURN(`sequence too small`):
fi:
 
F:=0:
var:={}:
for i1 from 0 to degx do
 for i2 from 0 to degP do
  F:=F+a[i1,i2]*x^i1*P^i2:
  var:=var union {a[i1,i2]}:
 od:
od:
 
cand:=0:
 
for i1 from 0 to (degx+1)*(degP+1)+3 do
 cand:=cand+op(i1+1,gu)*x^i1:
od:
 
lu:=subs(P=cand,F):
lu:=taylor(lu,x=0,nops(gu)-1):
 
eq:={}:
 
for i1 from 0 to (degx+1)*(degP+1)+3 do
eq:=eq union {coeff(lu,x,i1)=0}
od:
 
mu:=solve(eq,var):
 
zak:={}:

for mu1 in mu do
 if op(1,mu1)=op(2,mu1) then
   zak:=zak union {op(1,mu1)}:
 fi:
od:

 
F:=subs(mu,F):
 
if F=0 then
 RETURN(FAIL):
fi:

if nops(zak)>1 then
 RETURN(FAIL):
else
yofee:=coeff(F,zak[1],1):
 cand:=add(op(i1+1,gu)*x^i1,i1=0..nops(gu)-1):

 if ldegree(expand(subs(P=cand,yofee)),x)<nops(gu)-4 then
    RETURN(FAIL):
   else
    RETURN(yofee):
 fi:
 

fi:

end:
 

Empir:=proc(gu1,x,P)
local degx,degP,L,lu,gu:
 
gu:=[1,op(gu1)]:

for L from 1 to (nops(gu)-3)/3 do
for degP from 1 to L do
for degx from 0 to min(trunc((nops(gu)-3)/(1+degP))-1,L-degP) do
 
lu:=empir(gu,degx,degP,x,P):
 
if lu<>FAIL then
RETURN(collect(lu,P)):
fi:
od:
od:
od:
FAIL:
end:





#Meda(Li,x,P,n,N,MaxC,K): Given a sequence Li, outputs, if possible
#(i) a (conjectured, but often proved [depending on the context])
# algebraic equation satisfied by its ordinary generating function
#denoted by P(x)=1+add(L[i]*x^i,i=1..infinity) (but of course only
#using nops(Li) terms) (it gives FAIL if it fails)
#
#(ii) a (conjectured, but often routinely provable
# [depending on the context]) linear recurrence operator 
#in n and the shift-operator N annihilating the sequence
#of complexity <=MaxC
#
#(iii) Its K-th term (from the operator, if found) otherwise FAIL
#(iv) empirically derived (but once again provable) asymptotics
#For example, try:
#Meda([seq(binomial(2*i,i)/(i+1),i=1..50)],x,P,n,N,10,1000);

Meda:=proc(Li,x,P,n,N,MaxC,K) local gu,gu1,ope,LiG,LiG2,d1,d2,ju1,ju2,
C1a,C1b,d:
Digits:=1000:
gu1:=Empir(Li,x,P):

if gu1=FAIL then
 RETURN(FAIL):
fi:

gu:=[gu1]:

ope:=Findrec(Li,n,N, MaxC):



if ope=FAIL then
 RETURN(gu):
fi:

gu:=[op(gu),ope]:



LiG:=SeqFromRec(ope,n,N,[op(1..degree(ope,N),Li)],K):
gu:=[op(gu),LiG[K]]:


LiG2:=SeqFromRec(ope,n,N,[op(1..degree(ope,N),Li)],2*K):



ju1:=Zinn(LiG):

if ju1=FAIL then
 RETURN([op(gu),FAIL]):
fi:

ju2:=Zinn(LiG2):

if ju2=FAIL then
 RETURN([op(gu),FAIL]):
fi:

if not (type(ju1[1],numeric) and  type(ju1[2],numeric) 
and type(ju2[1],numeric) and type(ju2[2],numeric)) then
RETURN([op(gu), FAIL]):
fi:

if (ju1[1]=0  or ju2[1]=0) then
 RETURN([op(gu),[0,evalf(ju1[2],10)]]):
fi:

if ju1=ju2 then
gu:=[op(gu),[ju1]]:


gu:=[op(gu), FAIL]:

else

if abs(ju1[1])<0.01 then

if abs(abs(ju1[2]/ju2[2])-1)<0.01 then

d2:=-trunc(evalf(ln(abs(ju1[2]-ju2[2]))/ln(10)))-1:
gu:=[op(gu),[0,evalf(ju2[2],min(d2,10))]]:

else

gu:=[op(gu), FAIL]:

fi:
else

if abs(abs(ju1[1]/ju2[1])-1)<0.01 and 
  abs(abs(ju1[2]/ju2[2])-1)<0.01 then

d1:=-trunc(evalf(ln(abs(ju1[1]-ju2[1]))/ln(10)))-1:

d2:=-trunc(evalf(ln(abs(ju1[2]-ju2[2]))/ln(10)))-1:


gu:=[op(gu),[evalf(ju2[1],d1),evalf(ju2[2],min(d2,10))]]:
else

gu:=[op(gu), FAIL]:
fi:

fi:

fi:

C1a:=LiG[K]/K^gu[4][1]/gu[4][2]^K:
C1b:=LiG2[2*K]/(2*K)^gu[4][1]/gu[4][2]^(2*K):

if C1a<>C1b then
d:=trunc(-log(evalf(abs(C1a-C1b)))/log(10.));
else
d:=20:
fi:

if d>=2 then
gu:=[op(gu),evalf(C1b,min(d,10))]:
else
 gu:=[op(gu),FAIL]:
fi:
gu:
end:

 



#PolH(L,H,t): the polynomial, P, in H and t such that
#P(H(t),t)=0, where H(t) is the ordinary generating function
#for enumerating towers of k-mers of sizes drawn from L
#For example, try:
#PolH([1,3],H,t);
PolH:=proc(L,H,t) local k:
sort(1-H+add(t^L[k]*H^L[k],k=1..nops(L))):
end:







#Mishpat1(L,t,H,P,M,n,MaxC,K,Godel1,Godel2,co):
#inputs:
#(i) a list, L, of increasing positive integers
#(ii) a symbol t
#(iii, iv, v) symbols H,P,M (for ordinary generating functions)
#(vi) symbol n
#(vii) a positive integer MaxC (for the maximal complexity of the recurrence)
#(viii) a large positive integer (for the asympotics and for displaying
#a specific value)
#(ix,x) positive integers Godel1 and Godel2 (Godel1<Godel2)
#(xi): a pos. integer co, for labelling the Theorem
#Outputs: a theorem about the sequences enumerating
#half-pyramids, pyramids, and towers consisting of k-mers
#whose sizes are drawn from the list L.
#For example, try:
#Mishpat1([1,2],t,H,P,M,n,10,1000,30,100,1):
Mishpat1:=proc(L,t,H,P,M,n,MaxC,K,Godel1,Godel2,co) local guH,guP,guM,N,
luH,luP,luM,ope,C,i,g1:

g1:=igcd(op(L)):

guH:=SeqH(L,Godel2*g1):

guP:=SeqP(L,Godel2*g1):

guM:=SeqM(L,Godel2*g1):

guH:=[seq(guH[g1*i],i=1..Godel2)]:
guP:=[seq(guP[g1*i],i=1..Godel2)]:
guM:=[seq(guM[g1*i],i=1..Godel2)]:

print(`Theorem Number`, co):

print(`This theorem is about the enumeration of half-pyramids,`):
print(`pyramids, and towers using k-mers of sizes drawn from the set `):
print(convert(L,set)):

print(`Part 1: Half-Pyramids`):

print(`Let  h(n) be the number of ways of building`):
print(`half-pyramids of area`, n*g1, `using k-mers of sizes from the set`, 
convert(L,set)):

print(`The values of h(n) for n from 1 to`, Godel1, `are `):
print([op(1..Godel1,guH)]):


luH:=Meda(guH,t,H,n,N,MaxC,K):

if luH<>FAIL and nops(luH)>=1 then
print(`Let H(t) be the ordinary generating function of h(n), i.e.`):
print(H(t)=1+Sum(h(n)*t^n,n=1..infinity)):
print(`then, using the clever Betrema-Penaud argument, as extended`):
print(`by Doron Zeilberger , it is readily seen that H(t) satisfies`):
print(`the algebaric equation`):
print(luH[1]=0):
print(`and in Maple input format`):
lprint(luH[1]=0):
fi:

if  nops(luH)>=3 and luH[2]<>FAIL then
ope:=luH[2]:
print(`The sequence h(n) satisfies the following linear recurrence`):
print(`equation with polynomial coefficients`):
print(add(coeff(ope,N,i)*h(n+i),i=0..degree(ope,N))=0):
print(``):
print(`just for kicks, here is the value of `, h(K)):
print(luH[3]):
print(``):
fi:



if nops(luH)>=4 and luH[4]<>FAIL then
print(`h(n) is, apprx., asymototically equal to`):
print(C*n^luH[4][1]*luH[4][2]^n):
print(`for some constant C`):
fi:

if nops(luH)>=5 and luH[5]<>FAIL then
print(`that is approximately`):
print(luH[5]):
fi:


print(`Part 2: Pyramids`):

print(`Let  p(n) be the number of ways of building`):
print(`pyramids of area`, n*g1, `using k-mers of sizes from the set`, 
convert(L,set)):

print(`The values of p(n) for n from 1 to`, Godel1, `are `):
print([op(1..Godel1,guP)]):


luP:=Meda(guP,t,P,n,N,MaxC,K):



if luP<>FAIL and nops(luP)>=1 then
print(`Let P(t) be the ordinary generating function of p(n), i.e.`):
print(P(t)=1+Sum(p(n)*t^n,n=1..infinity)):
print(`then, using the clever Betrema-Penaud argument, as extended`):
print(`by Doron Zeilberger , it is readily seen that P(t) satisfies`):
print(`the algebaric equation`):
print(luP[1]=0):
print(`and in Maple input format`):
lprint(luP[1]=0):
fi:

if  nops(luP)>=3 and luP[2]<>FAIL then
ope:=luP[2]:
print(`The sequence p(n) satisfies the following linear recurrence`):
print(`equation with polynomial coefficients`):
print(add(coeff(ope,N,i)*p(n+i),i=0..degree(ope,N))=0):
print(``):
print(`just for kicks, here is the value of `, p(K)):
print(luP[3]):
print(``):
fi:



if nops(luP)>=4 and luP[4]<>FAIL then
print(`p(n) is, apprx., asymototically equal to`):
print(C*n^luP[4][1]*luP[4][2]^n):
print(`for some constant C`):
fi:

if nops(luP)>=5 and luP[5]<>FAIL then
print(`that is approximately`):
print(luP[5]):
fi:




print(`Part 3: Towers`):

print(`Let  m(n) be the number of ways of building`):
print(`TOWERS of area`, n*g1, `using k-mers of sizes from the set`, 
convert(L,set)):

print(`The values of m(n) for n from 1 to`, Godel1, `are `):
print([op(1..Godel1,guM)]):


luM:=Meda(guM,t,M,n,N,MaxC,K):



if luM<>FAIL and nops(luM)>=1 then
print(`Let M(t) be the ordinary generating function of m(n), i.e.`):
print(M(t)=1+Sum(m(n)*t^n,n=1..infinity)):
print(`then, using the clever Betrema-Penaud argument, as extended`):
print(`by Doron Zeilberger , it is readily seen that M(t) satisfies`):
print(`the algebaric equation`):
print(luM[1]=0):
print(`and in Maple input format`):
lprint(luM[1]=0):
fi:

if  nops(luM)>=3 and luM[2]<>FAIL then
ope:=luM[2]:
print(`The sequence m(n) satisfies the following linear recurrence`):
print(`equation with polynomial coefficients`):
print(add(coeff(ope,N,i)*m(n+i),i=0..degree(ope,N))=0):
print(``):
print(`just for kicks, here is the value of `, m(K)):
print(luM[3]):
print(``):
fi:



if nops(luM)>=4 and luM[4]<>FAIL then
print(`m(n) is, apprx., asymototically equal to`):
print(C*n^luM[4][1]*luM[4][2]^n):
print(`for some constant C`):
fi:

if nops(luM)>=5 and luM[5]<>FAIL then
print(`that is approximately`):
print(luM[5]):
fi:

print(`this ends this exciting theorem, whose proof is routine.`):

end:




#Mishpat2(L,t,H,P,M,n,MaxC,K,Godel1,Godel2,co):
#inputs:
#(i) a pos. integer L, 
#(ii) a symbol t
#(iii, iv, v) symbols H,P,M (for ordinary generating functions)
#(vi) symbol n
#(vii) a positive integer MaxC (for the maximal complexity of the recurrence)
#(viii) a large positive integer (for the asympotics and for displaying
#a specific value)
#(ix,x) positive integers Godel1 and Godel2 (Godel1<Godel2)
#(xi): a pos. integer co, for labelling the Theorem
#Outputs: a theorem about the sequences enumerating
#half-pyramids, pyramids, and towers consisting of L-mers
#without any restrictions about the interfaces
#For example, try:
#Mishpat2(2,t,H,P,M,n,10,1000,30,100,1):
Mishpat2:=proc(L,t,H,P,M,n,MaxC,K,Godel1,Godel2,co) local guH,guP,guM,N,
luH,luP,luM,ope,C,i:

guH:=SeqH([L],L*Godel2):

guH:=[seq(guH[i*L],i=1..Godel2)]:

guP:=SeqP([L],L*Godel2):

guP:=[seq(guP[i*L],i=1..Godel2)]:

guM:=SeqM([L],L*Godel2):

guM:=[seq(guM[i*L],i=1..Godel2)]:

print(`Theorem Number`, co):

print(`This theorem is about the enumeration of half-pyramids,`):
print(`pyramids, and towers using L-mers of the same size`, L):
print(`Where it is allowed to place a brick anywhere on top of`):
print(`a brick below it`):

print(`Part 1: Half-Pyramids`):

print(`Let  h(n) be the number of ways of building`):
print(`half-pyramids with n bricks, all of size`, L):
print(`where all the interfaces are allowed.`):

print(`The values of h(n) for n from 1 to`, Godel1, `are `):
print([op(1..Godel1,guH)]):


luH:=Meda(guH,t,H,n,N,MaxC,K):

if luH<>FAIL and nops(luH)>=1 then
print(`Let H(t) be the ordinary generating function of h(n), i.e.`):
print(H(t)=1+Sum(h(n)*t^n,n=1..infinity)):
print(`then, using the clever Betrema-Penaud argument, as extended`):
print(`by Doron Zeilberger , it is readily seen that H(t) satisfies`):
print(`the algebaric equation`):
print(luH[1]=0):
print(`and in Maple input format`):
lprint(luH[1]=0):
fi:

if  nops(luH)>=3 and luH[2]<>FAIL then
ope:=luH[2]:
print(`The sequence h(n) satisfies the following linear recurrence`):
print(`equation with polynomial coefficients`):
print(add(coeff(ope,N,i)*h(n+i),i=0..degree(ope,N))=0):
print(``):
print(`just for kicks, here it the value of `, h(K)):
print(luH[3]):
print(``):
fi:



if nops(luH)>=4 and luH[4]<>FAIL then
print(`h(n) is, apprx., asymototically equal to`):
print(C*n^luH[4][1]*luH[4][2]^n):
print(`for some constant C`):
fi:

if nops(luH)>=5 and luH[5]<>FAIL then
print(`that is approximately`):
print(luH[5]):
fi:


print(`Part 2: Pyramids`):

print(`Let  p(n) be the number of ways of building`):
print(`pyramids with n bricks, all of size`, L):
print(`where all the interfaces are allowed.`):


print(`The values of p(n) for n from 1 to`, Godel1, `are `):
print([op(1..Godel1,guP)]):


luP:=Meda(guP,t,P,n,N,MaxC,K):


if luP<>FAIL and nops(luP)>=1 then
print(`Let P(t) be the ordinary generating function of p(n), i.e.`):
print(P(t)=1+Sum(p(n)*t^n,n=1..infinity)):
print(`then, using the clever Betrema-Penaud argument, as extended`):
print(`by Doron Zeilberger , it is readily seen that P(t) satisfies`):
print(`the algebaric equation`):
print(luP[1]=0):
print(`and in Maple input format`):
lprint(luP[1]=0):
fi:

if  nops(luP)>=3 and luP[2]<>FAIL then
ope:=luP[2]:
print(`The sequence p(n) satisfies the following linear recurrence`):
print(`equation with polynomial coefficients`):
print(add(coeff(ope,N,i)*p(n+i),i=0..degree(ope,N))=0):
print(``):
print(`just for kicks, here it the value of `, p(K)):
print(luP[3]):
print(``):
fi:



if nops(luP)>=4 and luP[4]<>FAIL then
print(`p(n) is, apprx., asymototically equal to`):
print(C*n^luP[4][1]*luP[4][2]^n):
print(`for some constant C`):
fi:

if nops(luP)>=5 and luP[5]<>FAIL then
print(`that is approximately`):
print(luP[5]):
fi:




print(`Part 3: Towers`):

print(`Let  m(n) be the number of ways of building`):
print(`TOWERS with n bricks, all of size`, L):
print(`where all the interfaces are allowed.`):


print(`The values of m(n) for n from 1 to`, Godel1, `are `):
print([op(1..Godel1,guM)]):



luM:=Meda(guM,t,M,n,N,MaxC,K):



if luM<>FAIL and nops(luM)>=1 then
print(`Let M(t) be the ordinary generating function of m(n), i.e.`):
print(M(t)=1+Sum(m(n)*t^n,n=1..infinity)):
print(`then, using the clever Betrema-Penaud argument, as extended`):
print(`by Doron Zeilberger , it is readily seen that M(t) satisfies`):
print(`the algebaric equation`):
print(luM[1]=0):
print(`and in Maple input format`):
lprint(luM[1]=0):
fi:

if  nops(luM)>=3 and luM[2]<>FAIL then
ope:=luM[2]:
print(`The sequence m(n) satisfies the following linear recurrence`):
print(`equation with polynomial coefficients`):
print(add(coeff(ope,N,i)*m(n+i),i=0..degree(ope,N))=0):
print(``):
print(`just for kicks, here it the value of `, m(K)):
print(luM[3]):
print(``):
fi:




if nops(luM)>=4 and luM[4]<>FAIL then
print(`m(n) is, apprx., asymototically equal to`):
print(C*n^luM[4][1]*luM[4][2]^n):
print(`for some constant C`):
fi:

if nops(luM)>=5 and luM[5]<>FAIL then
print(`that is approximately`):
print(luM[5]):
fi:

print(`this ends this exciting theorem, whose proof is routine.`):

end:




#Mishpat3(L,t,H,P,M,n,MaxC,K,Godel1,Godel2,co):
#inputs:
#(i) a pos. integer L, 
#(ii) a symbol t
#(iii, iv, v) symbols H,P,M (for ordinary generating functions)
#(vi) symbol n
#(vii) a positive integer MaxC (for the maximal complexity of the recurrence)
#(viii) a large positive integer (for the asympotics and for displaying
#a specific value)
#(ix,x) positive integers Godel1 and Godel2 (Godel1<Godel2)
#(xi): a pos. integer co, for labelling the Theorem
#Outputs: a theorem about the sequences enumerating
#half-pyramids, pyramids, and towers consisting of L-mers
#where it is not allowed to place a brick exactly alligned on it
#(i.e. it must jut-out either to the left or to the right)
#For example, try:
#Mishpat3(2,t,H,P,M,n,10,1000,30,100,1):
Mishpat3:=proc(L,t,H,P,M,n,MaxC,K,Godel1,Godel2,co) local guH,guP,guM,N,
luH,luP,luM,ope,C,i:

guH:=KamaOneBHnzFe(L,Godel2):
guP:=KamaOneBPnzFe(L,Godel2):
guM:=KamaOneBMnzFe(L,Godel2):


print(`Theorem Number`, co):

print(`This theorem is about the enumeration of half-pyramids,`):
print(`pyramids, and towers using L-mers of the same size`, L):
print(`Where it is allowed to place a brick anywhere on top of`):
print(`a brick below it, EXCEPT exactly lined-up`):

print(`Part 1: Half-Pyramids`):

print(`Let  h(n) be the number of ways of building`):
print(`half-pyramids with n bricks, all of size`, L):
print(`where all the interfaces are allowed.`):
print(`EXCEPT exactly lined-up`):

print(`The values of h(n) for n from 1 to`, Godel1, `are `):
print([op(1..Godel1,guH)]):


luH:=Meda(guH,t,H,n,N,MaxC,K):

if luH<>FAIL and nops(luH)>=1 then
print(`Let H(t) be the ordinary generating function of h(n), i.e.`):
print(H(t)=1+Sum(h(n)*t^n,n=1..infinity)):
print(`then, using the clever Betrema-Penaud argument, as extended`):
print(`by Doron Zeilberger , it is readily seen that H(t) satisfies`):
print(`the algebaric equation`):
print(luH[1]=0):
print(`and in Maple input format`):
lprint(luH[1]=0):
fi:

if  nops(luH)>=3 and luH[2]<>FAIL then
ope:=luH[2]:
print(`The sequence h(n) satisfies the following linear recurrence`):
print(`equation with polynomial coefficients`):
print(add(coeff(ope,N,i)*h(n+i),i=0..degree(ope,N))=0):
print(``):
print(`just for kicks, here it the value of `, h(K)):
print(luH[3]):
print(``):
fi:



if nops(luH)>=4 and luH[4]<>FAIL then
print(`h(n) is, apprx., asymototically equal to`):
print(C*n^luH[4][1]*luH[4][2]^n):
print(`for some constant C`):
fi:

if nops(luH)>=5 and luH[5]<>FAIL then
print(`that is approximately`):
print(luH[5]):
fi:


print(`Part 2: Pyramids`):

print(`Let  p(n) be the number of ways of building`):
print(`pyramids with n bricks, all of size`, L):
print(`where all the interfaces are allowed.`):
print(`EXCEPT exactly lined-up`):

print(`The values of p(n) for n from 1 to`, Godel1, `are `):
print([op(1..Godel1,guP)]):


luP:=Meda(guP,t,P,n,N,MaxC,K):


if luP<>FAIL and nops(luP)>=1 then
print(`Let P(t) be the ordinary generating function of p(n), i.e.`):
print(P(t)=1+Sum(p(n)*t^n,n=1..infinity)):
print(`then, using the clever Betrema-Penaud argument, as extended`):
print(`by Doron Zeilberger , it is readily seen that P(t) satisfies`):
print(`the algebaric equation`):
print(luP[1]=0):
print(`and in Maple input format`):
lprint(luP[1]=0):
fi:

if  nops(luP)>=3 and luP[2]<>FAIL then
ope:=luP[2]:
print(`The sequence p(n) satisfies the following linear recurrence`):
print(`equation with polynomial coefficients`):
print(add(coeff(ope,N,i)*p(n+i),i=0..degree(ope,N))=0):
print(``):
print(`just for kicks, here it the value of `, p(K)):
print(luP[3]):
print(``):
fi:



if nops(luP)>=4 and luP[4]<>FAIL then
print(`p(n) is, apprx., asymototically equal to`):
print(C*n^luP[4][1]*luP[4][2]^n):
print(`for some constant C`):
fi:

if nops(luP)>=5 and luP[5]<>FAIL then
print(`that is approximately`):
print(luP[5]):
fi:




print(`Part 3: Towers`):

print(`Let  m(n) be the number of ways of building`):
print(`TOWERS with n bricks, all of size`, L):
print(`where all the interfaces are allowed.`):
print(`EXCEPT exactly lined-up`):

print(`The values of m(n) for n from 1 to`, Godel1, `are `):
print([op(1..Godel1,guM)]):



luM:=Meda(guM,t,M,n,N,MaxC,K):



if luM<>FAIL and nops(luM)>=1 then
print(`Let M(t) be the ordinary generating function of m(n), i.e.`):
print(M(t)=1+Sum(m(n)*t^n,n=1..infinity)):
print(`then, using the clever Betrema-Penaud argument, as extended`):
print(`by Doron Zeilberger , it is readily seen that M(t) satisfies`):
print(`the algebaric equation`):
print(luM[1]=0):
print(`and in Maple input format`):
lprint(luM[1]=0):
fi:

if  nops(luM)>=3 and luM[2]<>FAIL then
ope:=luM[2]:
print(`The sequence m(n) satisfies the following linear recurrence`):
print(`equation with polynomial coefficients`):
print(add(coeff(ope,N,i)*m(n+i),i=0..degree(ope,N))=0):
print(``):
print(`just for kicks, here it the value of `, m(K)):
print(luM[3]):
print(``):
fi:




if nops(luM)>=4 and luM[4]<>FAIL then
print(`m(n) is, apprx., asymototically equal to`):
print(C*n^luM[4][1]*luM[4][2]^n):
print(`for some constant C`):
fi:

if nops(luM)>=5 and luM[5]<>FAIL then
print(`that is approximately`):
print(luM[5]):
fi:

print(`this ends this exciting theorem, whose proof is routine.`):

end:




#Mamar2(L1,t,H,P,M,n,MaxC,K,Godel1,Godel2):
#outputs a paper about Mishpat2(L,t,H,P,M,n,MaxC,K,Godel1,Godel2):
#for L from 2  to L1, try:
#Mamar2(3,t,H,P,M,n,15,1000,30,100):
Mamar2:=proc(L1,t,H,P,M,n,MaxC,K,Godel1,Godel2) local L, co,t0:
t0:=time():
co:=1:

print(`Enumeration of Half-Pyramids, Pyramids, and, most interestingly,`):
print(`Towers, With One Allowed Brick for Brick-Sizes from 2 to`, L1):
print(`Where all interfaces are allowed. `):
print(``):
print(`By Shalosh B. Ekhad`):

for L from 2 to L1 do
 Mishpat2(L,t,H,P,M,n,MaxC,K,Godel1,Godel2,co):
 co:=co+1:
od:

print(`This ends this exciting article. It took`, time()-t0, `to create it.`):

end:


#Mamar3(L1,t,H,P,M,n,MaxC,K,Godel1,Godel2):
#outputs a paper about Mishpat3(L,t,H,P,M,n,MaxC,K,Godel1,Godel2):
#for L from 2  to L1, try:
#Mamar3(3,t,H,P,M,n,15,1000,30,100):
Mamar3:=proc(L1,t,H,P,M,n,MaxC,K,Godel1,Godel2) local L, co,t0:
t0:=time():
co:=1:

print(`Enumeration of Half-Pyramids, Pyramids, and, most interestingly,`):
print(`Towers, With One Allowed Brick for Brick-Sizes from 2 to`, L1):
print(`Where all interfaces are allowed, EXCEPT EXACTLY LINED-UP `):
print(``):
print(`By Shalosh B. Ekhad`):

for L from 2 to L1 do
 Mishpat3(L,t,H,P,M,n,MaxC,K,Godel1,Godel2,co):
 co:=co+1:
od:

print(`This ends this exciting article. It took`, time()-t0, `to create it.`):

end:



#Mamar1(L1,t,H,P,M,n,MaxC,K,Godel1,Godel2):
#outputs a paper about Mishpat1(L,t,H,P,M,n,MaxC,K,Godel1,Godel2):
#for all subsets L of {1, ..., L1} with at least two members. Try:
#Mamar1(3,t,H,P,M,n,15,1000,30,100):
Mamar1:=proc(L1,t,H,P,M,n,MaxC,K,Godel1,Godel2) local k, co,t0,k1,L:
t0:=time():
co:=1:

print(`Enumeration of Half-Pyramids, Pyramids, and, most interestingly,`):
print(`Towers, With Several Allowed Brick for All Sets of Brick Sizes`):
print(`with at least two members, up to brick-size`, L1):
print(``):
print(`By Shalosh B. Ekhad`):

for k from 2 to L1 do
L:=choose(L1,k):
for k1 from 1 to nops(L) do
 Mishpat1(L[k1],t,H,P,M,n,MaxC,K,Godel1,Godel2,co):
 co:=co+1:
od:
od:

print(`This ends this exciting article. It took`, time()-t0, `to create it.`):

end: