#####################################################################
##MDD.txt: Save this file as  MDD.txt                                #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read MDD.txt <Enter>                               #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#zeilberg at math dot rutgers dot edu                                #
######################################################################
 
#Created: Nov. 2016
 
print(`Created:  2016`):
print(` This is MDD.txt `):
print(`It is one of the packages that accompany the article `):
print(`The Number of Monomer-Dimer Tilings of an m by n Rectangle with a Fixed Proportion of Monomers `):
print(`by Gleb Pogudin  and Doron Zeilberger`):
print(`It also needs the Maple package TILINGS available from Zeilberger's website`):
print(`and also available from Pogudin's and Zeilberger's websites `):
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(`---------------------------------------`):
 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(`---------------------------------------`):

with(combinat):


Ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: Empir `):
 print(``):

else
ezra(args):
fi:

end:

Ezra:=proc()

if args=NULL then
 print(`The main procedures are: CC, fmkNdata, fmkNodata,fmkE, fmkN, Gleb, GFmd, GFmdPC, MDalg, PCdata, SeqT , SeqTo, Story`):
 print(` `):

elif nops([args])=1 and op(1,[args])=CC then
print(`CC(m,k,N): The "connective constant" for the sequence SeqT(m,k). Try:`):
print(`CC(2,1/2,40); `):


elif nops([args])=1 and op(1,[args])=Empir then
print(`Empir(L,x,P): inputs a sequence of numbers L, starting at n=0 and guesses an algebratic equation for its ogf P(x)`):


elif nops([args])=1 and op(1,[args])=fmkNdata then
print(`fmkNdata(m,A,N): a table of [k,fmkN(m,k,N)] for all k=a/b with gcd(a,b)=1, 1<=a<b<=A. For even m. Try:`):
print(`fmkNdata(2,3,100);`):

elif nops([args])=1 and op(1,[args])=fmkNodata then
print(`fmkNodata(m,A,N): a table of [k,fmkN(m,k,N)] for all k=a/b with gcd(a,b)=1, 1<=a<b<=A. For odd m. Try:`):
print(`fmkNodata(3,3,100);`):

elif nops([args])=1 and op(1,[args])=fmkE then
print(`fmkE(m,k,N): Let A_m(k,n) be the number of monomer-dimer tilings of an m by b*n rectangle (m even)`):
print(`with monomer density k, finds (EXACTLY) the limit of log(A_m(k,n))/(n*m*b) as n goes to infinity`):
print(`Try:`):
print(`fmkNE(2,1/2,100);`):

elif nops([args])=1 and op(1,[args])=fmkN then
print(`fmkN(m,k,N): Let A_m(k,n) be the number of monomer-dimer tilings of an m by b*n rectangle (m even)`):
print(`with monomer density k, finds (NUMERICALLY) the limit of log(A_m(k,n))/(n*m*b) as n goes to infinity`):
print(`Try:`):
print(`fmkN(2,1/2,100);`):

elif nops([args])=1 and op(1,[args])=fmkNo then
print(`fmkNo(m,k,N): Let A_m(k,n) be the number of monomer-dimer tilings of an m by b*n rectangle (m odd)`):
print(`with monomer density k, finds (NUMERICALLY) the limit of log(A_m(k,n))/(n*m*b) as n goes to infinity`):
print(`Try:`):
print(`fmkNo(2,1/2,100);`):

elif nops([args])=1 and op(1,[args])=GFmd then
print(`GFmd(m,t,z,w): inputs a positive integer m, and variables t,z,w, outputs the rational function of t,z,w`):
print(`whose coefficient of t^n*z^i*w^j is the number of tilings of an m by n rectangle with exactly i dimers and j monomers.`):
print(`Try:`):
print(`[seq(GFmd(i,t,z,w),i=1..3)];`):

elif nops([args])=1 and op(1,[args])=GFmdPC then
print(`GFmdPC(m,t,z,w): a pre-computed version of GFmd(m,t,z,w) (q.v.), for 1<=m<=6, otherwise does GFmd(m,t,z,w).`):
print(`Try:`):
print(`[seq(GFmdPC(i,t,z,w),i=1..3)];`):

elif nops([args])=1 and op(1,[args])=Gleb then
print(`Gleb(p): the entroy function for the plane up to 43 terms`):

elif nops([args])=1 and op(1,[args])=MDalg then
print(`MDalg(m,k,N,x,P): inputs an even integer m, a rational number k (=a/b,say), a large pos. integer N, and symbols x,P`):
print(`outputs a guessed (but rigorizable) algebraic equation for the ordinary generating function of the sequence`):
print(`a(n):=Number of monomer-dimer tilings of an m by b*n rectangle with m*b*n*k= m*n*a monomers and the rest m*n*(b-a)/2 dimers`):
print(`Try:`):
print(`MDalg(2,1/2,40,x,P);`):

elif nops([args])=1 and op(1,[args])=PCdata then
print(`PCdata(): pre-computed data for curves using fmkN(m,k,N). It returns the curves for m=2,m=3,m=4,m=5,m=6, and m=infinity, respectively`):
print(`Try: `):
print(` PCdata(); `):
print(` Try: `):
print(` plot(PCdata()); `):

elif nops([args])=1 and op(1,[args])=SeqT then
print(`SeqT(m,k,N): inputs a positive even integer m, a rational number k=a/b, and outputs the sequence of the number of tilings of a m by (n*b) rectangle`):
print(`with a/b*(m*n*b)=m*n*a monomers and (m*n*b-m*n*a)/2=m/2*n*(b-a) dimers for n from 1 to N. Try:`):
print(`SeqT(2,1/2,40);`):

elif nops([args])=1 and op(1,[args])=SeqTo then
print(`SeqTo(m,k,N): inputs a positive odd integer m, a rational number k=a/b, and outputs the sequence of the number of tilings of a m by (2*n*b) rectangle`):
print(`with a/b*(2*m*n*b)=2*m*n*a monomers and (m*2*n*b-m*2*n*a)/2=m*n*(b-a) dimers for n from 1 to N. Try:`):
print(`SeqTo(3,1/2,40);`):

elif nops([args])=1 and op(1,[args])=Story then
print(`Story(m,k,N,x,P): tells the story of MDalg(m,k,N,x,P). Try:`):
print(`Story(2,1/2,50,x,P):`):

print(``):


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



###start data for curves
#PCdata(): pre-computed data for curves using fmkN(m,k,N). It returns the curves for m=2,m=3,m=4,m=5,m=6, and m=infinity, respectively
#Try:
#PCdata();
#Try:
#plot(PCdata());

PCdata:=proc() local L2,L3,L4,L5,L6:
L2:=
[[0, .2406059125298017237488794567121842115675921671928302598305090844200819338\
041108872060047145613617376], [1/10, .42650368658969872599672070312651041728923\
51339886941714806917593596926746434035329247958399361918063], [1/9, .4399363030\
7924764563508508838900991656863579697785238837194632735040025019750730999280540\
94339743753], [1/8, .4555743404545550560118341530775027606997118109572899298009\
970934162074502917938064105154650944522290], [1/7, .473948562716645529691692499\
9565587576240164297530874640967327903280634300493460177967052294623050209], [1/
6, .495685444217973992557276838530346102525426885913557797187860078892852819616\
7369825433470614290573167], [1/5, .52134594685357671652206911030894959041378112\
40899880808243912262045684941656795844559578634616781264], [2/9, .5356447022167\
7885062036992284440987866323915748072114420401811007299787287635286452227092574\
36822053], [1/4, .5506156462572977438590883854722104780476244145160566392702280\
592678699883975326891580782317811576260], [2/7, .565480079144925085160659229759\
9185018796866745180358282048811007064182395328851748855317324833221551], [3/10,
.570126099385350161833456596912715079218081982472351383192627816312494906259208\
7647374505143807019615], [1/3, .57823605455057866746386439348314670082938420700\
87450478961299951982409426471101545541161575613863140], [3/8, .5832897725813679\
6548506938326300675886393794979898221172843845633156470043732023039450180914725\
41536], [2/5, .5837424237803130346139265777460058374789016442150311044109662019\
630087615904225851076939958257601860], [3/7, .581991731568947393104550210176441\
2197854066696120359748108738237558614816432879846852327876258214986], [4/9, .57\
9999573633936950700658639722841721380357691865387976940214303043572175567169152\
5090987262160166928], [1/2, .56745385754479755848181177919006432491741143932793\
30879284878067172575363721346302858308208459180160], [5/9, .5464842282633104931\
6297035725506673742859211465992509250040975532877342289937039006137040119505910\
69], [4/7, .5389663364868319006105600515431498680067604808484853264803991187004\
251528801085625372551666142574747], [3/5, .523731350057588077745872567080984798\
9449143139766062055073554984633481184366625837478457266395576818], [5/8, .50859\
9819045100191216213420587888632006637780993502497487667351385192958595293128742\
2633997095894664], [2/3, .47959034537435507313547436971326834890814021775705151\
88499071738329463704609621076304586736807007890], [7/10, .452899605116386432856\
1359196163252966603583373803785731173005993432219491442616458909780399814316902
], [5/7, .440480894454208380265240145012855774073381593203747826224500914943989\
6732108945192047393160333671600], [3/4, .40676554128069869711308446141973778017\
17760504084820015519806896627740035383220367047929728550469038], [7/9, .3777833\
9812311090326739175535811679609739657803189233419495807707316523479902335084106\
26165946539876], [4/5, .3527403861940681820922239086160021573380497565066695738\
924783856292027278306506830525380108631621148], [5/6, .311828412649164380022457\
2636055345790916723951742129151836057930718231098758235072385387782372298620],
[6/7, .279896123574213749771239966826472116847919281643480435573854799187869672\
2868386965190474708764509963], [7/8, .25428908687848307492211326103516397570797\
99016236248439218508039778199970929714562558705365857286005], [8/9, .2332839546\
3472837131503006016761823220042349367992994204447282776097319845619617120981939\
77851821142], [9/10, .215726518498798241072562545681290480476148402626396233974\
2713022853769905017254386016421968848995240]]:




L3:=
[[1/10, .4476004424774455177421813638\
686560992832513835641504679802480925851540224866671688543392974512504570], [1/9
, .4622497868556369906832211583400204469781960952638033855889958532355704587975\
234257120450776080135760], [1/8, .479128892544286041594676655475127577517111464\
2807944320173956036257330930413165397045325994591060812], [1/7, .49873122022443\
8378702854495228038810954082544126365022142514152259078747393497682361673001322\
1575386], [1/6, .52160998837032422850421284194147085588537658952266790975631742\
07882221724280664719878139772881110610], [1/5, .5481794529333348527592234221696\
332675555452363960966312217710051531213437539149128492144092214974325], [2/9, .\
5627583821801587211800611482907953141488014172489029511067053613660974351318919\
236157644535874167970], [1/4, .577813473551306079703471588361659507084290696096\
0193989808056961373278754627379384912703028098277000], [2/7, .59245502061701061\
6078504380161119356980895401874125209176990520542527743045291337875193359870751\
4605], [3/10, .5969310230726701659659049397328265667183643533088857304807223416\
915084962494577637017896174210798250], [1/3, .604506018628375263016332955877333\
6138110023292038143116507503168623219902154419701858240736517701620], [3/8, .60\
8666445526841819262119454093500797834154065578346909354839732868763493144032573\
8950408809617885300], [2/5, .60849256027727895465534670079792168043075130392776\
54406831122997888743405283462205638027770898285755], [3/7, .6059557312722338779\
1044683893267882697427785681213529788393634141567570438589819903305488205696929\
70], [4/9, .6034991002217855068539445245337565826794035254302307608458431862415\
495714781171371203234154638530190], [1/2, .589200212973265290423586848012104192\
3933597574586813032813057089301687272798638570559412816547023370], [5/9, .56629\
9438834095525985079686553417714784723497874176744542957279790180979733412442249\
4791102449092660], [4/7, .55820390914227050542560206431744452652590127053624477\
71184667125705773410648156963982008912482241460], [3/5, .5419014735123281910361\
569771254852868011330245145239287040537886340651156681007873846722415501373460]
, [5/8, .5258070538286288124294946293512042225749792263988656413676399371631707\
867599077228597727503736725035], [2/3, .495140884432851979021656767090225060733\
2179442321513944901898483932333399209259034045279224625746976], [7/10, .4670741\
2888538011074036844634783621383151151233885297059005515875089620734285723598464\
92281861971402], [5/7, .4540543672290432036847080954566980838043972585890893940\
917598690800523354917160505399305575593625334], [3/4, .418803966339132991830855\
2168164363068971912820955723863106555214071218126055762016514509438849104250],
[7/9, .388594276710182687605878959593262624444960776504122604022031343784008064\
5584090313685878629836020163], [4/5, .36255101676017805483351230510563848870676\
52826344501107429701279120028538557103272413227255265958447], [5/6, .3201041284\
9673871739976507768313651281115981604707203058138736340331610156508758798253990\
69094287939], [6/7, .2870510925462744280218177095082498845862665958732331645247\
323465364646226135577536823425786897182950], [7/8, .260590127024239120279406349\
9567128506439772305209098951710937121674703949600211077395817404991733996], [8/
9, .238912885069502485364839241060930492778992197559051126323572305120149736474\
7368076168389481256424515], [9/10, .2208127461237711077771935686598529173310460\
780681925018484877567965782893348324959141595948780033134]]:

L4:=
[[0, .2609982087725492799745257988388437406076138228067772291479345229153926571\
718380835456288491951577142], [1/10, .46662152230657006946327786185453003858620\
88371445635773603234534213876157759551738806779427788555678], [1/9, .4809449173\
0172713627194576459935982966927458901450230059530790380040860999856415536804570\
43779577136], [1/8, .4974752192008457159605584361655703317319817906572606239233\
735888129747549439148296020280721268504744], [1/7, .516689880996980273077791972\
6338299940628398353640448970344594946019474680375543249062408560685068926], [1/
6, .539110813969269548785686203960956422986995846435679104442120534719011440326\
6725526389099593790214970], [1/5, .56508812072496493714414570523453307597293158\
61155300495556719034815961345928543843245546464151989600], [2/9, .5792841777737\
2987214554958060351901715580312329014988766355411372962898315288044333469017883\
73075102], [1/4, .5938537347695697421023857012075258272091128926179395226555326\
411187477411901668758496950489175066945], [2/7, .607866126031015568735758832304\
9879537462605636692455326323090306319731628855194630400151933347729629], [3/10,
.612085890581359888488632713640663311949153277496019641764837371760047063689450\
6920821801150473240815], [1/3, .61904480013225881992187382339132559979612972058\
58484145273948627211293402010119957006939795877969803], [3/8, .6224246037838868\
5865486716010515465984580550175873188834148389801919566722720363457987307783304\
30302], [2/5, .6217671644314974928640830448523669964210067846977583712989120758\
691741300139439239033048138404664696], [3/7, .618671474175862166858984649587386\
9964056421799882141197874866121300818081803787819191320149425886057], [4/9, .61\
5900274592779747548874329015750387946111400344548204987177146027017940180875962\
7233470534693663958], [1/2, .60046813217347261143390021566876289753952208274832\
15226179749725246937943681231597107483428693040575], [5/9, .5764202118314540762\
8751288422014789428405957329171146560425152996902985851609480952834545729087393\
13], [4/7, .5679886155998489335755783267923591498775535539986000004490378176324\
133157187481962176228179644729700], [3/5, .551075725741410868831316155666218541\
0576033669963520089482691749627518967514323816886617494705037860], [5/8, .53444\
4230618870743901651241832861630082614217946401907311000554107277289188166466810\
9496938808427928], [2/3, .50286448268769985128385420461463178612349530761636935\
08696432271078262862891310390114809996037083357], [7/10, .474065816348599940575\
1108127622408520271128946617321372119260470733909188474875494042478273367916768
], [5/7, .460726714659093883332958526751227585625054915025984494231230272857916\
3017939922453080299735262678136], [3/4, .42467139875688569712341869366025090021\
69644495509797790346322041258558410722898745997572532998462085], [7/9, .3938341\
0058652424738165993282479471614593517279378075289889533616702203612057897566641\
94754121365739], [4/5, .3672804549052901819323440068528623840560345322430134871\
527137247615322646112440561549956936324368918], [5/6, .324065061149796888000353\
3030590876292180868783758774712548685443304709971832088076430333368592283892],
[6/7, .290458305086027478781102816701153212042292772861764049724781981955971755\
3449093341307132845351143717], [7/8, .26357943060540359521190290972239633410720\
39683513301418967418402334039269943352761729157259273100734], [8/9, .2415755921\
0978104232119149102553520416516611675138055217629236045044842128985320420969202\
64260036752], [9/10, .223213191983881112256895403705670606501380822416514713968\
4190633823730328852706077987085694645816633]]:




L5:=
[[0, .2529222887091617973659574689429223542927155178542752355333047338343710952\
770733662987330235130462685], [1/10, .47625296317370050023190969092534341062219\
18931475233297018683965055543339605062790764567563937435055], [1/9, .4906847515\
0446907847278142250707024516463349990241680739911015150648972948400013236592122\
99242712749], [1/8, .5072913750751801267173382288118972766435330757054280089152\
857332127095746042457019951168513214964860], [1/7, .526535659973609733220696391\
7307329592678236234355531628095543169597590080303375776731797181881785210], [1/
6, .548915510701869451932047066705847332238513467473009764376577805194181462733\
3531349695092010730757860], [1/5, .57473779700391979544876017716289449540802640\
57389079635419461629604365287787357854538733524790847235], [2/9, .5887882097591\
5841539340318610058777764234155669067183219107051728872708392474867625539674526\
40332350], [1/4, .6031462780547753709198718051648586565111814770156947697445880\
311886622945391850602823575231704460350], [2/7, .616846192322129400851607452364\
0218306680283627128588166189962398047989388948379880721941117378935155], [3/10,
.620931694956083657323810776968664661065557839764364014284786417226114439268695\
4058010756468250369340], [1/3, .62756410331785226102956945258839603662600151471\
65555545186418119946036794909254723345467338729947705], [3/8, .6305047261504800\
9988424828312215785926786870896720118377204714673217877775732638399083775270303\
38240], [2/5, .6295742233874080086560795297249091526589975067393068149112724964\
525125674314555836663351446488762940], [3/7, .626156202673373593240874076927957\
2864338646433525400895977494693316215072866800804675943017605092675], [4/9, .62\
3202250846655564615692578180485446190272929871063351932033528174050933134243593\
4760168201816106080], [1/2, .60711718309196926586577831545410539051239771986812\
63927278221582776341912326660485185750110876217455], [5/9, .5823850168328119312\
0453790635290437308426640105875536589552396643520435277541975334537453071689642\
10], [4/7, .5737554646736524203280594104986570303142526474733815070162679855336\
231380742218885235612610639429485], [3/5, .556482413800514220545908158796049099\
4603754759328524659822062710271805922836036707005155759968614890], [5/8, .53953\
0808544561882632461037046535118603143915313268210147501091894223374245632238058\
6454252253138690], [2/3, .50741355330526549304370938116916249791896437519287403\
02101436151209847455318929265317945487042281645], [7/10, .478175226912208433516\
4603547766126400992121295373323984361568270432338872481274172323369774563491361
], [5/7, .464647586157420058962878482889971976171849633032176253744046412651080\
2578837365219731909291201044368], [3/4, .42811733577289220333816202357563981296\
16687653558634201433480302110246406173828272259095900608621150], [7/9, .3969052\
3792336105909402123204469943927100922170372164226459583096967728216564917888902\
35861382694372], [4/5, .3700517228010491700297122118327966936150161338767968399\
169435877791728745644025257878649836550200576], [5/6, .326381955876331308394520\
1822607549324561958314069244247278180113625216938463204987902846416763990234],
[6/7, .292448596475395863066475534711435190604108393662185079966323742248855034\
0010698040520668849415662596], [7/8, .26532371638502818363188195423988273867183\
00020264277479191322529405094798366686023340733329842572962], [8/9, .2431279439\
1899472153193587236429597961739764205538053730212398396564424035929067596137156\
54165716104], [9/10, .224611633633503104454956693714594368034713582087280416253\
8673127335153419507405685011239705142021506]]:



L6:=
[[0, .2698623053679531284318890227584687741865021994212489579854461929873566288\
330244237745850426475123728], [1/10, .48319870898025216848293467121724050010658\
85897058581262170436238968902835285393024066238717996541555], [1/9, .4975922146\
7082195165652121594287471784468728194949116509353635427088332566349701145802965\
43693381280], [1/8, .5141497957958156209186581019428241155796500497728398814060\
316278583987310316868370194237487367215952], [1/7, .533325835061081189892867811\
5938648365817559024259900908582811563253264661331423189181381664826030661], [1/
6, .555603264632997286845281890402206605568261789008223072629721924363372597380\
7390342227006692989699120], [1/5, .58126058838421800863903410867219961841964329\
94914518477739547208744485712993997505670796785077696796], [2/9, .5951888046712\
0898350683312591646067880833583630140080979348701248192504903609709983865809428\
81850939], [1/4, .6093822149713012060046616534728505110546621163671865242744703\
331738673966508336410628929773514208912], [2/7, .622854672276189168786015692862\
9768809083457477440289444246846713740947211258122377442780695069970293], [3/10,
.626844884132152739939088215309490157173814348346346919466600453600929976156342\
6598995200469910338103], [1/3, .63324647159474683432324010872966056011442965343\
77629096360806288291968301486997352524814426333438783], [3/8, .6358851528384504\
0119823223827210643396711120927342810629892683212402414396792182594314233794525\
75272], [2/5, .6347670233169747565106933063780331570960901801329876130130948011\
024115681938495022861989438462352360], [3/7, .631129676823314723942058670328443\
3059698806613444597876771793943255888863533167620470457264772891000], [4/9, .62\
8051775176884079255669792440782089434837749535909869169262708150073515553543771\
6261469554259008339], [1/2, .61152205839447821163943037796675872960447275911880\
59242313498204468162305144569516600887389213418535], [5/9, .5863319685183704672\
0771882454917119220363244616540429421944775156747850066332126870196594591344333\
97], [4/7, .5775692601217849198437561817751000751299594926377333136799914848545\
254187170412349961184641631526731], [3/5, .560054326164583388078905399079997875\
6481226529334557806073612712166806331332693149875338884429403246], [5/8, .54288\
9030114227117017774352411269236189850101351122293459314282630622037639413739235\
2639976484769591], [2/3, .51041132628401870405526125445207438357892867786812847\
96000611152555957385524006793300985047017217543], [7/10, .480881749633614993599\
4614471404612631840197274027619875251655964788565684148268712483879716736159898
], [5/7, .467228373245920263819085473033932092105130015459839068503267711886364\
7262818901271095935676349910757], [3/4, .43038190256586459283808127089527441595\
82593991729799209718504677322339155670919119185446781440766218], [7/9, .3989223\
5378264456379057390512841934314661844749754879049335016893906546340057731656005\
61056838703152], [4/5, .3718697511699974668619672522624081364941832052930861958\
604310565263633651586099024883655073125523750], [5/6, .327900073399094768456671\
2279178515401186243483922971920928169140713140086003585689929938486152628567],
[6/7, .293751557539277625529166928950256436280543173896823933316489805673737159\
6505779255503109765386595810], [7/8, .26646485348854792599481929929452627151811\
89562938497209444988732355215667754923067094385890550785956], [8/9, .2441429700\
7446467320530471822489950197185721328324561642573350776801690066987723917081203\
35674677320], [9/10, .225525625501572697140359770599465397885057954332583830951\
1111360759572014152896539430000789975254852]]:

Linf:=evalf([seq([L2[i][1],Gleb(1-L2[i][1])],i=2..nops(L2))]):

[L2,L3,L4,L5,L6,Linf]:
end:
###end data for curves

#########FROM TILINGS
######################################################################
## TILINGS: Save this file as TILINGS to use it,                     #
# stay in the                                                        #
## same directory, get into Maple (by typing: maple <Enter> )        #
#then type: read TILINGS: <Enter>                                    #
## Then follow the instructions given there                          #
##                                                                   #
## Written by Doron Zeilberger, Rutgers University ,                 #
##  zeilberg@math.rutgers.edu.                                       # 
######################################################################

with(combinat):

Digits:=100:

print(`First Written: Jan. 2006: tested for Maple 10  `):
print(`Version of Jan. , 2006:  `):
print():
print(`This is TILINGS, A Maple package`):
print(`accompanying Shalsoh B. Ekhad and Doron Zeilberger's  article: `):
print(`"Automatic CounTilings" `):
print(`-------------------------------`):
print(`Version of Sept. 30, 2015, procedure Dimer added`):
print(`-------------------------------`):
print():
print(`The most current version is available on WWW at:`):
print(` http://www.math.rutgers.edu/~zeilberg/tokhniot/TILINGS .`):
print(`Please report all bugs to: zeilberg at math dot rutgers dot edu .`):
print():
print(`For general help, and a list of the MAIN functions,`):
print(` type "ezra();". For specific help type "ezra(procedure_name);" `):
print(`For a list of the supporting functions type: ezra1();`):
print():
 

ezra1:=proc()
if args=NULL then

print(` The supporting procedures are: BuildTree1, Kmn`):
print(` `):
print(`For help with a specific procedure, type "ezra(procedure_name);"`):
else
 ezra(args):
fi:
end:

ezra:=proc()
if args=NULL then

print(` TILINGS: A Maple package for counting tilings `):
print(`of rectangles by any set of tiles`):
print(`The MAIN procedures are: AaveT, AllTilings, AllTilingsNice, Astat, AstatV `):
print(`BuildTree, BuildTreeH, Dimer, GFt, GFtH , GFtMD, ProveKasteleyn,  ProveKasteleyn1, RandTiling,RandTilingNice,  Rt`):
print(` Sidra, SidraW, Sipur, SipurArokh `):

elif nargs=1 and args[1]=AllTilings then
print(`AllTilings(Tiles,k,n0): The set of tilings of the k by n0 board `):
print(`using the tiles of Tiles, in set-partition form`):
print(`For example, try: AllTilings([Rt(1,2),Rt(2,1)],2,2);`):


elif nargs=1 and args[1]=AllTilingsNice then
print(`AllTilingsNice(Tiles,k,n0): The set of tilings of the k by n0 board `):
print(`using the tiles of Tiles, in matrix form.`):
print(`The numbers do not have any significance, but cells labelled with the`):
print(`same number belong to the same tile.`):
print(`The tilings  are listed one by one.`):
print(`For example, try: AllTilingsNice([Rt(1,2),Rt(2,1)],2,2);`):

elif nargs=1 and args[1]=AaveT then
print(`AaveT(k,ListOfTiles): Given a postive integer k and a symbol n,`):
print(` and list of tiles ListOfTiles , `):
print(` returns the list asymptotic averages (divided by n) of the number of tiles`):
print(`of tiling a k by n rectangular board by the  tiles in`):
print(`followed by the average total number of total number of tiles (divided by n)`):
print(`For example, try: AaveT(2,[Rt(1,2),Rt(2,1)]) `):


elif nargs=1 and args[1]=Astat then
print(`Astat(k,ListOfTiles,n): Given a postive integer k and a symbol n`):
print(` and given a List of tiles`):
print(`ListOfTiles`):
print(`returns: (1) The lists of for the`):
print(`asymptotics (in n) for pairs [average,variance] for the`):
print(`r.v. number of tiles of type ListOfTiles[i], in a random tiling by the rectangular tiles`):
print(`of Tiles (in order), followed by the pair [av,variance] for the`):
print(`Total number of tiles, followed by the asympto. correlation`):
print(`matrix for example try:`):
print(`Astat(4,[Rt(1,2),Rt(2,1)],n):`):

elif nargs=1 and args[1]=AstatV then
print(`AstatV(k,ListOfTiles,n): Verbose version of Astat (q.v.) `):
print(`For example, try:`):
print(`AstatV(4,[Rt(1,2),Rt(2,1)],n):`):

elif nargs=1 and args[1]=AstatVs then
print(`AstatVs(K,ListOfTiles,n): AstatV(k,ListOfTiles,n) for k from 1 to K `):
print(`For example, try:`):
print(`AstatVs(4,[[1,2],[2,1]],n):`):

elif nargs=1 and args[1]=Asy1 then
print(`Asy1(f,t,n): the Asymptotics of the coeff. of the rational function f`):
print(`try for example: Asy1(1/(1-t-t^2),t,n);`):

elif nargs=1 and args[1]=BuildTree then
print(`BuildTree(k,Tiles,t): The tree together with the generating functions in t only.`):
print(`For example, try: BuildTree(2,[Rt(1,2),Rt(2,1)],t);`):

elif nargs=1 and args[1]=BuildTreeH then
print(`BuildTreeH(k,Tiles,t,H): The tree together with the generating functions in H.`):
print(`For example, try: BuildTreeH(2,[Rt(1,2),Rt(2,1)],t,H);`):

elif nargs=1 and args[1]=BuildTree1 then
print(`BuildTree1(k,Tiles): Builds a tree without the generating functions`):
print(`For example, try: BuildTree1(2,[Rt(1,2),Rt(2,1)];`):


elif nargs=1 and args[1]=Dimer then
print(`Dimer(n,z): the C-finite representation of the generating function for the dimer problem`):
print(`in an n by m strip (m general, n specific) with vertical variable z and horiz. variable 1, try:`):
print(`Dimer(4,z);`):

elif nargs=1 and args[1]=GFt then
print(`GFt(k,Tiles,t): The generating function, in the variable t, whose coeffs. of`):
print(`t^n is the number of tilings of the rectangular k by n board `):
print(`(k-side parallel to the y-axis, and n-side parallel to the x-axis), with the tiles of the list Tiles.`):
print(`The tiles of Tiles are described as normalized sets of lattice points`):
print(`Or in case of rectangular tiles, one can use Rt(a,b) (q.v.) for tiles `):
print(`For example, try: GFt(2,[Rt(1,2),Rt(2,1)],t);`):

elif nargs=1 and args[1]=GFtH then
print(`GFtH(k,Tiles,t,H): The generating function, in the variable t,`):
print(`whose coeffs. of`):
print(`t^n is the weight-enumerator of tilings of the rectangular k by n board with`):
print(`(k-side parallel to the y-axis, n-side parallel to the x-axis) with the tiles of the list Tiles.`):
print(`The tiles are described as normalized sets of lattice points.`):
print(`(For rectangular tiles, use Rt(a,b) for a b by a tile.)`):
print(`The weight of a tiling is the product of H[tile] over all tiles in the tilings.`):
print(`For example, try: GFtH(2,[Rt(1,2),Rt(2,1)],t,H);`):

elif nargs=1 and args[1]=GFtMD then
print(`GFtMD(k,t,h,v,m): The generating function, in the variable t,`):
print(`whose coeffs. of`):
print(`t^n is the weight-enumerator of tilings of the rectangular k by n board with`):
print(`(k-side parallel to the y-axis, n-side parallel to the x-axis) with the horizontal dimers, vertical dimers, and monomers`):
print(`where the weight is h^(NumberOfHorizontalDimers)*v^(NumberOfVerticalDimers)*m^(NumberofMonomers)`):
print(`For example, try: GFtMD(2,t,h,v,m); `):

elif nargs=1 and args[1]=Growth1 then
print(`Growth1(f,t): the growth of the coeff. of the rational function f`):
print(`try for example: Growth1(1/(1-t-t^2),t);`):

elif nargs=1 and args[1]=Kmn then
print(`Kmn(m,n,z): Given a positive integer m and an even positive integers n, computes`):
print(`the Kasteleyn polynomial for tiling`):
print(`an m by n board by domino tiles computed`):
print(`according to Kasteleyn's fromula as it appears in`):
print(`his paper "The statistics of Dimers on a Lattice I...."`):
print(`Physica 27 (1961) Eq. (13) p. 1215`):
print(`For example, try: Kmn(3,4,z);`):


elif nargs=1 and args[1]=ProveKasteleyn then
print(`ProveKasteleyn(m): Given an even integer m,`):
print(`proves Kasteleyn's celebrated formula RIGOROUSLY for any specific even m`):
print(`but ALL n. For m>=8 it may take a very long time.`):
print(`For example, try: ProveKasteleyn(4);`):

elif nargs=1 and args[1]=ProveKasteleyn1 then
print(`ProveKasteleyn1(m): Given an even integer m,`):
print(`proves the straight-enumeration version (z=1)`):
print(`of Kasteleyn's celebrated formula RIGOROUSLY for any specific even m`):
print(`but ALL n. For m>=10 it may take a very long time.`):
print(`For example, try: ProveKasteleyn(4);`):

elif nargs=1 and args[1]=RandTiling then
print(`RandTiling(Tiles,k,n0): a (unif.) random tiling of the k by n0 board `);
print(`using the tiles of Tiles, in set-partition form.`):
print(`For example, try: RandTiling([Rt(1,2),Rt(2,1)],2,2); .`):

elif nargs=1 and args[1]=RandTilingNice then
print(`RandTilingNice(Tiles,k,n0): `):
print(`a (unif.) random tiling of the k by n0 board using the tiles of Tiles.`):
print(`The tiling is given in matrix form.`):
print(`For example, try: RandTilingNice([Rt(1,2),Rt(2,1)],2,2); .`):

elif nargs=1 and args[1]=Rt then
print(`Rt(a,b): Given two positive integers, outputs the rectangular b by a tile`):
print(`{seq(seq([i,j],i=0..b-1),j=0..a-1)}:`):
print(`It is here just to save you typing`):
print(`For example, try: Rt(2,1);`):





elif nargs=1 and args[1]=Sidra then
print(`Sidra(k,Tiles,N): Given a positive integer k, and a list of tiles `):
print(` Tiles, and a positive integer N `):
print(`outputs a list whose (n+1)-th entry is`):
print(`the number of ways of enumerating a k by n rectangle `):
print(` the k-side vetrical, the n-side horiz.`):
print(`The tiles are inputed as sets of lattice points.`):
print(` You can use Rt(a,b) for a b by a rectangular tile `):
print(` For example, try: Sidra(2,[Rt(1,2),Rt(2,1)],20); `):

elif nargs=1 and args[1]=SidraW then
print(`SidraW(k,H,Tiles,N): the first N+1 terms in the`):
print(`weight-enumerating sequence for the number of tilings `):
print(`of an n by k rectangle using tiles from the list `):
print(`of Rectangular tiles Tiles.`):
print(`Here tiles are represented as normalized sets of lattice points.`):
print(` You can use Rt(a,b) for a b by a rectangle `):
print(` and the corresponding weight is H[tile]`):
print(` For example, try: SidraW(2,H,[Rt(1,2),Rt(2,1)],20); `):

elif nargs=1 and args[1]=Sipur then
print(`Sipur(Tiles,t,K,N,Di1): Tells the story for the`);
print(`enumeration of rectangular boards of width up to K`):
print(`using the rectangular tiles Tiles. t is the variable for`):
print(`the length of the board `):
print(`N is a positive integer for the lengths of the`):
print(`desired beginnings for the enumetaring sequences. For example, try: `):
print(` Di1 is the desired number of digits in the floating point (<=100) `):
print(`Sipur([Rt(1,2),Rt(2,1)],t,4,10,10);`):

elif nargs=1 and args[1]=SipurArokh then
print(`SipurArokh(Tiles,t,N,Di1,K1,K2,K3): Tells the long story for the`);
print(`enumeration of rectangular boards of width between K1 and K2, by increment of K3`):
print(`using the rectangular tiles Tiles. t is the variable for`):
print(`the length of the board `):
print(`N is a positive integer for the lengths of the`):
print(`desired beginnings for the enumetaring sequences. For example, try: `):
print(`It also does asymptoic statistics (average nuber of tiles and correlation matrices)`):
print(` Di1 is the desired number of digits in the floating point (<=100) `):
print(`SipurArokh([Rt(1,2),Rt(2,1)],t,10,10,2,4,2):`):



else
 print(`There is no such thing as`, args):

fi:



end:




###############Begin Stuff from StatGF


Digits:=100:
#Asy(f,t,n): Given a rational functionf of the variab;e t, finds the
#asymptotics (in floating point) if the coeff. of t^n
#For example, try: Asy(1/(1-t-^2),t,n);

Asy:=proc(f,t,n) local gu,gu1,ku,i,alufim,shi,muam,lu1,gu2:

if f=0 then
 RETURN(0):
fi:

if normal(subs(t=-t,f)/f)=1 then
 RETURN(FAIL):
fi:

gu:=convert(f,parfrac,t,complex):


if not type(gu,`+`) then
 RETURN(Asy1(gu,t,n)):
fi:

ku:={seq({solve(denom(op(i,gu)),t)}[1],i=1..nops(gu))}:



alufim:={ku[1]}:
shi:=abs(ku[1]):

for i from 2 to nops(ku) do
muam:=abs(ku[i]):

if muam<shi then
  alufim:={ku[i]}:
  shi:=muam:

elif muam=shi then
 alufim:=alufim union {ku[i]}:
fi:
od:

gu1:=0:

for i from 1 to nops(gu) do
lu1:=op(i,gu):

if member({solve(denom(lu1),t)}[1],alufim) then

gu1:=gu1+ lu1:
fi:

od:


if not type(gu1,`+`)  then
 RETURN(Asy1(gu1,t,n)):
fi:

gu2:=0:


for i from 1 to nops(gu1) do

 gu2:=gu2+Asy1(op(i,gu1),t,n):
od:

gu2:

end:


#Asy1(f,t,n): The asymptotics of a single term
Asy1:=proc(f1,t,n) local top1,bot1,a,m:


top1:=numer(f1):
bot1:=denom(f1):


a:={solve(bot1,t)}[1]:



if type(bot1,`^`) then
 m:=op(2,bot1):
if coeff(op(1,bot1),t,1)<>0 then
 top1:=top1/coeff(op(1,bot1),t,1):
else
 RETURN(FAIL):
fi:

else
 m:=1:

if coeff(bot1,t,1)<>0 then
top1:=top1/coeff(bot1,t,1):
 else
 RETURN(FAIL):
fi:

fi:


top1*(-1)^m/a^m*simplify((m+n-1)!/(m-1)!/n!)*(1/a)^n:
end:

#CheckAsy(f,t): checks the correctness of Asy(f,t,n)
#For example, try: CheckAsy(1/(1-t),t);
CheckAsy:=proc(f,t)
local gu,lu,n,i:
gu:=taylor(f,t=0,100):
lu:=Asy(f,t,n):

{seq(coeff(gu,t,i)/subs(n=i,lu),i=80..99)}:

end:


#AsyM(F,t,VarList,ExpList,n): Given F, a rational function of t whose
#coeffs. of t^n are polynomials in the variables of VarList
#corresponding to r.v.
#copmutes the asymptotic value of the Expectation of 
#the product of VarList[i]^ExpList[i]
#For example, to find the asymptotic second moment
#of the number of Heads in tossing n coins type:
#AsyM(1/(1-x*t-y*t^2)),t,[x,y],[2,0],n);
AsyM:=proc(F,t,VarList,ExpList,n)
local gu,i,j,lu,ku:
gu:=F:


for i from 1 to nops(VarList) do
 for j from 1 to ExpList[i] do
   gu:=VarList[i]*diff(gu,VarList[i]):
   gu:=normal(gu):
 od:
od:


gu:=subs({seq(VarList[i]=1,i=1..nops(VarList))},gu):


lu:=subs({seq(VarList[i]=1,i=1..nops(VarList))},F):

ku:=normal(Asy(gu,t,n)/Asy(lu,t,n)):

if degree(ku,n)=FAIL then
 RETURN(FAIL):
fi:

ku:
end:


#Trim2(f,n,L): trims a floating-point expression to L-digit appx.
Trim2:=proc(f,n,L)
local gu,i:

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

gu:=0:

for i from 0 to degree(f,n) do
if abs(coeff(f,n,i))>10^(-L-1) then
 gu:=gu+evalf(coeff(f,n,i),L)*n^i:
fi:
od:
gu:
end:

Trim1:=proc(f,n,L) Trim2(coeff(f,I,1),n,L)*I+ Trim2(coeff(f,I,0),n,L): end:
#Amoment(F,t,VarList,i,k,n): The asymptotic k^th moment about the
#mean  of the r.v. corresponding to the i^th entry in VarList
#where F is a rational function in t whose coeffs. of t^n are
#weight-enumerators according to statistics represented by VarList
#the asymtotic (in n) variance of the number of horizontal tiles
#in a domino-tiling of a 4 by n board, type:
#Amoment(GFtH(4,[Rt(1,2),Rt(2,1)],t,H),t,[H[1,2],H[2,1]],1,2,n):
Amoment:=proc(F,t,VarList,i,k,n) local gu,j,lu,mu,i1:
gu:=[seq(AsyM(F,t,VarList,[0$(i-1),j,0$(nops(VarList)-i)],n),j=1..k)]:

mu:=gu[1]:

lu:=(-1)^k*mu^k:

for i1 from 1 to k do
 lu:=lu+(-1)^(k-i1)*binomial(k,i1)*mu^(k-i1)*gu[i1]:
od:

Trim1(expand(lu),n,40):


end:


#Akurtois(F,t,VarList,i,k,n): The asymptotic Kurtosis 
#of the r.v. corresponding to the i^th entry in VarList
#where F is a rational function in t whose coeffs. of t^n are
#weight-enumerators according to statistics represented by VarList
#For example (assuming that you have CountTilings), to get
#the asymtotic (in n) Kurtosis of the number of horizontal tiles
#in a domino-tiling of a 4 by n board, type:
#AKurtosis(GFtH(4,[Rt(1,2),Rt(2,1)},t,H),t,[H[1,2],H[2,1]],1,n):
Akurtosis:=proc(F,t,VarList,i,n) local m2,m4:
m2:=Amoment(F,t,VarList,i,2,n):
m4:=Amoment(F,t,VarList,i,4,n):

Trim1(normal(coeff(m4,n,degree(m4,n))*n^degree(m4,n)/
(coeff(m2,n,degree(m2,n))*n^degree(m2,n))^2),n,40):

end:





#Aaverage(F,t,VarList,i,n): The asymptotic expectation
#of the r.v. corresponding to the i^th entry in VarList
#where F is a rational function in t whose coeffs. of t^n are
#weight-enumerators according to statistics represented by VarList
#For example (assuming that you have CountTilings), to get
#the asymptotic (in n) average of the number of horizontal tiles
#in a domino-tiling of a 4 by n board, type:
#Aaverage(GFtH(4,[Rt(1,2),Rt(2,1)],t,H),t,[H[1,2],H[2,1]],1,n):
Aaverage:=proc(F,t,VarList,i,n) :
Trim1(AsyM(F,t,VarList,[0$(i-1),1,0$(nops(VarList)-i)],n),n,40):
end:


#Cov(F,VarList,i,j): Given a polynomial F in the variable-list VarList
#finds the covariance between the r.v. associated with VarList[i] and
#VarList[j]. For example, try:
#Cov(x+y^2,[x,y],1,1):
Cov:=proc(F,VarList,i,j) local x,y,gu,av1,av2,var1,var2,lu,i1:
x:=VarList[i]: y:=VarList[j]:
lu:={seq(VarList[i1]=1,i1=1..nops(VarList))}:
gu:=subs(lu,y*diff(x*diff(F,x),y))/subs(lu,F):
av1:=subs(lu,x*diff(F,x))/subs(lu,F):
av2:=subs(lu,y*diff(F,y))/subs(lu,F):
var1:=expand(subs(lu,x*diff(x*diff(F,x),x))/subs(lu,F)-av1^2):
var2:=expand(subs(lu,y*diff(y*diff(F,y),y))/subs(lu,F)-av2^2):
normal((gu-av1*av2)/sqrt(var1*var2)):
end:

#Acovariance(F,t,VarList,i,j,n): The asymptotic covariance
#of the r.v. corresponding to the i^th entry in VarList
#and the j^th entry (1<=i<j<=nops(varList))
#where F is a rational function in t whose coeffs. of t^n are
#weight-enumerators according to statistics represented by VarList
#For example (assuming that you have CountTilings), to get
#the asymptotic (in n) covariance of the number of horizontal tiles
#and the number of vertical tiles
#in a domino-tiling of a 4 by n board, type:
#Acovariance(GFtH(4,[Rt(1,2),Rt(2,1)],t,H),t,[H[1,2],H[2,1]],1,2,n):
Acovariance:=proc(F,t,VarList,i,j,n) local gu,gu1,gu2:
gu:=AsyM(F,t,VarList,[0$(i-1),1,0$(j-i-1),1,0$(nops(VarList)-j)],n):

gu1:=Aaverage(F,t,VarList,i,n) :
gu2:=Aaverage(F,t,VarList,j,n) :
Trim1(expand(gu-gu1*gu2),n,40):
end:


#Acorrelation(F,t,VarList,i,j): The asymptotic correlation
#of the r.v. corresponding to the i^th entry in VarList
#and the j^th entry (1<=i<j<=nops(varList))
#where F is a rational function in t whose coeffs. of t^n are
#weight-enumerators according to statistics represented by VarList
#For example (assuming that you have CountTilings), to get
#the asymptotic  correlation of the number of horizontal tiles
#and the number of vertical tiles
#in a domino-tiling of a 4 by n board, type:
#Acorrelation(GFtH(4,[Rt(1,2),Rt(2,1)],t,H),t,[H[1,2],H[2,1]],1,2):
Acorrelation:=proc(F,t,VarList,i,j) local gu,gu1,gu2,n,guC:
gu:=Acovariance(F,t,VarList,i,j,n):



gu1:=Amoment(F,t,VarList,i,2,n):
gu2:=Amoment(F,t,VarList,j,2,n):



guC:=coeff(gu,n,degree(gu,n)):

if guC=FAIL then
 RETURN(FAIL):
fi:
if evalf(abs(guC))<10^(-Digits/4) then
guC:=coeff(gu,n,degree(gu,n)-1):

if degree(gu1,n)+degree(gu2,n)-2*(degree(gu,n)-1)<>0 then
   RETURN(FAIL):
fi:
fi:



gu1:=sqrt(coeff(gu1,n,degree(gu1,n))):
gu2:=sqrt(coeff(gu2,n,degree(gu2,n))):


if gu1=0 or gu2=0 then
 RETURN(FAIL):
fi:

Trim1(evalf(guC/gu1/gu2,20),n,20):
end:


#aCorMatrix(F,t,VarList): The asymptotic correlation Matrix
#of the random variables corresponding to the r.v. corresponding
#to those in  VarList
#where F is a rational function in t whose coeffs. of t^n are
#weight-enumerators according to statistics represented by VarList
#For example (assuming that you have CountTilings), to get
#the asymptotic correlation of the number of horizontal tiles
#and the number of vertical tiles
#in a domino-tiling of a 4 by n board, type:
#aCorMatrirx(GFtH(4,[Rt(1,2),Rt(2,1)],t,H),t,[H[1,2],H[2,1]],n):
aCorMatrix:=proc(F,t,VarList) local T,i,j:
for i from 1 to nops(VarList) do
 for j from i+1 to nops(VarList) do
  T[i,j]:=Acorrelation(F,t,VarList,i,j):
 od:
od:
[seq(  
[seq(T[j,i],j=1..i-1),1,seq(T[i,j],j=i+1..nops(VarList))]
 , i=1..nops(VarList))]:

end:




###############End Stuff from StatGF





#NormalT1(S): The normal form of a "tile"
NormalT1:=proc(S)
local a,b,i,S1:
a:=min(seq(S[i][1],i=1..nops(S))):

S1:={}:
for i from 1 to nops(S) do
 if S[i][1]=a then
   S1:=S1 union {S[i]}:
 fi:
od:

b:=min(seq(S1[i][2],i=1..nops(S1))):

{seq([S[i][1]-a,S[i][2]-b],i=1..nops(S))}:

end:


NormalT:=proc(S) local i:
{seq(NormalT1(S[i]),i=1..nops(S))}:
end:


#FreeCells(L,N): Given a list of 0-1 lists describing the
#jagged edge of a board of width nops(L) finds the
#set of unoccopied cells in its extension to an N board
#(assumed to be in the positive quadrant of the discrete rectangular lattice)
#for example, try: FreeCells([[0,1],[1,0]],3);
FreeCells:=proc(L,N)
local gu,i,j,k:
k:=nops(L):
gu:={seq(seq([i,j],i=0..N-1),j=0..k-1)}:

for i from 1 to nops(L) do
for j from 1 to nops(L[i]) do
 if L[i][j]=1 then
  gu:=gu minus {[j-1,i-1]}:
 fi:
od:
od:
gu:
end:



#FundamentalCell(S): The fundamental cell of the set of lattice points S
#Try for example: FundamentalCell({[1,0],[0,1],[0,2]});
FundamentalCell:=proc(S)
local a,b,i,S1:
a:=min(seq(S[i][1],i=1..nops(S))):

S1:={}:
for i from 1 to nops(S) do
 if S[i][1]=a then
   S1:=S1 union {S[i]}:
 fi:
od:

b:=min(seq(S1[i][2],i=1..nops(S1))):

[a,b]:

end:

ExtentX:=proc(tile) local i:
max(seq(tile[i][1],i=1..nops(tile)))-
min(seq(tile[i][1],i=1..nops(tile))):
end:

Rt:=proc(a,b) local i,j:
{seq(seq([i,j],i=0..b-1),j=0..a-1)}:
end:



#SetToL1(S): Given a set of points with the same
#y-coordinate 
#to a list of 0 and 1's of the support
#For example, try: SetToL1({[0,3],[1,3],[4,3]});
SetToL1:=proc(S)
local i,m,y,gu:

if S={} then
 RETURN([]):
fi:

if nops({seq(S[i][2],i=1..nops(S))})<>1 then
 RETURN(FAIL):
fi:

y:=S[1][2]:

m:=max(seq(S[i][1],i=1..nops(S))):

gu:=[]:

for i from 0 to m do
 if member([i,y],S) then
  gu:=[op(gu),0]:
 else
  gu:=[op(gu),1]:
 fi:
od:
Chop1(gu):
end:

#Chop1(resh): chops the trailing zeroes of a list
Chop1:=proc(resh) local i:

for i from 1 to nops(resh) while resh[nops(resh)-i+1]=0 do od:

[op(1..nops(resh)-i+1,resh)]:
end:


#SetToL(S,k): Given a set of lattice points in 0<=y<=k-1
#converts it to list notation
SetToL:=proc(S,k) local gu,T,i,K:

for i from 0 to k-1 do
 T[i]:={}:
od:

for i from 1 to nops(S) do
T[S[i][2]]:=T[S[i][2]] union {S[i]}:
od:

gu:=[seq(SetToL1(T[i]),i=0..k-1)]:

K:=max(seq(nops(gu[i]),i=1..nops(gu))):
[seq([op(gu[i]),0$(K-nops(gu[i]))],i=1..nops(gu))]:
end:




ExtractOnes1:=proc(L) local i:

if member(0,{seq(nops(L[i]),i=1..nops(L))}) then
  RETURN(FAIL):
fi:

if {seq(L[i][1],i=1..nops(L))}={1} then
 RETURN([seq([op(2..nops(L[i]),L[i])],i=1..nops(L))]):
else
 RETURN(FAIL):
fi:
end:

ExtractOnes:=proc(L) local L1,co:
L1:=L:
co:=0:
while ExtractOnes1(L1)<>FAIL do
co:=co+1:
L1:=ExtractOnes1(L1):
od:
co,L1:
end:




#Eq1(L,t,F,H,Tiles): Finds the equation for the set of jagged rectangular
#boards of width k (say) encoded by the k-list L, and expressed
#as F[L] in terms of its "children" and the weight enumerator where
#the weight is t^(length) times the product of H[tile] ranging over all
#tiles that participate
#Tiles is the list of tiles (each a set of lattice points)
#it also returns the set of variables used
#For example, try:
#Eq1([[],[]],t,F,H,[ {[0,0],[1,0]},{[0,0],[0,1]} ] );
Eq1:=proc(L,t,F,H,Tiles) local gu,k,eq,N,fc,tile,lu,i,j,pt,gu1,ku,mu:
k:=nops(L):
 if L=[[]$k] then
  eq:=F[L]-1:
 else
  eq:=F[L]:
 fi:

N:=max(seq(ExtentX(Tiles[i]),i=1..nops(Tiles)))+
max(seq(nops(L[i]),i=1..nops(L)))+10:

gu:=FreeCells(L,N):
fc:=FundamentalCell(gu):


mu:={}:

for i from 1 to nops(Tiles) do
tile:=Tiles[i]:
 for j from 1 to nops(tile) do
   pt:=tile[j]:
  
  lu:={seq([tile[k][1]-pt[1]+fc[1],tile[k][2]-pt[2]+fc[2]],k=1..nops(tile))}:
 
 if lu minus  gu={} then
  gu1:=gu minus lu:
  gu1:=SetToL(gu1,k):
  ku:=ExtractOnes(gu1):
  eq:=eq-H[tile]*t^ku[1]*F[ku[2]]:
   mu:= mu union {ku[2]}:
 fi:
 
 od:
od:

eq,mu:
end:








#Yeladim(L,Tiles): Finds the children of a jagged board L
#of width k (say) encoded by the k-list L, For example, try:
#The children are of the form [L1,tile,pivot,generation_below]
#hwere L1 is the new jagged board, tile is the ti that is used
#pivot is the celln the tile coinciding with the Fundamental Cell of L
#and generation_below is a non-negative integer describing  how many
#generations below it is (usually 0). For example, try:
#Yeladim([[],[]],[ Rt(1,2),Rt(2,1) ] );
Yeladim:=proc(L,Tiles) local gu,k,yel,N,fc,tile,lu,i,j,pt,gu1,ku:
k:=nops(L):

yel:={}:

N:=max(seq(ExtentX(Tiles[i]),i=1..nops(Tiles)))+
max(seq(nops(L[i]),i=1..nops(L)))+10:

gu:=FreeCells(L,N):
fc:=FundamentalCell(gu):



for i from 1 to nops(Tiles) do
tile:=Tiles[i]:
 for j from 1 to nops(tile) do
   pt:=tile[j]:
  
  lu:={seq([tile[k][1]-pt[1]+fc[1],tile[k][2]-pt[2]+fc[2]],k=1..nops(tile))}:
 
 if lu minus  gu={} then
  gu1:=gu minus lu:
  gu1:=SetToL(gu1,k):
  ku:=ExtractOnes(gu1):
  yel:=yel union {[ku[2],tile,pt,ku[1]]}:
 fi:
 
 od:
od:

convert(yel,list):
end:


#IniTree(k,Tiles): the initial partial tree for tilings of a rectangular borad
#of width k using the tiles from the list Tiles
#trees are represented as a LIST of the form [parent, SetOfChildrenWithInfo]
IniTree:=proc(k,Tiles)
local L,yel:
L:=[[]$k]:
yel:=Yeladim(L,Tiles):
[[L,yel]]:
end:


#NotYet1(ets): Given a partial tree, finds the set of children that are not yet parents
NotYet1:=proc(ets) local i,j:
{seq(
seq(ets[i][2][j][1],j=1..nops(ets[i][2])),
i=1..nops(ets))}
minus {seq(ets[i][1],i=1..nops(ets))}:
end:

#BuildTree1(k,Tiles): Builds a tree
BuildTree1:=proc(k,Tiles1)
local T,khadas,Tiles:
option remember:
Tiles:=NormalT(Tiles1):
T:=IniTree(k,Tiles):


while NotYet1(T)<>{} do
 khadas:=NotYet1(T)[1]:
  T:=[op(T),[khadas, Yeladim(khadas,Tiles)]]:
od:
T:
end:


#BuildTreeH(k,Tiles,t,H): The tree together with the generating functions in H
BuildTreeH:=proc(k,Tiles,t,H)
local T1,eq,var,L,F,i,j,yel:
option remember:
T1:=BuildTree1(k,Tiles):

L:=T1[1][1]:
yel:=T1[1][2]:
eq:={F[L]-1-add(F[yel[i][1]]*H[yel[i][2]]*t^yel[i][4],i=1..nops(yel))}:

for j from 2 to nops(T1) do
L:=T1[j][1]:
yel:=T1[j][2]:
eq:=eq union {F[L]-add(F[yel[i][1]]*H[yel[i][2]]*t^yel[i][4],i=1..nops(yel))}:
od:

var:={seq(F[T1[i][1]],i=1..nops(T1))}:

var:=solve(eq,var):

if var=NULL then
 RETURN(FAIL):
fi:

[seq([T1[i][1],T1[i][2] , subs(var,F[T1[i][1]]) ],i=1..nops(T1)) ]
end:

#BuildTree(k,Tiles,t): The tree together with the generating functions in t only
BuildTree:=proc(k,Tiles,t)
local T1,eq,var,L,F,i,j,yel:
option remember:
T1:=BuildTree1(k,Tiles):

L:=T1[1][1]:
yel:=T1[1][2]:
eq:={F[L]-1-add(F[yel[i][1]]*t^yel[i][4],i=1..nops(yel))}:

for j from 2 to nops(T1) do
L:=T1[j][1]:
yel:=T1[j][2]:
eq:=eq union {F[L]-add(F[yel[i][1]]*t^yel[i][4],i=1..nops(yel))}:
od:

var:={seq(F[T1[i][1]],i=1..nops(T1))}:

var:=solve(eq,var):

if var=NULL then
 RETURN(FAIL):
fi:


[seq([T1[i][1],T1[i][2] , subs(var,F[T1[i][1]]) ],i=1..nops(T1)) ]

end:



#GFt(k,Tiles,t): The generating function, in the variable t, whose coeffs. of
#t^n is the number of tilings of the rectangular k by n board with
#(k parallel to y-axis, n parallel to x-axis) with the tiles of the list Tiles
#where ar tiles are described as normalized sets of lattice points
#For example, try: GFt(2,[Rt(1,2),Rt(2,1)],t);

GFt:=proc(k,Tiles,t):  BuildTree(k,Tiles,t)[1][3]: end:

#GFtH(k,Tiles,t,H): The generating function, in the variable t, whose coeffs. of
#t^n is the weight-enumerator of tilings of the rectangular k by n board with
#(k parallel to y-axis, n parallel to x-axis) with the tiles of the list Tiles
#where ar tiles are described as normalized sets of lattice points
#the weight of a tiling is the product of H[tile] over all tiles in the tilings
#For example, try: GFtH(2,[Rt(1,2),Rt(2,1)],t,H);
GFtH:=proc(k,Tiles,t,H):  BuildTreeH(k,Tiles,t,H)[1][3]: end:

#GFtMD(k,t,h,v,m):  the generating function for monomer-dimer tilings where h=horiz. dimer, v=vertical dimer, m=monomer
GFtMD:=proc(k,t,h,v,m) local gu: 
gu:=GFtH(k,[Rt(1,2),Rt(2,1),Rt(1,1)],t,H): 
normal(subs({H[Rt(1,2)]=h,H[Rt(2,1)]=v, H[Rt(1,1)]=m},gu)):
end:



#LoadedDie(resh): picks an integer from 1 to nops(resh) according to the loads of the
#list of positive integers resh:
LoadedDie:=proc(resh)
local L,ra,j:
L:=convert(resh,`+`):
if L=0 then
 RETURN(FAIL):
fi:
ra:=rand(1..L)():


if ra<=resh[1] then
 RETURN(1):
fi:

for j from 2 to nops(resh)-1 do
if ra>convert([op(1..j-1,resh)],`+`) and ra<=convert([op(1..j,resh)],`+`) then
 RETURN(j):
fi:
od:
nops(resh):

end:


#RandTi1(T,t,v,Tiles,k,n0): Given a tree T a vertex v and an integer n
#outputs a (unif.) random tiling of the k by n rect. board
#for example, try: RandTi1( T2,t,[[],[]],[Rt(1,2),Rt(2,1)],2,3);
RandTi1:=proc(T,t,v,Tiles,k,n0) local gu,i,Tab1,yel,gf1,mu,j,mu1,v1,tile1,pivot1,shift1,RandT1,fc,N:

N:=max(seq(ExtentX(Tiles[i]),i=1..nops(Tiles)))+
max(seq(nops(v[i]),i=1..nops(v)))+10:
gu:=FreeCells(v,N):
fc:=FundamentalCell(gu):


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

for i from 1 to nops(T) do
Tab1[T[i][1]]:=T[i][3]:
od:


if not member(v,gu) then
 RETURN(FAIL):
fi:

if n0=0 then
if v=[[]$k] then
 RETURN({}):
else

 RETURN(FAIL):
fi:
fi:

for i from 1 to nops(T) do
if T[i][1]=v then
 break:
fi:
od:
yel:=T[i][2]:
gf1:=T[i][3]:


mu:=[]:
for j from 1 to nops(yel) do

 mu:=[op(mu),
 coeff(taylor(Tab1[yel[j][1]],t=0,n0+1),t,n0-yel[j][4])]:
od:

mu1:=coeff(taylor(gf1,t=0, n0+1),t,n0):

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

if convert(mu,`+`)<>mu1 then
 print(`Something is wrong`):
 RETURN(FAIL):
fi:

j:=LoadedDie(mu):

v1:=yel[j][1]:
tile1:=yel[j][2]:
pivot1:=yel[j][3]:
shift1:=yel[j][4]:


RandT1:=RandTi1(T,t,v1,Tiles,k,n0-shift1):


RandT1:={seq(
{seq([RandT1[i][j][1]+shift1,RandT1[i][j][2]],j=1..nops(RandT1[i]))}
,i=1..nops(RandT1))}:

tile1:={seq([tile1[i][1]-pivot1[1]+fc[1],tile1[i][2]-pivot1[2]+fc[2]],i=1..nops(tile1))}:


{op(RandT1),tile1}:

end:




#RandTiling(Tiles,k,n0): a (unif.) random tiling of the k by n0 board 
#using the tiles of Tiles, in set-partition form
#For example, try: RandTiling([Rt(1,2),Rt(2,1)],2,2);
RandTiling:=proc(Tiles,k,n0) local T,t:
T:=BuildTree(k,Tiles,t):
RandTi1(T,t,[[]$k],Tiles,k,n0):
end:


#TileToMat(Ti,k0,n0): Given a tiling of the k0 by n0 board as a set partition
#returns it as a matrix
TileToMat:=proc(Ti,k0,n0) local M,i,j,vu,T,M1,co:
if Ti=FAIL then
 RETURN(FAIL):
fi:
M:=matrix(k0,n0):

for i from 1 to nops(Ti) do
 for j from 1 to nops(Ti[i]) do
  M[Ti[i][j][2]+1,Ti[i][j][1]+1]:=i:
 od:
od:

vu:={}:
co:=0:
for j from 1 to n0 do
for i from 1 to k0 do
 if not member(M[i,j],vu) then
   vu:=vu union {M[i,j]}:
  co:=co+1:
  T[M[i,j]]:=co:
 fi:
od:
od:

M1:=matrix(k0,n0):

for i from 1 to k0 do
 for j from 1 to n0 do
  M1[i,j]:=T[M[i,j]]:
 od:
od:
op(M1):
end:




#RandTilingNice(Tiles,k,n0): 
#a (unif.) random tiling of the k by n0 board using the tiles of Tiles
#in Matrix form
#For example, try: RandTilingNice([Rt(1,2),Rt(2,1)],2,2);
RandTilingNice:=proc(Tiles,k,n0) :
TileToMat(RandTiling(Tiles,k,n0),k,n0):
end:









#AllTilings(Tiles,k,n0): The set of tilings of the k by n0 board 
#using the tiles of Tiles, in set-partition form
#For example, try: AllTilings([Rt(1,2),Rt(2,1)],2,2);
AllTilings:=proc(Tiles,k,n0) local T,t:
T:=BuildTree(k,Tiles,t):
AllTi1(T,t,[[]$k],Tiles,k,n0):
end:


#AllTilingsNice(Tiles,k,n0): The set of tilings of the k by n0 board 
#using the tiles of Tiles, in matrix form and listed one-by-one
#For example, try: AllTilingsNice([Rt(1,2),Rt(2,1)],2,2);
AllTilingsNice:=proc(Tiles,k,n0) local gu,i:
gu:=AllTilings(Tiles,k,n0):

print(`There are`, nops(gu), `tilings using the tiles`, Tiles, `of the`, k, ` by `, n0 , ` rectangular board `):
print(`Here they are: `):
for i from 1 to nops(gu) do

print(TileToMat(gu[i],k,n0)):
od:
end:



#AllTi1(T,t,v,Tiles,k,n0): Given a tree T a vertex v and an integer n
#outputs all tilings of the k by n rect. board
#for example, try: AllTi1( T2,t,[[],[]],[Rt(1,2),Rt(2,1)],2,3);
AllTi1:=proc(T,t,v,Tiles,k,n0) local gu,i,yel,mu,j,mu1,v1,tile1,pivot1,shift1,fc,N,
i3,j3,k1,kha:


if n0<0 then
 RETURN({}):
fi:

if n0=0 then
if v=[[]$k] then
 RETURN({{}}):
else

 RETURN({}):
fi:
fi:

N:=max(seq(ExtentX(Tiles[i]),i=1..nops(Tiles)))+
max(seq(nops(v[i]),i=1..nops(v)))+10:
gu:=FreeCells(v,N):
fc:=FundamentalCell(gu):


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


if not member(v,gu) then
 RETURN(FAIL):
fi:


for i from 1 to nops(T) do
if T[i][1]=v then
 break:
fi:
od:
yel:=T[i][2]:

mu:={}:

for j from 1 to nops(yel) do

v1:=yel[j][1]:
tile1:=yel[j][2]:
pivot1:=yel[j][3]:
shift1:=yel[j][4]:


mu1:=AllTi1(T,t,v1,Tiles,k,n0-shift1):

tile1:={seq([tile1[i][1]-pivot1[1]+fc[1],tile1[i][2]-pivot1[2]+fc[2]],i=1..nops(tile1))}:

for k1 from 1 to nops(mu1) do
 kha:=mu1[k1]:
kha:={seq(
{seq([kha[i3][j3][1]+shift1,kha[i3][j3][2]],j3=1..nops(kha[i3]))}
,i3=1..nops(kha))}:

mu:= mu union {{op(kha),tile1}}:
od:

od:
mu:
end:


#T2:=BuildTree(2,[Rt(1,2),Rt(2,1)],t):

ez:=proc():print(` NormalT1(S), NormalT(S), Eq1(L,t,F,H,Tiles), FreeCells(L,N), FundamentalCell(S) `):
print(`Eq1([[],[]],t,F,H,[Rt(1,2),Rt(2,1) ] );`):
print(`Rt(a,b), SetToL1({[0,3],[1,3],[4,3]}); `):
print(`SetToL(S,k), Chop1(resh) , ExtractOnes1(L) , Eq(k,t,F,H,Tiles) `):
print(`Yeladim([[],[]],[ Rt(1,2),Rt(2,1) ] `);
print(`IniTree(k,Tiles), NotYet1(ets)`):
print(`BuildTree1(k,Tiles)`):
print(`BuildTreeH(k,Tiles,t,H), BuildTree(k,Tiles,t) , GFt(k,Tiles,t) `):
print(`RandTi1( T2,t,[[],[]],[Rt(1,2),Rt(2,1)],2,3);`):
print(`LoadedDie(resh), RandTi(Tiles,k,n0);`):
print(`TileToMat(Ti,k0,n0) `):
print(`RandTiNice(Tiles,k,n0):`): 
print(`GFtH(k,Tiles,t,H)`):
print(`AllTilings(Tiles,k,n0), AllTi1(T,t,v,Tiles,k,n0) `):
print(`AllTilingsNice(Tiles,k,n0)`):
end:




####Salvaged Stuff


#Astat(k,ListOfTiles,n): Given a List of rectangular Tiles
#ListOfTiles
#returns: (1) The lists of for the
#asymptotics (in n) for pairs [average,variance] for a random
#tiling by the rectangular tiles
#of Tiles (in order), 
#followed by the av. and variance for the
#Total number of tiles, followed by the asympto. correlation
#matrix for example try:
#Astat(4,[Rt(1,2),Rt(2,1)],n):
Astat:=proc(k,Tiles,n) local t,H,F,VarList,i,j,s,Fs,mu,gu:

VarList:=[seq(H[Tiles[i]],i=1..nops(Tiles))]:

F:=GFtH(k,Tiles,t,H):


Fs:=subs({seq(H[Tiles[i]]=s,i=1..nops(Tiles))},F):


mu:=aCorMatrix(F,t,VarList):


for i from 1 to nops(mu) do
 for j from 1 to nops(mu[i]) do
 if mu[i][j]=FAIL then
   RETURN(FAIL):
 fi:
 od:
od:


gu:=[[seq([Aaverage(F,t,VarList,i,n) ,Amoment(F,t,VarList,i,2,n)],i=1..nops(VarList))],
[Aaverage(Fs,t,[s],1,n) ,Amoment(Fs,t,[s],1,2,n) ],mu]:


if member(FAIL, {seq(gu[1][i][1],i=1..nops(VarList)),seq(gu[1][i][2],i=1..nops(VarList))}) then
 RETURN(FAIL):
fi:

if member(FAIL, convert(gu[2],set)) then
 RETURN(FAIL):
fi:

gu:
end:

#The list of asymptotic averages of the number of tiles
#of tiling a k by n rectangular board by the  tiles in
#the list of rectangular tiles Tiles ([a,b] means a b by a rectangle)
#For example, try: AaveT(2,[[1,2],[2,1]]) 
#followed by the total number of tiles
AaveT:=proc(k,Tiles) local t,H,F,VarList,i,s,Fs,n,gu:

VarList:=[seq(H[Tiles[i]],i=1..nops(Tiles))]:

F:=GFtH(k,Tiles,t,H):

Fs:=subs({seq(H[Tiles[i]]=s,i=1..nops(Tiles))},F):

gu:=[seq(coeff(Aaverage(F,t,VarList,i,n),n,1) ,i=1..nops(VarList)), coeff(Aaverage(Fs,t,[s],1,n),n,1) ]:

if abs(convert([op(1..nops(gu)-1,gu)],`+`)-gu[nops(gu)])>10^(-10) then
 print(`Couldn't compute averages`):
 RETURN(FAIL,gu):
fi:

gu:

end:


#AstatV(k,ListOfTiles,n): verbose version of Astat
#For example try:
#AstatV(4,[Rt(1,2),Rt(2,1)],n):
AstatV:=proc(k,Tiles,n) local gu,i:

gu:=Astat(k,Tiles,n):


print(`Consider a random tiling of the rectangular board of dimension`, k, ` by n `):
print(k, `  is the vertical dimension and n the horizontal one`):
print(`Using the tiles in the list of tiles`):
print(Tiles):
print(`We are interested in the asymptotics as n goes to infinity`):
print(`of the statistics for the number of tiles`):

for i from 1 to nops(Tiles) do
 print(`The asymptotic (in n)  average number of tiles of type`, Tiles[i] , `in such a random tiling is`):
 print(evalf(coeff(gu[1][i][1],n,1 ),10)*n, `[variance=`, evalf(coeff(gu[1][i][2],n,1),10)*n, `]` ):
# print(`and the variance is`, evalf(coeff(gu[1][i][2],n,1),10)*n):
od:

print(`The average total number of tiles is`):
print(evalf(coeff(gu[2][1],n,1),10)*n, `[variance=`, evalf(coeff(gu[2][2],n,1),10)*n, `]` ):
#print(`and the variance is`, evalf(coeff(gu[2][2],n,1),10)*n ):
print(`Finally, the asymptotic correlation matrix, in the order of appearace of`, Tiles, ` is `):
print(evalf(gu[3],10)):
end:



AstatVs:=proc(K,Tiles,n) local i:
print(`This is info about statistics for radnom tilings of rectangular board`):
print(`of width from 2 to`, K):

for i from 2 to K do
 AstatV(i,Tiles,n):
od:
end:



#Sidra(k,Tiles,N): the first N+1 terms in the
#enumerating sequence for the number of tilings using
#of an n by k rectangle using tiles from the set
Sidra:=proc(k,Tiles,N) local gu,t,i:
gu:=GFt(k,Tiles,t):
[seq(coeff(taylor(gu,t=0,N+1),t,i),i=0..N)]:
end:

#SidraW(k,H,Tiles,N): the first N+1 terms in the
#weight-enumerating sequence for the number of tilings using
#of an n by k rectangle using tiles from the set
#For example, try: SidraW(2,H,{[1,2],[2,1]},20);
SidraW:=proc(k,H,Tiles,N) local gu,t,i:
gu:=GFtH(k,Tiles,t,H):
[seq(expand(coeff(taylor(gu,t=0,N+1),t,i)),i=0..N)]:
end:

Sipur:=proc(Tiles,t,K,N,Di1) local k,lu,gu,vu,Di2:
if Di1>100 then
 print(`I only do precison with <=100 digits `):
Di2:=100:
else
 Di2:=Di1:
fi:

print(`This is the story of the number of tilings of  rectangular`):
print(`boards of width <=`,K, `using rectangular tiles of dimensions`):
print(Tiles):
vu:=[]:

for k from 1 to K do
 print(`The number of tilings using the tiles`):
 print(Tiles):
 print(`of a `, k, ` by n rectangular board is the coeff. of t^n in the rational function`):
 gu:=GFt(k,Tiles,t):
 lprint(gu):
 print(`The first `, N+1, `terms are `):
 print(Sidra(k,Tiles,N)):
 lu:=Growth1(gu,t):
 print(`The asymptotic rate of growth is`,evalf(lu,Di2)):
 print(`and adjusted (i.e. the `, k, ` -th root) it is`, evalf(lu^(1/k),Di2)):
 vu:=[op(vu),lu^(1/k)]:
 print(``):
 print(``):
od:
 print(`To summarize, the sequence of adjusted growth-rates `):
 print(`For tilings with the set of tiles`, Tiles , ` is `):
 print(evalf(vu,Di2)):
vu:
end:


SipurArokh:=proc(Tiles,t,N,Di1,K1,K2,K3) local k,lu,gu,vu,Di2,n:
if Di1>100 then
 print(`I only do precison with <=100 digits `):
Di2:=100:
else
 Di2:=Di1:
fi:

print(`This is the long story of the number of tilings of  rectangular`):
print(`boards of width bewteen`,K1, `and  `, K2, ` by inrcements of`, K3):
print( `using tiles from the following list:`):
print(Tiles):
vu:=[]:
print(`-------------------------------------------------------`):
for k from K1 to K2 by K3 do
 print(`The number of tilings using the tiles`):
 print(Tiles):
 print(`of a `, k, ` by n rectangular board is the coeff. of t^n in the Maclaurin series of the the rational function`):
 gu:=GFt(k,Tiles,t):
 lprint(gu):
 print(`The first `, N+1, `terms are `):
 print(Sidra(k,Tiles,N)):
 lu:=Growth1(gu,t):
 print(`The asymptotic rate of growth is`,evalf(lu,Di2)):
 print(`and adjusted (i.e. the `, k, ` -th root)  is`, evalf(lu^(1/k),Di2)):
 vu:=[op(vu),lu^(1/k)]:
 print(``):
print(`The asymptotics for the number of tilings as n goes to infinity is`):
 print(coeff(evalf(Asy(gu,t,n) ,Di1),I,0)):
print(``):

 if Astat(k,Tiles,n)<>FAIL then
print(` We now pause for asymptotics statistics for tilings of the set of the tiles in`):
print(Tiles):
print(`For rectangles of width`, k):
print(`------------------------------------------------------------`):
 AstatV(k,Tiles,n):
print(`------------------------------------------------------------`):
print(` this ends the  asymptotics statistics`):
print(` for tilings of the set of the tiles in`):
print(Tiles):
print(`For rectangles of width`, k):
 print(``):
print(`----------------------------------------------`):
fi:
od:
 print(`To summarize, the sequence of adjusted growth-rates `):
 print(`for tilings with the set of tiles`, Tiles , ` is `):
 print(evalf(vu,Di2)):
vu:
end:



#Growth1(f,t): the growth of the coeff. of the rational function f
#try for example: Growth1(1/(1-t-t^2),t);
Growth1:=proc(f,t) local i,aluf,ku,shi,P:
if normal(subs(t=-t,f)/f)=1 then
 RETURN(FAIL):
fi:
P:=denom(normal(f)):
if degree(P,t)=0 then
 RETURN(FAIL):
fi:
ku:={fsolve(P,t)}:
aluf:=ku[1]:
shi:=abs(aluf):

for i from 2 to nops(ku) do
 if abs(ku[i])<shi then
    aluf:=ku[i]:
    shi:=abs(ku[i]):
 fi:
od:
1/shi:
end:



#Kmn(m,n,z): The Kasteleyn polynomial for tiling
#an m by n board by domino tiles computed
#according to Kasteleyn's fromula as it appears in
#his paper "The statistics of Dimers on a Lattice I...."
#Physica 27 (1961) Eq. (13) p. 1215
#Using Floating point
#For example, try: Kmn(3,4,z);
Kmn:=proc(m,n,z) local gu,K,L,i:
Digits:=100:
gu:=1:


if n mod 2=1 then
 RETURN(FAIL):
fi:

gu:=2^(n*trunc(m/2)):


 for K from 1 to m/2 do
 for L from 1 to n/2 do
   gu:=gu*(z*cos(K*Pi/(m+1))^2+cos(L*Pi/(n+1))^2):
 od:
 od:

gu:=add(round(evalf(coeff(gu,z,i)))*z^i,i=0..degree(gu,z)):

RETURN(gu):

end:

#Kmn1(m,n): The Kasteleyn polynomial for tiling (with z=1)
#an m by n board by domino tiles computed
#according to Kasteleyn's fromula as it appears in
#his paper "The statistics of Dimers on a Lattice I...."
#Physica 27 (1961) Eq. (13) p. 1215
#Using Floating point
#For example, try: Kmn1(3,4);
Kmn1:=proc(m,n) local gu,K,L,i:
Digits:=100:
gu:=1:


if n mod 2=1 then
 RETURN(FAIL):
fi:

gu:=2^(n*trunc(m/2)):


 for K from 1 to m/2 do
 for L from 1 to n/2 do
   gu:=gu*(cos(K*Pi/(m+1))^2+cos(L*Pi/(n+1))^2):
 od:
 od:

gu:=round(gu):

RETURN(gu):

end:


#ProveKasteleyn(m): Given an even integer m,
#proves Kasteleyn's celebrated formula RIGOROUSLY for any specific even m
#but ALL n. For m>=10 it may take a very long time.
#For example, try: ProveKasteleyn(4);
ProveKasteleyn:=proc(m) local z,H,i:
if m mod 2 =1 then
 print(`Input must be an even integer`):
 RETURN(FAIL):
fi:
evalb([seq(Kmn(i,m,z),i=0..2^(m/2+1)+1 )]=
subs({H[Rt(1,2)]=sqrt(z),H[Rt(2,1)]=1},SidraW(m,H,[Rt(1,2),Rt(2,1)],2^(m/2+1)+1)));
end:
                                                            

#ProveKasteleyn1(m): Given an even integer m,
#proves Kasteleyn's celebrated formula (with z=1) RIGOROUSLY for any specific even m
#but ALL n. For m>=10 it may take a very long time.
#For example, try: ProveKasteleyn(4);
ProveKasteleyn1:=proc(m) local i:
if m mod 2 =1 then
 print(`Input must be an even integer`):
 RETURN(FAIL):
fi:
evalb([seq(Kmn1(i,m),i=0..2^(m/2+1)+1 )]=
Sidra(m,[Rt(1,2),Rt(2,1)],2^(m/2+1)+1));
end:
                                                            





#Dimer(n,z): the C-finite representation of the generating function for the dimer problem
#in an n by m strip (m general, n specific) with vertical variable z and horiz. variable 1, try:
#Dimer(4,z);
Dimer:=proc(n,z) local lu,H,t,i:
lu := subs({H[{[0, 0], [0, 1]}] = z, H[{[0, 0], [1, 0]}] = 1}, GFtH(n, [Rt(1, 2), Rt(2, 1)], t, H)):
lu:=denom(lu):
lu:=lu/coeff(lu,t,0):
[seq(-coeff(lu,t,i),i=1..degree(lu,t))]:
end:





#########END FROM TILINGS
#SeqT(m,k,N): inputs a positive even integer m, a rational number k=a/b, and outputs the sequence of the number of tilings of a m by (n*b) rectangle
#with a/b*(m*n*b)=m*n*a monomers and (m*n*b-m*n*a)/2=m/2*n*(b-a) dimers for n from 1 to N. Try:
#SeqT(2,1/2,40);

SeqT:=proc(m,k,N) local z,w,f,t,a,b,i,n,L,H:
option remember:
a:=numer(k):
b:=denom(k):


f:=GFmdPC(m,t,z,w) :

f:=taylor(f,t=0,N*b+1):

L:=[seq(coeff(f,t,i),i=1..N*b)]:



[seq(coeff(coeff(L[b*n],w,m*n*a),z,m/2*n*(b-a)),n=1..nops(L)/b)]:


end:

#SeqTo(m,k,N): inputs a positive odd integer m, a rational number k=a/b, and outputs the sequence of the number of tilings of a m by (2*n*b) rectangle
#with a/b*(2*m*n*b)=2*m*n*a monomers and (m*2*n*b-m*2*n*a)/2=m*n*(b-a) dimers for n from 1 to N. Try:
#SeqTo(2,1/2,40);

SeqTo:=proc(m,k,N) local z,w,f,t,a,b,i,n,L,H:
option remember:
a:=numer(k):
b:=denom(k):


f:=GFmdPC(m,t,z,w) :

f:=taylor(f,t=0,N*b*2+1):

L:=[seq(coeff(f,t,i),i=1..2*N*b)]:



[seq(coeff(coeff(L[2*b*n],w,2*m*n*a),z,m*n*(b-a)),n=1..nops(L)/b/2)]:


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,flo,pip:
if (1+degx)*(1+degP) > nops(gu)-3 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 nops(gu)-1 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 nops(gu)-2 do
eq:=eq union {coeff(lu,x,i1)=0}
od:
 
mu:=solve(eq,var):
 
 
F:=subs(mu,F):
 
if F=0 then
 RETURN(0):
fi:
 
flo:=degree(F,P):
pip:=coeff(F,P,flo):
flo:=degree(pip,x):
pip:=coeff(pip,x,flo):
F:=F/pip:
normal(F):
end:

#Empir(L,x,P): inputs a sequence of numbers L, starting at n=0 and guesses an algebratic equation for its ogf P(x)
Empir:=proc(gu,x,P)
local degx,degP,L,lu:
 
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<>0 then
RETURN(lu):
fi:
od:
od:
od:
0:
end:

sn:=proc(resh,n1):
-1/log(op(n1+1,resh)*op(n1-1,resh)/op(n1,resh)^2):
end:
 
#Zinn(resh): Zinn-Justin's method to estimate
#the C,mu, and theta such that
#resh[i] is appx. Const*mu^i*i^theta
#For example, try:
#Zinn([seq(5*i*2^i,i=1..30)]);
Zinn:=proc(resh)
local s1,s2,theta,mu,n1,i:
if nops({seq(sign(resh[i]),i=1..nops(resh))})<>1 then
 RETURN(FAIL):
fi:

n1:=nops(resh)-1:
s1:=sn(resh,n1):
s2:=sn(resh,n1-1):
theta:=evalf(2*(s1+s2)/(s1-s2)^2):
mu:=evalf(sqrt(op(n1+1,resh)/op(n1-1,resh))*exp(-(s1+s2)/((s1-s2)*s1))):
[theta,mu]:
end:
 
#MDalg(m,k,N,x,P): inputs an even integer m, a rational number k (=a/b,say), a large pos. integer N, and symbols x,P
#outputs a guessed (but rigorizable) algebraic equation for the ordinary generating function of the sequence
#a(n):=Number of monomer-dimer tilings of an m by b*n rectangle with m*b*n*k= m*n*a monomers and the rest m*n*(b-a)/2 dimers
#Try:
#MDalg(2,1/2,40,x,P);
MDalg:=proc(m,k,N,x,P) local L,eq,i:
L:=SeqT(m,k,N):

L:=[1,op(L)]:

eq:=Empir(L,x,P):


if eq=FAIL then
 print(`Make N bigger`):
 RETURN(FAIL):
fi:

add(factor(coeff(eq,P,i))*P^i,i=0..degree(eq,P)):

end:


#CC(m,k,N): The "connective constant" for the sequence SeqT(m,k). Try:
#CC(2,1/2,40):
CC:=proc(m,k,N) local eq,x,P,min1,sol, cham,i,sol1:
eq:=MDalg(m,k,N,x,P):
eq:=discrim(eq,P):
sol:=[solve(eq,x)]:

sol1:=[]:

for i from 1 to nops(sol) do
 if abs(coeff(evalf(sol[i]),I,1)) <10^(-Digits/2) and coeff(evalf(sol[i]),I,0)>0 then
  sol1:=[op(sol1),sol[i]]:
 fi:
od:


min1:=abs(evalf(sol1[1])):
cham:=1:

for i from 2 to nops(sol1) do
 if abs(evalf(sol1[i]))<min1  then
    cham:=i:
    min1:=abs(evalf(sol1[i])):
 fi:
od:

1/sol1[cham]:

end:


#Story(m,k,N,x,P): tells the story of MDalg(m,k,N,x,P). Try:
#Story(2,1/2,50,x,P):
Story:=proc(m,k,N,x,P) local eq,L,a,b,lu,mu,theta,t0,n,s:

t0:=time():
L:=SeqT(m,k,N):
a:=numer(k):
b:=denom(k):
print(`On the Sequence enuerating Monomer-Dimer Tilings of an `, m ,` by  `, n*b , ` rectangle `):
print(`with Monomer Density`, k):
print(``):
print(`Definition: Let s(n) be the number of Monomer-Dimer Tilings of an `, m ,` by  `, n*b , ` rectangle with monomer density`,k):
print(`in other words, there are `, m*n*a, `monomers and `, m*n*(b-a)/2, `dimers . `):
print(`The first`, N, `terms are `):
print(``):
print(L):

lu:=Zinn(L):

theta:=lu[1]:
mu:=evalf(lu[2],10):
theta:=1/round(1/theta):

print(`The asymptotics of the sequence, numerically is`, n^theta*mu^n):

print(``):

eq:=MDalg(m,k,N,x,P):

 if eq=FAIL then
   print(`This ends this article`):
   RETURN(L):
 fi:

print(`Theorem: Let`):
print(f(x)=Sum( s(n)*x^n,n=0..infinity)):
print(``):
print(`Then f(x)=P satisfies the following algebraic equation`):
print(``):
print(eq=0):
print(``):
print(`and in Maple format `):
print(``):
lprint(eq=0):
print(``):


lu:=CC(m,k,N):

print(`The exact value of the connective constant, derived numerically above is`):
print(lu):
print(``):
print(`and in floating-point`):
print(``):
print(evalf(lu,20)):
print(``):
print(`This ends this article, that took`, time()-t0, `to generate.`):

end:

#fmkN(m,k,N): Let A_m(k,n) be the number of monomer-dimer tilings of an m by b*n rectangle (m even)
#with monomer density k, finds (numerically) the limit of log(A_m(k,n))/(n*m*b) as n goes to infinity
#Try:
#fmkN(2,1/2,100);
fmkN:=proc(m,k,N) local L,gu,a,b:
b:=denom(k):
L:=SeqT(m,k,N):
gu:=Zinn(L)[2]:
log(gu)/m/b:
end:

#fmkNo(m,k,N): Let A_m(k,n) be the number of monomer-dimer tilings of an m by b*n rectangle (m odd)
#with monomer density k, finds (numerically) the limit of log(A_m(k,n))/(n*m*b) as n goes to infinity
#Try:
#fmkNo(2,1/2,100);
fmkNo:=proc(m,k,N) local L,gu,a,b:
b:=denom(k):
L:=SeqTo(m,k,N):
gu:=Zinn(L)[2]:
log(gu)/m/b/2:
end:



#fmkE(m,k,N): Let A_m(k,n) be the number of monomer-dimer tilings of an m by b*n rectangle (m even)
#with monomer density k, finds (numerically) the limit of log(A_m(k,n))/(n*m*b) as n goes to infinity
#Try:
#fmkE(2,1/2,100);
fmkE:=proc(m,k,N) local L,gu,a,b:
b:=denom(k):
L:=SeqT(m,k,N):
gu:=CC(m,k,N):
log(gu)/m/b:
end:


#fmkNdata(m,A,N): a table of [k,fmkN(m,k,N)] for all k=a/b with gcd(a,b)=1, 1<=a<b<=A. Try:
#fmkNdata(2,3,100);
fmkNdata:=proc(m,A,N) local a,b,gu,T,i:

gu:=[0]:
if m=2 then
T[0]:=evalf(log((sqrt(5)+1)/2)/2):
else
T[0]:=fmkN(m,0,N):
fi:


for b from 2 to A do
for a from 1 to b-1 do
if gcd(a,b)=1 then 
T[a/b]:=fmkN(m,a/b,N):
 gu:=[op(gu),a/b]:
fi:
od:
od:

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

end:


#fmkNodata(m,A,N): a table of [k,fmkN(m,k,N)] for all k=a/b with gcd(a,b)=1, 1<=a<b<=A. Try:
#fmkNodata(2,3,100);
fmkNodata:=proc(m,A,N) local a,b,gu,T,i:

gu:=[0]:
if m=2 then
T[0]:=evalf(log((sqrt(5)+1)/2)/2):
else
T[0]:=fmkNo(m,0,N):
fi:


for b from 2 to A do
for a from 1 to b-1 do
if gcd(a,b)=1 then 
T[a/b]:=fmkNo(m,a/b,N):
 gu:=[op(gu),a/b]:
fi:
od:
od:

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

end:


#GFmdPC(m,t,z,w): Pre-computed version of GFmd(m,t,z,w) (q.v) for 1<=m<=6  otherwise does GFmd(m,t,z,w): Try
#GFmdPC(4,t,z,w);

GFmdPC:=proc(m,t,z,w) local L,i:
if m>6 then
 RETURN(GFmd(m,t,z,w)):
else
L:=
[-1/(t^2*z+t*w-1), -(t*z-1)/(t^3*z^3-t^2*w^2*z-t*w^2-2*t*z+1), -(t^4*z^6+t^3*w*
z^4-2*t^2*w^2*z^2-2*t^2*z^3-t*w*z+1)/(t^6*z^9+t^5*w^3*z^6-t^5*w*z^7-2*t^4*w^4*z
^4-3*t^4*w^2*z^5-5*t^4*z^6+t^3*w^5*z^2+t^3*w^3*z^3-2*t^3*w*z^4+2*t^2*w^4*z+7*t^
2*w^2*z^2+5*t^2*z^3+t*w^3+3*t*w*z-1), -(t^7*z^14-2*t^6*w^2*z^11-3*t^5*w^4*z^8-3
*t^5*w^2*z^9-5*t^5*z^10+3*t^4*w^4*z^6+5*t^4*w^2*z^7+2*t^4*z^8+3*t^3*w^4*z^4+11*
t^3*w^2*z^5+6*t^3*z^6-3*t^2*w^4*z^2-9*t^2*w^2*z^3-3*t^2*z^4-2*t*w^2*z-2*t*z^2+1
)/(t^9*z^18-t^8*w^4*z^14+t^8*w^2*z^15-t^8*z^16-2*t^7*w^6*z^11-3*t^7*w^4*z^12-9*
t^7*w^2*z^13-t^6*w^8*z^8-9*t^7*z^14+6*t^6*w^4*z^10+19*t^6*w^2*z^11+3*t^5*w^8*z^
6+5*t^6*z^12+14*t^5*w^6*z^7+29*t^5*w^4*z^8+24*t^5*w^2*z^9-3*t^4*w^8*z^4+21*t^5*
z^10-18*t^4*w^6*z^5-41*t^4*w^4*z^6-40*t^4*w^2*z^7+t^3*w^8*z^2-9*t^4*z^8+4*t^3*w
^6*z^3-4*t^3*w^4*z^4-27*t^3*w^2*z^5-15*t^3*z^6+2*t^2*w^6*z+13*t^2*w^4*z^2+21*t^
2*w^2*z^3+5*t^2*z^4+t*w^4+5*t*w^2*z+3*t*z^2-1), -(t^18*z^45+2*t^17*w^3*z^41-7*t
^16*w^6*z^37-2*t^16*w^4*z^38-19*t^16*w^2*z^39+4*t^15*w^9*z^33-17*t^16*z^40-12*t
^15*w^7*z^34-10*t^15*w^5*z^35-44*t^15*w^3*z^36+14*t^14*w^10*z^30-6*t^15*w*z^37+
29*t^14*w^8*z^31+117*t^14*w^6*z^32+188*t^14*w^4*z^33+10*t^13*w^11*z^27+195*t^14
*w^2*z^34+28*t^13*w^9*z^28+111*t^14*z^35+154*t^13*w^7*z^29+364*t^13*w^5*z^30-20
*t^12*w^12*z^24+258*t^13*w^3*z^31-147*t^12*w^10*z^25+56*t^13*w*z^32-427*t^12*w^
8*z^26-778*t^12*w^6*z^27-40*t^11*w^13*z^21-1224*t^12*w^4*z^28-260*t^11*w^11*z^
22-865*t^12*w^2*z^29-730*t^11*w^9*z^23-359*t^12*z^30-1200*t^11*w^7*z^24-26*t^10
*w^14*z^18-1402*t^11*w^5*z^25-109*t^10*w^12*z^19-884*t^11*w^3*z^26+35*t^10*w^10
*z^20-184*t^11*w*z^27+1055*t^10*w^8*z^21-6*t^9*w^15*z^15+2620*t^10*w^6*z^22+36*
t^9*w^13*z^16+3043*t^10*w^4*z^23+480*t^9*w^11*z^17+2005*t^10*w^2*z^24+1708*t^9*
w^9*z^18+632*t^10*z^25+2840*t^9*w^7*z^19+26*t^8*w^14*z^13+2528*t^9*w^5*z^20+133
*t^8*w^12*z^14+1550*t^9*w^3*z^21+37*t^8*w^10*z^15+312*t^9*w*z^22-1155*t^8*w^8*z
^16-3228*t^8*w^6*z^17-40*t^7*w^13*z^11-3927*t^8*w^4*z^18-324*t^7*w^11*z^12-2429
*t^8*w^2*z^19-1118*t^7*w^9*z^13-632*t^8*z^20-2108*t^7*w^7*z^14-2314*t^7*w^5*z^
15+20*t^6*w^12*z^9-1328*t^7*w^3*z^16+163*t^6*w^10*z^10-300*t^7*w*z^17+679*t^6*w
^8*z^11+1746*t^6*w^6*z^12+2476*t^6*w^4*z^13+10*t^5*w^11*z^7+1513*t^6*w^2*z^14+
124*t^5*w^9*z^8+359*t^6*z^15+466*t^5*w^7*z^9+720*t^5*w^5*z^10+514*t^5*w^3*z^11-\
14*t^4*w^10*z^5+160*t^5*w*z^12-133*t^4*w^8*z^6-437*t^4*w^6*z^7-648*t^4*w^4*z^8-\
451*t^4*w^2*z^9+4*t^3*w^9*z^3-111*t^4*z^10+20*t^3*w^7*z^4+2*t^3*w^5*z^5-72*t^3*
w^3*z^6-42*t^3*w*z^7+7*t^2*w^6*z^2+38*t^2*w^4*z^3+51*t^2*w^2*z^4+17*t^2*z^5+2*t
*w^3*z+4*t*w*z^2-1)/(t^20*z^50+t^19*w^5*z^45-2*t^19*w^3*z^46+3*t^19*w*z^47-3*t^
18*w^8*z^41-t^18*w^6*z^42-18*t^18*w^4*z^43+3*t^17*w^11*z^37-18*t^18*w^2*z^44-2*
t^17*w^9*z^38-25*t^18*z^45-t^16*w^14*z^33-45*t^17*w^5*z^40+10*t^16*w^12*z^34-65
*t^17*w^3*z^41+39*t^16*w^10*z^35-31*t^17*w*z^42+142*t^16*w^8*z^36-6*t^15*w^15*z
^30+229*t^16*w^6*z^37-9*t^15*w^13*z^31+265*t^16*w^4*z^38+21*t^15*w^11*z^32+381*
t^16*w^2*z^39+235*t^15*w^9*z^33-15*t^14*w^16*z^27+216*t^16*z^40+610*t^15*w^7*z^
34-84*t^14*w^14*z^28+809*t^15*w^5*z^35-278*t^14*w^12*z^29+777*t^15*w^3*z^36-584
*t^14*w^10*z^30-20*t^13*w^17*z^24+101*t^15*w*z^37-1110*t^14*w^8*z^31-145*t^13*w
^15*z^25-2141*t^14*w^6*z^32-553*t^13*w^13*z^26-2838*t^14*w^4*z^33-1365*t^13*w^
11*z^27-15*t^12*w^18*z^21-2421*t^14*w^2*z^34-2751*t^13*w^9*z^28-102*t^12*w^16*z
^22-895*t^14*z^35-4727*t^13*w^7*z^29-252*t^12*w^14*z^23-5329*t^13*w^5*z^30+84*t
^12*w^12*z^24-6*t^11*w^19*z^18-2609*t^13*w^3*z^31+1910*t^12*w^10*z^25-15*t^11*w
^17*z^19-105*t^13*w*z^32+5588*t^12*w^8*z^26+219*t^11*w^15*z^20+10196*t^12*w^6*z
^27+1756*t^11*w^13*z^21-t^10*w^20*z^15+12124*t^12*w^4*z^28+5921*t^11*w^11*z^22+
16*t^10*w^18*z^16+6997*t^12*w^2*z^29+11744*t^11*w^9*z^23+221*t^10*w^16*z^17+
2023*t^12*z^30+15003*t^11*w^7*z^24+953*t^10*w^14*z^18+11581*t^11*w^5*z^25+1322*
t^10*w^12*z^19+6*t^9*w^19*z^13+3931*t^11*w^3*z^26-2629*t^10*w^10*z^20+27*t^9*w^
17*z^14+24*t^11*w*z^27-13637*t^10*w^8*z^21-215*t^9*w^15*z^15-24259*t^10*w^6*z^
22-2280*t^9*w^13*z^16-22171*t^10*w^4*z^23-8597*t^9*w^11*z^17-15*t^8*w^18*z^11-\
10798*t^10*w^2*z^24-17024*t^9*w^9*z^18-162*t^8*w^16*z^12-2640*t^10*z^25-18451*t
^9*w^7*z^19-600*t^8*w^14*z^13-10465*t^9*w^5*z^20-236*t^8*w^12*z^14-3387*t^9*w^3
*z^21+4742*t^8*w^10*z^15+20*t^7*w^17*z^9-96*t^9*w*z^22+15824*t^8*w^8*z^16+265*t
^7*w^15*z^10+23972*t^8*w^6*z^17+1453*t^7*w^13*z^11+19900*t^8*w^4*z^18+4209*t^7*
w^11*z^12+9217*t^8*w^2*z^19+6979*t^7*w^9*z^13-15*t^6*w^16*z^7+2023*t^8*z^20+
6939*t^7*w^7*z^14-204*t^6*w^14*z^8+4541*t^7*w^5*z^15-1150*t^6*w^12*z^9+1905*t^7
*w^3*z^16-3548*t^6*w^10*z^10+253*t^7*w*z^17-6914*t^6*w^8*z^11+6*t^5*w^15*z^5-\
9425*t^6*w^6*z^12+69*t^5*w^13*z^6-8710*t^6*w^4*z^13+279*t^5*w^11*z^7-4245*t^6*w
^2*z^14+409*t^5*w^9*z^8-895*t^6*z^15-78*t^5*w^7*z^9-t^4*w^14*z^3-717*t^5*w^5*z^
10-2*t^4*w^12*z^4-657*t^5*w^3*z^11+75*t^4*w^10*z^5-209*t^5*w*z^12+514*t^4*w^8*z
^6+1337*t^4*w^6*z^7+1649*t^4*w^4*z^8-3*t^3*w^11*z^2+969*t^4*w^2*z^9-30*t^3*w^9*
z^3+216*t^4*z^10-88*t^3*w^7*z^4-43*t^3*w^5*z^5+113*t^3*w^3*z^6+67*t^3*w*z^7-3*t
^2*w^8*z-29*t^2*w^6*z^2-86*t^2*w^4*z^3-86*t^2*w^2*z^4-25*t^2*z^5-t*w^5-6*t*w^3*
z-7*t*w*z^2+1), -(t^34*z^102-3*t^33*w^4*z^97-4*t^33*z^99-9*t^32*w^8*z^92+6*t^32
*w^6*z^93-31*t^32*w^4*z^94-5*t^31*w^12*z^87-43*t^32*w^2*z^95+36*t^31*w^10*z^88-\
32*t^32*z^96+31*t^31*w^8*z^89+180*t^31*w^6*z^90+30*t^30*w^14*z^83+286*t^31*w^4*
z^91+16*t^30*w^12*z^84+197*t^31*w^2*z^92+200*t^30*w^10*z^85+136*t^31*z^93+355*t
^30*w^8*z^86-60*t^29*w^16*z^79+589*t^30*w^6*z^87-184*t^29*w^14*z^80+858*t^30*w^
4*z^88-1090*t^29*w^12*z^81+1067*t^30*w^2*z^89-3380*t^29*w^10*z^82+406*t^30*z^90
-7046*t^29*w^8*z^83-214*t^28*w^16*z^76-8716*t^29*w^6*z^84-248*t^28*w^14*z^77-\
8500*t^29*w^4*z^85+674*t^28*w^12*z^78+210*t^27*w^20*z^71-5329*t^29*w^2*z^86+
2174*t^28*w^10*z^79+1904*t^27*w^18*z^72-1968*t^29*z^87-2310*t^28*w^8*z^80+7806*
t^27*w^16*z^73-7357*t^28*w^6*z^81+22116*t^27*w^14*z^74-420*t^26*w^22*z^67-9861*
t^28*w^4*z^82+51170*t^27*w^12*z^75-3556*t^26*w^20*z^68-9710*t^28*w^2*z^83+
102470*t^27*w^10*z^76-14266*t^26*w^18*z^69-2525*t^28*z^84+150382*t^27*w^8*z^77-\
36013*t^26*w^16*z^70+420*t^25*w^24*z^63+158060*t^27*w^6*z^78-66266*t^26*w^14*z^
71+2984*t^25*w^22*z^64+115973*t^27*w^4*z^79-98162*t^26*w^12*z^72+7285*t^25*w^20
*z^65+60530*t^27*w^2*z^80-109012*t^26*w^10*z^73-7160*t^25*w^18*z^66-240*t^24*w^
26*z^59+15896*t^27*z^81-80770*t^26*w^8*z^74-101310*t^25*w^16*z^67-841*t^24*w^24
*z^60-8405*t^26*w^6*z^75-352554*t^25*w^14*z^68+6000*t^24*w^22*z^61+37543*t^26*w
^4*z^76-771052*t^25*w^12*z^69+61509*t^24*w^20*z^62+75*t^23*w^28*z^55+34511*t^26
*w^2*z^77-1233550*t^25*w^10*z^70+256319*t^24*w^18*z^63-436*t^23*w^26*z^56+6776*
t^26*z^78-1531300*t^25*w^8*z^71+657875*t^24*w^16*z^64-8869*t^23*w^24*z^57-\
1392464*t^25*w^6*z^72+1164960*t^24*w^14*z^65-48384*t^23*w^22*z^58-10*t^22*w^30*
z^51-873096*t^25*w^4*z^73+1535046*t^24*w^12*z^66-122897*t^23*w^20*z^59+380*t^22
*w^28*z^52-377572*t^25*w^2*z^74+1598531*t^24*w^10*z^67-63063*t^23*w^18*z^60+
3138*t^22*w^26*z^53-79330*t^25*z^75+1264385*t^24*w^8*z^68+623746*t^23*w^16*z^61
+2825*t^22*w^24*z^54+681329*t^24*w^6*z^69+2467870*t^23*w^14*z^62-75479*t^22*w^
22*z^55-80*t^21*w^30*z^48+182528*t^24*w^4*z^70+5312120*t^23*w^12*z^63-521014*t^
22*w^20*z^56+280*t^21*w^28*z^49+29436*t^24*w^2*z^71+7928789*t^23*w^10*z^64-\
1892095*t^22*w^18*z^57+9908*t^21*w^26*z^50+6972*t^24*z^72+8636031*t^23*w^8*z^65
-4605912*t^22*w^16*z^58+70778*t^21*w^24*z^51+6915613*t^23*w^6*z^66-8137433*t^22
*w^14*z^59+267106*t^21*w^22*z^52-270*t^20*w^30*z^45+3897672*t^23*w^4*z^67-\
10765441*t^22*w^12*z^60+556371*t^21*w^20*z^53-2120*t^20*w^28*z^46+1428552*t^23*
w^2*z^68-10901014*t^22*w^10*z^61+249601*t^21*w^18*z^54-1586*t^20*w^26*z^47+
256782*t^23*z^69-8735052*t^22*w^8*z^62-2478206*t^21*w^16*z^55+64769*t^20*w^24*z
^48-5411154*t^22*w^6*z^63-9644166*t^21*w^14*z^56+519479*t^20*w^22*z^49-480*t^19
*w^30*z^42-2315755*t^22*w^4*z^64-20212525*t^21*w^12*z^57+2337876*t^20*w^20*z^50
-6265*t^19*w^28*z^43-694331*t^22*w^2*z^65-28432563*t^21*w^10*z^58+7335535*t^20*
w^18*z^51-41224*t^19*w^26*z^44-112080*t^22*z^66-28536123*t^21*w^8*z^59+16991757
*t^20*w^16*z^52-177421*t^19*w^24*z^45-20650337*t^21*w^6*z^60+29686577*t^20*w^14
*z^53-506640*t^19*w^22*z^46-420*t^18*w^30*z^39-10591279*t^21*w^4*z^61+39851459*
t^20*w^12*z^54-750813*t^19*w^20*z^47-6200*t^18*w^28*z^40-3429222*t^21*w^2*z^62+
41756404*t^20*w^10*z^55+782387*t^19*w^18*z^48-52524*t^18*w^26*z^41-551964*t^21*
z^63+34179438*t^20*w^8*z^56+7760534*t^19*w^16*z^49-330277*t^18*w^24*z^42+
21408308*t^20*w^6*z^57+23183242*t^19*w^14*z^50-1578432*t^18*w^22*z^43+9579383*t
^20*w^4*z^58+43901515*t^19*w^12*z^51-5717385*t^18*w^20*z^44+1400*t^17*w^28*z^37
+2722210*t^20*w^2*z^59+58473361*t^19*w^10*z^52-15938821*t^18*w^18*z^45+13720*t^
17*w^26*z^38+390098*t^20*z^60+55515336*t^19*w^8*z^53-35060029*t^18*w^16*z^46+
37076*t^17*w^24*z^39+37365655*t^19*w^6*z^54-61750572*t^18*w^14*z^47-121420*t^17
*w^22*z^40+420*t^16*w^30*z^33+17618709*t^19*w^4*z^55-86595130*t^18*w^12*z^48-\
1315020*t^17*w^20*z^41+9280*t^16*w^28*z^34+5284981*t^19*w^2*z^56-94984535*t^18*
w^10*z^49-5526330*t^17*w^18*z^42+87484*t^16*w^26*z^35+791368*t^19*z^57-79654653
*t^18*w^8*z^50-15730668*t^17*w^16*z^43+497349*t^16*w^24*z^36-49301752*t^18*w^6*
z^51-33834528*t^17*w^14*z^44+2031660*t^16*w^22*z^37+480*t^15*w^30*z^30-21439711
*t^18*w^4*z^52-55508292*t^17*w^12*z^45+6659353*t^16*w^20*z^38+9305*t^15*w^28*z^
31-5749764*t^18*w^2*z^53-68189138*t^17*w^10*z^46+18316491*t^16*w^18*z^39+78444*
t^15*w^26*z^32-759282*t^18*z^54-61296184*t^17*w^8*z^47+41890001*t^16*w^16*z^40+
393825*t^15*w^24*z^33-39246766*t^17*w^6*z^48+77756660*t^16*w^14*z^41+1385336*t^
15*w^22*z^34+270*t^14*w^30*z^27-17422960*t^17*w^4*z^49+114774026*t^16*w^12*z^42
+3856429*t^15*w^20*z^35+4060*t^14*w^28*z^28-5095010*t^17*w^2*z^50+131424941*t^
16*w^10*z^43+9103311*t^15*w^18*z^36+20506*t^14*w^26*z^29-741912*t^17*z^51+
112602589*t^16*w^8*z^44+18196666*t^15*w^16*z^37+1155*t^14*w^24*z^30+69196836*t^
16*w^6*z^45+29833990*t^15*w^14*z^38-499829*t^14*w^22*z^31+80*t^13*w^30*z^24+
29073171*t^16*w^4*z^46+39247013*t^15*w^12*z^39-3073060*t^14*w^20*z^32+400*t^13*
w^28*z^25+7543524*t^16*w^2*z^47+40779621*t^15*w^10*z^40-11152997*t^14*w^18*z^33
-8908*t^13*w^26*z^26+941714*t^16*z^48+32814096*t^15*w^8*z^41-29401425*t^14*w^16
*z^34-122158*t^13*w^24*z^27+20041995*t^15*w^6*z^42-59613115*t^14*w^14*z^35-\
724278*t^13*w^22*z^28+10*t^12*w^30*z^21+8844251*t^15*w^4*z^43-93277559*t^14*w^
12*z^36-2617399*t^13*w^20*z^29-280*t^12*w^28*z^22+2735557*t^15*w^2*z^44-\
110626824*t^14*w^10*z^37-6468775*t^13*w^18*z^30-5818*t^12*w^26*z^23+423272*t^15
*z^45-96653962*t^14*w^8*z^38-11512254*t^13*w^16*z^31-36029*t^12*w^24*z^24-\
59826504*t^14*w^6*z^39-14770698*t^13*w^14*z^32-56671*t^12*w^22*z^25-24866103*t^
14*w^4*z^40-12808971*t^13*w^12*z^33+451946*t^12*w^20*z^26-75*t^11*w^28*z^19-\
6397918*t^14*w^2*z^41-6245919*t^13*w^10*z^34+3329109*t^12*w^18*z^27-304*t^11*w^
26*z^20-774090*t^14*z^42-713657*t^13*w^8*z^35+11676092*t^12*w^16*z^28+8881*t^11
*w^24*z^21+451003*t^13*w^6*z^36+27034151*t^12*w^14*z^29+110344*t^11*w^22*z^22-\
224285*t^13*w^4*z^37+44995177*t^12*w^12*z^30+594449*t^11*w^20*z^23-359670*t^13*
w^2*z^38+55084242*t^12*w^10*z^31+1840373*t^11*w^18*z^24+240*t^10*w^26*z^17-\
108168*t^13*z^39+49309544*t^12*w^8*z^32+3418202*t^11*w^16*z^25+3161*t^10*w^24*z
^18+31438358*t^12*w^6*z^33+3296462*t^11*w^14*z^26+11416*t^10*w^22*z^19+13451971
*t^12*w^4*z^34-395488*t^11*w^12*z^27-41745*t^10*w^20*z^20+3515685*t^12*w^2*z^35
-5888263*t^11*w^10*z^28-552661*t^10*w^18*z^21+424744*t^12*z^36-8430207*t^11*w^8
*z^29-2472159*t^10*w^16*z^22-420*t^9*w^24*z^15-6205243*t^11*w^6*z^30-6578068*t^
10*w^14*z^23-6904*t^9*w^22*z^16-2497980*t^11*w^4*z^31-11790238*t^10*w^12*z^24-\
47285*t^9*w^20*z^17-508932*t^11*w^2*z^32-15066505*t^10*w^10*z^25-168080*t^9*w^
18*z^18-26454*t^11*z^33-14117465*t^10*w^8*z^26-283582*t^9*w^16*z^19-9601487*t^
10*w^6*z^27+66786*t^9*w^14*z^20+420*t^8*w^22*z^13-4437744*t^10*w^4*z^28+1402928
*t^9*w^12*z^21+7196*t^8*w^20*z^14-1219984*t^10*w^2*z^29+3260910*t^9*w^10*z^22+
54250*t^8*w^18*z^15-153260*t^10*z^30+4136976*t^9*w^8*z^23+237989*t^8*w^16*z^16+
3227448*t^9*w^6*z^24+677938*t^8*w^14*z^17+1487736*t^9*w^4*z^25+1331814*t^8*w^12
*z^18-210*t^7*w^20*z^11+359088*t^9*w^2*z^26+1882340*t^8*w^10*z^19-3304*t^7*w^18
*z^12+31970*t^9*z^27+1966598*t^8*w^8*z^20-23318*t^7*w^16*z^13+1509743*t^8*w^6*z
^21-99000*t^7*w^14*z^14+801157*t^8*w^4*z^22-283226*t^7*w^12*z^15+250203*t^8*w^2
*z^23-568326*t^7*w^10*z^16+34832*t^8*z^24-789870*t^7*w^8*z^17-346*t^6*w^16*z^10
-718980*t^7*w^6*z^18-4184*t^6*w^14*z^11-391253*t^7*w^4*z^19-20406*t^6*w^12*z^12
-108622*t^7*w^2*z^20-52110*t^6*w^10*z^13-11612*t^7*z^21-77490*t^6*w^8*z^14+60*t
^5*w^16*z^7-75889*t^6*w^6*z^15+1064*t^5*w^14*z^8-55543*t^6*w^4*z^16+7354*t^5*w^
12*z^9-25422*t^6*w^2*z^17+26296*t^5*w^10*z^10-4499*t^6*z^18+53778*t^5*w^8*z^11+
65168*t^5*w^6*z^12-30*t^4*w^14*z^5+45980*t^5*w^4*z^13-396*t^4*w^12*z^6+16575*t^
5*w^2*z^14-1920*t^4*w^10*z^7+2172*t^5*z^15-4295*t^4*w^8*z^8-4447*t^4*w^6*z^9-\
1618*t^4*w^4*z^10+5*t^3*w^12*z^3+415*t^4*w^2*z^11+24*t^3*w^10*z^4+234*t^4*z^12-\
115*t^3*w^8*z^5-856*t^3*w^6*z^6-1638*t^3*w^4*z^7-1103*t^3*w^2*z^8-208*t^3*z^9+9
*t^2*w^8*z^2+66*t^2*w^6*z^3+139*t^2*w^4*z^4+89*t^2*w^2*z^5+8*t^2*z^6+3*t*w^4*z+
12*t*w^2*z^2+8*t*z^3-1)/(t^36*z^108-t^35*w^6*z^102+2*t^35*w^4*z^103-6*t^35*w^2*
z^104-3*t^34*w^10*z^97-3*t^35*z^105+t^34*w^8*z^98-24*t^34*w^6*z^99-3*t^33*w^14*
z^92-23*t^34*w^4*z^100+9*t^33*w^12*z^93-41*t^34*w^2*z^101+9*t^33*w^10*z^94-t^32
*w^18*z^87-48*t^34*z^102+121*t^33*w^8*z^95+18*t^32*w^16*z^88+261*t^33*w^6*z^96+
47*t^32*w^14*z^89+375*t^33*w^4*z^97+210*t^32*w^12*z^90+9*t^31*w^20*z^83+284*t^
33*w^2*z^98+333*t^32*w^10*z^91-16*t^31*w^18*z^84+168*t^33*z^99+397*t^32*w^8*z^
92-158*t^31*w^16*z^85+579*t^32*w^6*z^93-993*t^31*w^14*z^86-36*t^30*w^22*z^79+
1661*t^32*w^4*z^94-3056*t^31*w^12*z^87-150*t^30*w^20*z^80+1698*t^32*w^2*z^95-\
6596*t^31*w^10*z^88-426*t^30*w^18*z^81+794*t^32*z^96-10913*t^31*w^8*z^89-266*t^
30*w^16*z^82+84*t^29*w^24*z^75-15761*t^31*w^6*z^90+1147*t^30*w^14*z^83+606*t^29
*w^22*z^76-14314*t^31*w^4*z^91+3023*t^30*w^12*z^84+2741*t^29*w^20*z^77-7776*t^
31*w^2*z^92+1107*t^30*w^10*z^85+8512*t^29*w^18*z^78-126*t^28*w^26*z^71-3300*t^
31*z^93-5765*t^30*w^8*z^86+22490*t^29*w^16*z^79-1092*t^28*w^24*z^72-20414*t^30*
w^6*z^87+57174*t^29*w^14*z^80-5288*t^28*w^22*z^73-30257*t^30*w^4*z^88+132583*t^
29*w^12*z^81-16720*t^28*w^20*z^74+126*t^27*w^28*z^67-23711*t^30*w^2*z^89+244070
*t^29*w^10*z^82-40136*t^28*w^18*z^75+1092*t^27*w^26*z^68-6048*t^30*z^90+335671*
t^29*w^8*z^83-84399*t^28*w^16*z^76+4350*t^27*w^24*z^69+326016*t^29*w^6*z^84-\
159136*t^28*w^14*z^77+6007*t^27*w^22*z^70-84*t^26*w^30*z^63+229506*t^29*w^4*z^
85-229223*t^28*w^12*z^78-23059*t^27*w^20*z^71-564*t^26*w^28*z^64+111392*t^29*w^
2*z^86-185885*t^28*w^10*z^79-161004*t^27*w^18*z^72-37*t^26*w^26*z^65+33023*t^29
*z^87+22430*t^28*w^8*z^80-541166*t^27*w^16*z^73+18845*t^26*w^24*z^66+36*t^25*w^
32*z^59+206888*t^28*w^6*z^81-1324760*t^27*w^14*z^74+126726*t^26*w^22*z^67+45*t^
25*w^30*z^60+234636*t^28*w^4*z^82-2553967*t^27*w^12*z^75+474909*t^26*w^20*z^68-\
2819*t^25*w^28*z^61+143410*t^28*w^2*z^83-3813048*t^27*w^10*z^76+1225703*t^26*w^
18*z^69-26513*t^25*w^26*z^62-9*t^24*w^34*z^55+20021*t^28*z^84-4232750*t^27*w^8*
z^77+2352346*t^26*w^16*z^70-117461*t^25*w^24*z^63+114*t^24*w^32*z^56-3458309*t^
27*w^6*z^78+3432661*t^26*w^14*z^71-275108*t^25*w^22*z^64+2097*t^24*w^30*z^57-\
2063618*t^27*w^4*z^79+3697857*t^26*w^12*z^72-143967*t^25*w^20*z^65+11062*t^24*w
^28*z^58+t^23*w^36*z^51-883350*t^27*w^2*z^80+2694571*t^26*w^10*z^73+1395868*t^
25*w^18*z^66+6604*t^24*w^26*z^59-60*t^23*w^34*z^52-194412*t^27*z^81+1012491*t^
26*w^8*z^74+6001405*t^25*w^16*z^67-226391*t^24*w^24*z^60-456*t^23*w^32*z^53-\
303265*t^26*w^6*z^75+14290038*t^25*w^14*z^68-1476209*t^24*w^22*z^61+2131*t^23*w
^30*z^54-617173*t^26*w^4*z^76+23878737*t^25*w^12*z^69-5264634*t^24*w^20*z^62+
45965*t^23*w^28*z^55+10*t^22*w^36*z^48-289646*t^26*w^2*z^77+30004484*t^25*w^10*
z^70-12752169*t^24*w^18*z^63+296436*t^23*w^26*z^56-102*t^22*w^34*z^49+9824*t^26
*z^78+29039527*t^25*w^8*z^71-22375271*t^24*w^16*z^64+1054265*t^23*w^24*z^57-\
2692*t^22*w^32*z^50+21217854*t^25*w^6*z^72-29284850*t^24*w^14*z^65+2060625*t^23
*w^22*z^58-18579*t^22*w^30*z^51+11277982*t^25*w^4*z^73-29422577*t^24*w^12*z^66+
558191*t^23*w^20*z^59-39453*t^22*w^28*z^52+45*t^21*w^36*z^45+4160854*t^25*w^2*z
^74-23383864*t^24*w^10*z^67-10152931*t^23*w^18*z^60+231940*t^22*w^26*z^53+288*t
^21*w^34*z^46+723680*t^25*z^75-14685395*t^24*w^8*z^68-37085051*t^23*w^16*z^61+
2260185*t^22*w^24*z^54-3020*t^21*w^32*z^47-6986027*t^24*w^6*z^69-78330452*t^23*
w^14*z^62+9890773*t^22*w^22*z^55-52069*t^21*w^30*z^48-2437087*t^24*w^4*z^70-\
118255803*t^23*w^12*z^63+28775761*t^22*w^20*z^56-345440*t^21*w^28*z^49+120*t^20
*w^36*z^42-1010546*t^24*w^2*z^71-135585322*t^23*w^10*z^64+61568761*t^22*w^18*z^
57-1333629*t^21*w^26*z^50+1737*t^20*w^34*z^43-335436*t^24*z^72-119384004*t^23*w
^8*z^65+101644578*t^22*w^16*z^58-3087228*t^21*w^24*z^51+8636*t^20*w^32*z^44-\
79543367*t^23*w^6*z^66+133371575*t^22*w^14*z^59-2826897*t^21*w^22*z^52-8975*t^
20*w^30*z^45-38409024*t^23*w^4*z^67+141650850*t^22*w^12*z^60+8751206*t^21*w^20*
z^53-364676*t^20*w^28*z^46+210*t^19*w^36*z^39-12193656*t^23*w^2*z^68+123635080*
t^22*w^10*z^61+47692357*t^21*w^18*z^54-2468634*t^20*w^26*z^47+3960*t^19*w^34*z^
40-1771516*t^23*z^69+89658539*t^22*w^8*z^62+127466721*t^21*w^16*z^55-10307250*t
^20*w^24*z^48+33796*t^19*w^32*z^41+52430694*t^22*w^6*z^63+239658681*t^21*w^14*z
^56-31760245*t^20*w^22*z^49+166665*t^19*w^30*z^42+23146042*t^22*w^4*z^64+
339567623*t^21*w^12*z^57-77911627*t^20*w^20*z^50+492510*t^19*w^28*z^43+252*t^18
*w^36*z^36+7722658*t^22*w^2*z^65+367406801*t^21*w^10*z^58-157581774*t^20*w^18*z
^51+669466*t^19*w^26*z^44+5292*t^18*w^34*z^37+1362430*t^22*z^66+302169976*t^21*
w^8*z^59-264854884*t^20*w^16*z^52-1349348*t^19*w^24*z^45+52208*t^18*w^32*z^38+
185485301*t^21*w^6*z^60-368197833*t^20*w^14*z^53-11708193*t^19*w^22*z^46+323228
*t^18*w^30*z^39+81755204*t^21*w^4*z^61-421714181*t^20*w^12*z^54-44135557*t^19*w
^20*z^47+1445010*t^18*w^28*z^40+210*t^17*w^36*z^33+22711950*t^21*w^2*z^62-\
396207778*t^20*w^10*z^55-120862268*t^19*w^18*z^48+5189844*t^18*w^26*z^41+4464*t
^17*w^34*z^34+2898554*t^21*z^63-300432012*t^20*w^8*z^56-259018461*t^19*w^16*z^
49+16227740*t^18*w^24*z^42+43948*t^17*w^32*z^35-177219840*t^20*w^6*z^57-\
438362897*t^19*w^14*z^50+45508458*t^18*w^22*z^43+266837*t^17*w^30*z^36-76425268
*t^20*w^4*z^58-582550428*t^19*w^12*z^51+112486410*t^18*w^20*z^44+1148292*t^17*w
^28*z^37+120*t^16*w^36*z^30-21849264*t^20*w^2*z^59-598680250*t^19*w^10*z^52+
238516436*t^18*w^18*z^45+3929700*t^17*w^26*z^38+2313*t^16*w^34*z^31-3033924*t^
20*z^60-464214848*t^19*w^8*z^53+427193732*t^18*w^16*z^46+11843130*t^17*w^24*z^
39+18626*t^16*w^32*z^32-264259094*t^19*w^6*z^54+639455064*t^18*w^14*z^47+
32893401*t^17*w^22*z^40+70605*t^16*w^30*z^33-106389641*t^19*w^4*z^55+786827772*
t^18*w^12*z^48+82189107*t^17*w^20*z^41+16096*t^16*w^28*z^34+45*t^15*w^36*z^27-\
26641826*t^19*w^2*z^56+776882916*t^18*w^10*z^49+177311542*t^17*w^18*z^42-\
1275484*t^16*w^26*z^35+612*t^15*w^34*z^28-3143968*t^19*z^57+596838548*t^18*w^8*
z^50+322276735*t^17*w^16*z^43-7970472*t^16*w^24*z^36+640*t^15*w^32*z^29+
342775504*t^18*w^6*z^51+484882525*t^17*w^14*z^44-30708649*t^16*w^22*z^37-45809*
t^15*w^30*z^30+139600414*t^18*w^4*z^52+589223474*t^17*w^12*z^45-90721561*t^16*w
^20*z^38-493592*t^15*w^28*z^31+10*t^14*w^36*z^24+35683278*t^18*w^2*z^53+
560213202*t^17*w^10*z^46-218994974*t^16*w^18*z^39-2829461*t^15*w^26*z^32-6*t^14
*w^34*z^25+4330900*t^18*z^54+402650580*t^17*w^8*z^47-435866700*t^16*w^16*z^40-\
11046052*t^15*w^24*z^33-2908*t^14*w^32*z^26+211283862*t^17*w^6*z^48-707611257*t
^16*w^14*z^41-32571185*t^15*w^22*z^34-37537*t^14*w^30*z^27+77650285*t^17*w^4*z^
49-920777389*t^16*w^12*z^42-77083682*t^15*w^20*z^35-227161*t^14*w^28*z^28+t^13*
w^36*z^21+18320150*t^17*w^2*z^50-936771324*t^16*w^10*z^43-150480235*t^15*w^18*z
^36-700318*t^14*w^26*z^29-48*t^13*w^34*z^22+2145708*t^17*z^51-720339076*t^16*w^
8*z^44-241656439*t^15*w^16*z^37-329487*t^14*w^24*z^30-996*t^13*w^32*z^23-\
402210356*t^16*w^6*z^45-312953567*t^15*w^14*z^38+7205905*t^14*w^22*z^31-5013*t^
13*w^30*z^24-155083962*t^16*w^4*z^46-318718819*t^15*w^12*z^39+38924791*t^14*w^
20*z^32+27345*t^13*w^28*z^25-36958030*t^16*w^2*z^47-248502359*t^15*w^10*z^40+
122591935*t^14*w^18*z^33+503706*t^13*w^26*z^26-9*t^12*w^34*z^19-4163732*t^16*z^
48-144915366*t^15*w^8*z^41+280182834*t^14*w^16*z^34+3269433*t^13*w^24*z^27-\
62832745*t^15*w^6*z^42+491014073*t^14*w^14*z^35+12979743*t^13*w^22*z^28+2641*t^
12*w^30*z^21-20129702*t^15*w^4*z^43+664542792*t^14*w^12*z^36+35608061*t^13*w^20
*z^29+34546*t^12*w^28*z^22-5077334*t^15*w^2*z^44+685133706*t^14*w^10*z^37+
70664285*t^13*w^18*z^30+203142*t^12*w^26*z^23-745654*t^15*z^45+524565059*t^14*w
^8*z^38+101758963*t^13*w^16*z^31+542375*t^12*w^24*z^24+36*t^11*w^32*z^17+
287879712*t^14*w^6*z^39+100523238*t^13*w^14*z^32-429297*t^12*w^22*z^25+525*t^11
*w^30*z^18+107548344*t^14*w^4*z^40+53603671*t^13*w^12*z^33-9453064*t^12*w^20*z^
26+607*t^11*w^28*z^19+24946198*t^14*w^2*z^41-10525496*t^13*w^10*z^34-40998909*t
^12*w^18*z^27-40937*t^11*w^26*z^20+2737322*t^14*z^42-43453684*t^13*w^8*z^35-\
108706281*t^12*w^16*z^28-431535*t^11*w^24*z^21-34354743*t^13*w^6*z^36-203586598
*t^12*w^14*z^29-2276228*t^11*w^22*z^22-84*t^10*w^30*z^15-13443994*t^13*w^4*z^37
-281508133*t^12*w^12*z^30-7415459*t^11*w^20*z^23-1740*t^10*w^28*z^16-2371306*t^
13*w^2*z^38-290186428*t^12*w^10*z^31-15419662*t^11*w^18*z^24-14005*t^10*w^26*z^
17-70564*t^13*z^39-221254609*t^12*w^8*z^32-18497443*t^11*w^16*z^25-42157*t^10*w
^24*z^18-121971339*t^12*w^6*z^33-4017086*t^11*w^14*z^26+128802*t^10*w^22*z^19-\
46149151*t^12*w^4*z^34+28956789*t^11*w^12*z^27+1755287*t^10*w^20*z^20+126*t^9*w
^28*z^13-10925802*t^12*w^2*z^35+58557244*t^11*w^10*z^28+8150001*t^10*w^18*z^21+
2940*t^9*w^26*z^14-1225424*t^12*z^36+60351097*t^11*w^8*z^29+23099098*t^10*w^16*
z^22+29550*t^9*w^24*z^15+37556808*t^11*w^6*z^30+44781755*t^10*w^14*z^23+164351*
t^9*w^22*z^16+13954630*t^11*w^4*z^31+62142013*t^10*w^12*z^24+527807*t^9*w^20*z^
17+2819746*t^11*w^2*z^32+63470229*t^10*w^10*z^25+811360*t^9*w^18*z^18-126*t^8*w
^26*z^11+208464*t^11*z^33+48785883*t^10*w^8*z^26-598728*t^9*w^16*z^19-3024*t^8*
w^24*z^12+28238197*t^10*w^6*z^27-6078320*t^9*w^14*z^20-31832*t^8*w^22*z^13+
11648751*t^10*w^4*z^28-15145013*t^9*w^12*z^21-193548*t^8*w^20*z^14+3003282*t^10
*w^2*z^29-22089108*t^9*w^10*z^22-755540*t^8*w^18*z^15+365084*t^10*z^30-21125050
*t^9*w^8*z^23-1996379*t^8*w^16*z^16+84*t^7*w^24*z^9-13420127*t^9*w^6*z^24-\
3687284*t^8*w^14*z^17+1950*t^7*w^22*z^10-5351700*t^9*w^4*z^25-4911343*t^8*w^12*
z^18+19697*t^7*w^20*z^11-1183400*t^9*w^2*z^26-4966381*t^8*w^10*z^19+114436*t^7*
w^18*z^12-103332*t^9*z^27-4100224*t^8*w^8*z^20+428230*t^7*w^16*z^13-2835964*t^8
*w^6*z^21+1101866*t^7*w^14*z^14-36*t^6*w^22*z^7-1490062*t^8*w^4*z^22+2043729*t^
7*w^12*z^15-750*t^6*w^20*z^8-477162*t^8*w^2*z^23+2825642*t^7*w^10*z^16-6522*t^6
*w^18*z^9-68735*t^8*z^24+2924085*t^7*w^8*z^17-30802*t^6*w^16*z^10+2158288*t^7*w
^6*z^18-86435*t^6*w^14*z^11+1023966*t^7*w^4*z^19-148721*t^6*w^12*z^12+9*t^5*w^
20*z^5+262832*t^7*w^2*z^20-156675*t^6*w^10*z^13+140*t^5*w^18*z^6+26761*t^7*z^21
-93623*t^6*w^8*z^14+670*t^5*w^16*z^7-7172*t^6*w^6*z^15-305*t^5*w^14*z^8+45671*t
^6*w^4*z^16-14332*t^5*w^12*z^9+34393*t^6*w^2*z^17-58530*t^5*w^10*z^10-t^4*w^18*
z^3+7268*t^6*z^18-118149*t^5*w^8*z^11-136527*t^5*w^6*z^12+191*t^4*w^14*z^5-\
91026*t^5*w^4*z^13+1762*t^4*w^12*z^6-31058*t^5*w^2*z^14+6895*t^4*w^10*z^7-3896*
t^5*z^15+13665*t^4*w^8*z^8+13779*t^4*w^6*z^9-3*t^3*w^14*z^2+6137*t^4*w^4*z^10-\
37*t^3*w^12*z^3+280*t^4*w^2*z^11-119*t^3*w^10*z^4-286*t^4*z^12+155*t^3*w^8*z^5+
1485*t^3*w^6*z^6+2645*t^3*w^4*z^7+1658*t^3*w^2*z^8-3*t^2*w^10*z+296*t^3*z^9-37*
t^2*w^8*z^2-152*t^2*w^6*z^3-241*t^2*w^4*z^4-131*t^2*w^2*z^5-12*t^2*z^6-t*w^6-8*
t*w^4*z-18*t*w^2*z^2-9*t*z^3+1)]:
RETURN(L[m]):
fi:

end:

#GFmd(m,t,z,w): inputs a positive integer m, and variables t,z,w, outputs the rational function of t,z,w
#whose coefficient of t^n*z^i*w^j is the number of tilings of an m by n rectangle with exactly i dimers and j monomers.
#Try:
#GFmd(2,t,z,w);
GFmd:=proc(m,t,z,w) local f,H:
option remember:
f:=GFtH(m,[Rt(1,2),Rt(2,1),Rt(1,1)],t,H);
subs({H[{[0, 0], [1, 0]}]=z,H[{[0, 0], [0, 1]}]=z, H[{[0, 0]}]=w},f):
end:


#Gleb(p): the entroy function for the plane up to 43 terms
Gleb:=proc(p)

1/2*(p*log(4) - p*log(p) - 2*(1-p)*log(1-p) - p)+
1/16*p^2 + 1/192*p^3 + 7/1536*p^4 + 41/10240*p^5 
+ 181/61440*p^6 + 757/344064*p^7 + 3291/1835008*p^8 + 14689/9437184*p^9 + 64771/47185920*p^10 
+ 276101/230686720*p^11 + 1132693/1107296256*p^12 + 4490513/5234491392*p^13 + 17337685/24427626496*p^14 
+ 65867621/112742891520*p^15 + 249437227/515396075520*p^16 
+ 955110593/2336462209024*p^17 + 3740591431/10514079940608*p^18 + 15039656569/47004122087424*p^19 + 61727254227/208907209277440*p^20 + 255640084561/923589767331840*p^21 + 50273131919/193514046488576*p^22 + 4312434281365/17803292276948992*p^23 + 5789230773063/25895697857380352*p^24 + 69044819053441/337769972052787200*p^25 + 272097812497681/1463669878895411200*p^26 + 1068966474984721/6323053876828176384*p^27 + 601281977474899/3891110078048108544*p^28 + 16672616519735441/117021532717594968064*p^29 + 66545602395606901/501520854503978434560*p^30 + 267471214350929957/2144433998568735375360*p^31 + 1080431496491179115/9149585060559937601536*p^32 + 4374403039126240385/38959523483674573012992*p^33 + 17705045340400677607/165577974805616935305216*p^34 + 71484177460946258777/702452014326859725537280*p^35 + 287529593953850293471/2975090884207876484628480*p^36 + 1151710503160001680385/12580384310364734849286144*p^37 + 4596336312298962012663/53117178199317769363652608*p^38 + 18298456303802689186745/223953508083610054614319104*p^39 + 72784234597284215364691/942962139299410756270817280*p^40 + 289698911730110389042529/3965276688335983693036257280*p^41 + 1155125274097244765650075/16654162091011131510752280576*p^42 + 4616317010648384103125561/69866240967168649264619323392*p^43:
end:


#MaxFun(L): Given a list of values [a,b], finds  its maximum
MaxFun:=proc(L) local i:
for i from 1 to nops(L) while L[i][2]-L[i][2]>0 do od:
(L[i-1][2]+L[i][2])/2:
end: