Table Of Content

POLYHEDRALOMEGA:ANEWALGORITHMFORSOLVINGLINEARDIOPHANTINESYSTEMS FELIXBREUERANDZAFEIRAKISZAFEIRAKOPOULOS ABSTRACT. PolyhedralOmegaisanewalgorithmforsolvinglinearDiophantinesystems(LDS),i.e.,for computingamultivariaterationalfunctionrepresentationofthesetofallnon-negativeintegersolutionsto asystemoflinearequationsandinequalities. PolyhedralOmegacombinesmethodsfrompartitionanal- ysiswithmethodsfrompolyhedralgeometry. Inparticular,wecombineMacMahon’siterativeapproach basedontheOmegaoperatorandexplicitformulasforitsevaluationwithgeometrictoolssuchasBrion decompositionsandBarvinok’sshortrationalfunctionrepresentations.Inthisway,weconnecttworecent branchesofresearchthathavesofarremainedseparate,unifiedbytheconceptofsymbolicconeswhich 5 weintroduce.TheresultingLDSsolverPolyhedralOmegaissignificantlyfasterthanprevioussolversbased 1 onpartitionanalysisanditiscompetitivewithstate-of-the-artLDSsolversbasedongeometricmethods. 0 Mostimportantly,thissynthesisofideasmakesPolyhedralOmegathesimplestalgorithmforsolvinglin- 2 earDiophantinesystemsavailabletodate.Moreover,weprovideanillustratedgeometricinterpretationof n partitionanalysis,withtheaimofmakingideasfrombothareasaccessibletoreadersfromawiderangeof a backgrounds. J 0 CONTENTS 3 1. Introduction 2 ] O 2. PartitionAnalysisanditsGeometricInterpretation 6 2.1. MacMahon 6 C 2.2. Elliott 12 . h 2.3. Andrews-Paule-Riese 13 at 2.4. QuestionsofConvergence 16 m 2.5. ApplyingOmegavs.SolvingLDS 17 2.6. Xin’sPartialFractionDecomposition 18 [ 3. PolyhedralOmega–theNewAlgorithm 19 1 3.1. MacMahonLifting 20 v 3.2. IterativeEliminationoftheLastCoordinateusingSymbolicCones 20 3 7 3.3. ConversiontoRationalFunctions 29 7 3.4. RationalFunctioninNormalForm 35 7 4. ComputationalComplexity 35 0 4.1. StructureTheorem 35 . 1 4.2. ComplexityofLiftingandElimination 38 0 4.3. ComplexityofRationalFunctionConversion 40 5 4.4. SummaryofComplexityAnalysis 42 1 5. ComparisonwithRelatedWork 43 : v 5.1. SpeedupviaSymbolicCones 43 i X 5.2. FirstPolynomial-TimeAlgorithmfromPartitionAnalysisFamily 45 5.3. ComparisonwithPolyhedralMethods 46 r a 6. SummaryandFutureResearch 46 Acknowledgements 47 References 48 Keywordsandphrases. LinearDiophantinesystem,linearinequalitysystem,integersolutions,partitionanalysis,partition theory,polyhedralgeometry,rationalfunction,symboliccone,generatingfunction,implementation,Omegaoperator. FelixBreuerwaspartiallysupportedbyGermanResearchFoundation(DFG)grantBR4251/1-1andbyAustrianScienceFund (FWF)specialresearchgroupAlgorithmicandEnumerativeCombinatoricsSFBF50-06. ZafeirakisZafeirakopouloswaspartiallysupportedbytheAustrianScienceFund(FWF)specialresearchgroupAlgorithmic andEnumerativeCombinatoricsSFBF50-06andW1214-N15projectDK06,theAustrianMarshallPlanFoundationandbytheEu- ropeanUnion(EuropeanSocialFund–ESF)andGreeknationalfundsthroughtheOperationalProgram"EducationandLifelong Learning"oftheNationalStrategicReferenceFramework(NSRF)-ResearchFundingProgram:THALIS–UOA(MIS375891). 1 2 FELIXBREUERANDZAFEIRAKISZAFEIRAKOPOULOS 1. INTRODUCTION InthisarticlewepresentPolyhedralOmega,anewalgorithmforcomputingamultivariaterational functionexpressionforthesetofallsolutionstoalinearDiophantinesystem.Thisalgorithmconnects twobranchesofresearch–partitionanalysisandpolyhedralgeometry–betweenwhichtherehasbeen littleinteractioninthepast. Tomakethispaperaccessibletoresearchersfromeitherfield,aswellas toreaderswithotherbackgrounds,wegiveanelementarypresentationofthealgorithmitselfandwe takecaretomotivatethekeyideasbehindit,inparticulartheirgeometriccontent.Tobegin,weusethis introductiontodefinetheproblemofcomputingrationalfunctionsolutionstolinearDiophantinesys- temsandtogiveanoverviewofthealgorithmsdevelopedinpartitionanalysisandpolyhedralgeometry tosolveit,beforepointingoutthebenefitsofPolyhedralOmega. LinearDiophantineSystemsandRationalFunctions. Let A∈(cid:90)m×n beanintegermatrixandletb∈ (cid:90)m be an integer vector. In this article, we are interested in finding the set of non-negative integer vectorsx∈(cid:90)n suchthat Ax(cid:202)b. Sincewearerestrictingourattentiontonon-negativesolutionsx∈ (cid:202)0 (cid:90)n ,itisequivalenttoconsidersystemsconsistingofanycombinationofequationsandinequalities. (cid:202)0 However,tostreamlinethisarticle,wearegoingtofocusonthecase Ax(cid:202)b. Wecallsuchasystemof constraints,givenbyAandb,alinearDiophantinesystem(LDS). LinearDiophantinesystemsareofgreatimportance,bothinpracticeandintheory. Forexample, theIntegerProgramming(IP)problem–whichisaboutcomputingasolutiontoanLDSthatmaximizes agivenlinearfunctional–playsapivotalroleinoperationsresearchandcombinatorialoptimization [55]. Inthisarticle,wearenotinterestedinfindingoneoptimalsolution,however. Instead,wewishto computearationalfunctionrepresentationofthesetofallsolutionstoanLDS. Of course, the set of solutions to a linear Diophantine system may be infinite. For example, the equation x −x =0, which is equivalent to the system x −x (cid:201)0 and x −x (cid:202)0, has the solution 1 2 1 2 1 2 set{(0,0),(1,1),(2,2),...}.Onewaytorepresentsuchinfinitesetsofsolutionsisviamultivariaterational functions.Weidentifyavectorx∈(cid:90)nwiththemonomialzx:=zx1·zx2·...·zxninnvariables.Then,using 1 2 n the geometric series expansion formula, we can represent the above set of solutions by the rational function ∞ 1 = 1 =(cid:88)z(i,i)=z(0,0)+z(1,1)+z(2,2)+.... 1−z(1,1) 1−z1z2 i=0 Ingeneral, werepresentasetS ⊂(cid:90)n ofnon-negativeintegervectorsbythemultivariategenerating (cid:202)0 function φ (z)=(cid:88)zx, S x∈S i.e.,thecoefficientzx inφ is1ifx∈S and0ifx(cid:54)∈S. Itisawell-knownfactthatwhenS isthesetof S solutionstoalinearDiophantinesystem, thenφ isalwaysarationalfunctionρ . Itisthisrational S S functionρ thatweseektocomputewhensolvingagivenlinearDiophantinesystem.Wearenotgoing S tobeinterestedinthenormalformofρ ,though,sincethenormalformmaybeunnecessarilylarge. S Takeforexamplethesystemx +x =100.Here,thesetofsolutionsis 1 2 S={(100,0),(99,1),...,(0,100)} whichcanberepresentedbytherationalfunction z(100,0)−z(−1,101) =z(100,0)+z(99,1)+...+z(0,100). 1−z(−1,1) Note that the polynomial on the right-hand side is the normal form of this rational function, which has101terms. (Infact,thenormalformofarationalfunctionrepresentingafinitesetwillalwaysbe apolynomial.) Therationalfunctionexpressionontheleft-handsideisnotinnormalformsincethe denominatordividesthenumerator,howeveritismuchshorter,havingonly4terms. Wearetherefore interestedincomputingmultivariaterationalfunctionexpressions,whicharenotuniquelydetermined, asopposedtotheuniquerationalfunctioninnormalform.Giventhisterminologyandnotation,wecan nowstatethecomputationalproblemofsolvinglinearDiophantinesystemsthatthisarticleisabout. POLYHEDRALOMEGA:SOLVINGLINEARDIOPHANTINESYSTEMS 3 Problem1.1. RationalFunctionSolutionofLinearDiophantineSystem(rfsLDS) Input: A∈(cid:90)m×n,b∈(cid:90)m Output:Anexpressionfortherationalfunctionρ∈Q(z ,...,z )representingthesetofallnon-negative 1 n integervectorsx∈(cid:90)n suchthatAx(cid:202)b. (cid:202)0 SuchrationalfunctionsolutionstoLDSareofgreatimportanceinmanyapplications.Forexample, theycanbeusedtoprovetheoremsinnumbertheoryandcombinatorics[4,5,50],computevolumes [16], count integer points in polyhedra [17], to maximize non-linear functions over lattice points in polyhedral [32], to compute Pareto optima in multi-criteria optimization [30], to integrate and sum functionsoverpolyhedra[12], tocomputeGröbnerbasesoftoricideals[29], toperformvariousop- erationsonrationalfunctionsandquasipolynomials[15],andtosampleobjectsfrompolyhedra[52], fromcombinatorialfamilies[35]andfromstatisticaldistributions[28]. Werecommendthetextbooks [13,20,31]foranintroduction. Note that, while above rfsLDS is stated in terms of pure inequality systems, this restriction is not essential. Themethodsandalgorithmwepresentinthisarticlecanbeeasilyextendedtohandlethe caseofmixedsystemsdirectly,asdoesourimplementation[23]. PartitionAnalysis. OneofthemajorlandmarksinthelonghistoryofProblem1.1isMacMahon’ssemi- nalworkCombinatoryAnalysis[50],publishedin1915,wherehepresentspartitionanalysisasageneral frameworkforsolvinglinearDiophantinesystemsintheabovesense,particularlyinthecontextofpar- titiontheory. Thegeneralapproachofpartitionanalysisistoemployasetofexplicitanalyticformulas formanipulatingrationalfunctionexpressionsthatcanbeiterativelyappliedtoobtainasolutiontoa givenlinearDiophantinesystem. WewillexaminepartitionanalysisindetailinSection2;fornow,we willconfineourselvestoaquickpreviewtoconveythegeneralflavor. ConsiderthelinearDiophantineinequality (1) 2x +3x −5x (cid:202) 4 1 2 3 Tosolvethissystem,MacMahonintroducestheOmegaoperator Ω(cid:202) definedonrationalfunctionsin (cid:75)(x ,x ,x )(λ)intermsoftheirseriesexpansionby 1 2 3 +∞ +∞ Ω(cid:202) (cid:88) csλs=(cid:88)cs s=−∞ s=0 wherecs ∈(cid:75)(x1,x2,x3)forall s. Inotherwords, Ω(cid:202) actsbytruncatingalltermsintheseriesexpan- sionwithnegativeλexponentanddroppingthevariableλ. UsingtheOmegaoperator,thegenerating functionφ ofthesetofsolutionsSof(1)canbewrittenas S λ−4 φS(z)=Ω(cid:202)(1−z λ2)(1−z λ3)(1−z λ−5). 1 2 3 Theexpressionontheright-handsideofthisequationisknownasthecrudegeneratingfunction. PartitionanalysisisaboutdevelopingcalculiofexplicitformulasforevaluatingtheOmegaoperator whenappliedtorationalfunctionsofacertainform.MacMahon’s1stRuleprovidesatypicalexample: 1 1 Ω(cid:202)(1−λx)(cid:161)1− y (cid:162)=(1−x)(1−xsy). λs Appliediteratively,suchformulasallowustosymbolicallytransformthecrudegeneratingfunctioninto arationalfunctionexpressionforφ ,asdesired. S Attheturnofthe20thcentury,Elliottwasthefirsttodevelopsuchacalculusthatyieldedanalgo- rithmicmethodforsolvingalinearDiophantineequation.Atthesametime,MacMahonwaspursuing theambitiousandvisionaryprojectofapplyingsuchmethodstopartitiontheory. Whilealgorithmic, Elliott’smethodwasnotpracticalforMacMahon’spurposes,asitwastoocomputationallyintensiveto carryoutbyhand. Therefore,MacMahondevelopedandcontinuallyextendedhisowncalculuswhich enabledhimtosolvemanypartitiontheoreticproblems, eventhoughitdidnotconstituteageneral algorithm. Despitemanysuccesses,MacMahondidnotachievethegoalhehadoriginallysethimself ofusinghispartitionanalysiscalculustoprovehisfamousformulaforplanepartitions[5].Itisperhaps forthisreasonthatpartitionanalysisfelloutofinterestformuchofthe20thcentury,beforeitstarted toattractrenewedattentioninrecentdecades,beginningwithStanley[56]. 4 FELIXBREUERANDZAFEIRAKISZAFEIRAKOPOULOS OMEGA MacMahon Andrews Paule Riese OMEGA2 Elliott Xin PolyhedralOmega Ell Ell2 CTEuclid FIGURE1. Anoverviewofthefieldofpartitionanalysis. It was Andrews who observed the computational potential of MacMahon’s partition analysis and waitedfortherightproblemtoapplyit. Inthelate1990s,theseminalpaperofBousquet–Mélouand Eriksson[21]onlecture-hallpartitionsappeared. AndrewssuggestedtoPauletoexplorethecapacity ofpartitionanalysiscombinedwithsymboliccomputation. Thiscollaborationgaverisetoanongo- ingseriesofover10papersandincludesOMEGA[4]andOMEGA2[5],twofullyalgorithmicversionsof partitionanalysispoweredbysymboliccomputation. Thislineofresearchhasalsoattractedrenewed interestinthisfield,includingforexample,asequenceofpapersbyXin[64,63,65]whodevelopsal- ternativepartitionanalysiscalculiandcorrespondingsoftwarepackages.Anillustrationoftheconnec- tionsbetweenthesedifferentsolutionstotheproblemisgiveninFigure1. PolyhedralGeometry. Independentlyfromtheworkonpartitionanalysis,thefieldofpolyhedralge- ometry has made major contributions to the solution of linear Diophantine systems. Integer linear programming isahugeandveryactivefieldofresearchconcernedwiththecomputationnotofara- tionalfunctionrepresentationofthesetofallnon-negativeintegervectorssatisfyingagivenLDS,but insteadwiththecomputationofasinglesuchvectorwhichisoptimalwithrespecttoalinearfunctional. Theintegerprogrammingproblemisdistinctfromtheoneweareinterestedininthisarticle.Nonethe- less,tworesultsonintegerprogrammingareimportanttomentionhere. Ontheonehand,deciding satisfiabilityoflinearDiophantinesystemsisNP-complete,asmanyimportantandpracticaldecision problemscanbemodeledeasilyusinglinearDiophantinesystems[54,55].Ontheotherhand,Lenstra hasshownin1983[47]thatifthenumberofvariablesisfixed, thenthesatisfiabilityofalinearDio- phantinesystemcanbedecidedand(ifitissatisfiable)anoptimalsolutioncanbefoundinpolynomial time. The problem of computing a rational function solution to a linear Diophantine system (rfsLDS), whichweareconcernedwithinthisarticle,isalsoatleastashardasanyprobleminNP,inthefollowing sense. ItiseasytoseethatNP-hardclassesofthesatisfiabilityproblemforLDScanbereducedtothe problemofcomputingtheuniquenormalformrationalfunctionsolutionoftheLDS.Interpolationar- gumentsthenshowthatmoreoveranyrationalfunctionexpressionofpolynomialsizewillreadilyyield thesolutiontotheNP-hardprobleminpolynomialtime. Thus,thequestionarises,whetherrfsLDSalsobecomespolynomialtimecomputableifthenumber ofvariablesisfixed.Lenstra’salgorithmdoesnothelphere.Thekeyobstacleisthataprioritheencoding lengthoftheoutputmaybeexponentialintheencodinglengthoftheinput.Aswehavealreadyseenin theexampleabove,thenormalformoftherationalfunctionsolutionofx+y=bhasb+1termswhile theLDSitselfhasencodinglengthO(logb).Whileinthisparticularcase,wecanfindarationalfunction expression, z(b,0)−z(−1,b+1),whoseencodinglengthisinO(logb),itisfarfromclearthatsuchshortrational 1−z(−1,1) functionexpressionsexistingeneral.Theydo,however. POLYHEDRALOMEGA:SOLVINGLINEARDIOPHANTINESYSTEMS 5 In1993Barvinokwasabletoshowinaseminalpaper[16]thatforeveryLDStheredoesexistaratio- nalfunctionexpressionforitssetofsolutionswhoseencodinglengthisboundedpolynomiallyinterms oftheencodinglengthoftheinput–providedthatthenumberofvariablesisfixed.Moreover,Barvinok gaveanalgorithmforcomputingsuchashortrationalfunctionrepresentation.Thekeyingredientthat makestheseshortdecompositionspossibleisBarvinok’smethodforcomputingashortsigneddecom- positionofasimplicialconeintounimodularsimplicialcones.TheseBarvinokdecompositionswillplay animportantroleinthisarticleaswell. Barvinok’salgorithmforsolvingrfsLDScombinesBarvinokdecompositionswithanumberofother sophisticatedtoolsfrompolyhedralgeometry,rangingfromthecomputationofverticesandedgedi- rectionsofpolyhedra[40,41]tothetriangulationofpolyhedralcones[34,53].Thiscombinationmakes Barvinok’salgorithmverycomplexandmuchmoreinvolvedthantherelativelystraightforwardcollec- tionofrulesbasedonexplicitformulasthatMacMahonhadenvisionedinthepartitionanalysisframe- work.AtestamentofthisdifficultyisthefactthateventhoughBarvinok’salgorithmquicklygrabbedthe attentionofscientificcommunity,ittook10yearsandasequenceofadditionalresearchpapers,before Barvinok’salgorithmwasfirstimplementedin2004byDeLoeraetal.[33]. The importance of Barvinok’s short rational function representations in this area cannot be over- stated,astheyhavegivenrisetoawholefamilyofalgorithmsandtheoreticalcomplexityresultsfora largerangeofproblems,manyofwhichwehavementionedabove.Ofparticularnoteinourcontextis thatBarvinokandWoodsprovedatheoremin[15]whichimpliesinparticularthatMacMahon’sOmega operatorcanbeevaluatedinpolynomialtimeinfixeddimension–notonlyoncrudegeneratingfunc- tionsastheyariseinpartitionanalysisbutonawiderclassofrationalfunctions. Toconcludeouroverviewofpriorliterature,webrieflymentionacoupleoffurtherrelatedpublica- tions[25,43,44],adetaileddiscussionofwhichisoutofscopeofthisarticle.[43,44]givetwodifferent algorithmsforsolvingrfsLDSviaanexplicitformula, phrasedintermsofasummationoverallsub- matricesoftheoriginalsystem. Thiscanbeviewedasabruteforceiterationoverallsimplicialcones formed by the hyperplane arrangement given by the system, together with a recursive algorithm for computingthenumeratorsofthecorrespondingrationalfunctions. Incontrasttoanexplicitformula givenin[25],theformulasfrom[43,44]arealgorithmicallytractable,thoughtherunningtimecanbe exponential, even if the dimension is fixed. Since the iterative application of Ω(cid:202) is not the focus of [43,44],wedonotviewthesealgorithmsaspartofthepartitionanalysisfamilyforthepurposesofthis article.Toourknowledgethesealgorithmshavenotbeenimplemented. OurContribution. Inthisarticle,webringpartitionanalysisandpolyhedralgeometrytogether.Some- whatsurprisingly,thesetwobranchesofresearchhavesofarhadrelativelylittleinteraction. Tobridge thisgap, weprovideadetailedandillustratedstudyofthesetwofieldsfromaunifiedpointofview. Inparticular,weinterpretthepartitionanalysiscalculiprovidedbyElliott,MacMahonandAndrews– Paule–Riesefromageometricperspective. Forexample, MacMahon’sansatzofsettingupthecrude generating function corresponds in the language of geometry to a clever way of writing an arbitrary polyhedronastheintersectionofasimplicialconewithacollectionofcoordinatehalf-spaces–acon- structionwhichwecalltheMacMahonlifting. Whiletheconnectionbetweenrationalfunctionsand polyhedralconesiswell-knownandadiscussionoftheElliott–MacMahonpartitionanalysisalgorithm fromthispointofviewcanbefoundin[56],someofthismaterialisnew,includingourgeometricdis- cussionofAndrews–Paule–RieseandXin. ThemaincontributionofthisarticleisPolyhedralOmega,anewalgorithmforsolvinglinearDio- phantinesystemswhichisasynthesisofkeyideasfrombothpartitionanalysisandpolyhedralgeom- etry. Justlikepartitionanalysisalgorithms,PolyhedralOmegaisaniterativealgorithmgiveninterms ofafewsimpleexplicitformulas. Thekeydifferenceisthattheformulasarenotgivenintermsofra- tionalfunctions,butinsteadintermsofsymboliccones. Thesearesymbolicexpressionsforsimplicial polyhedral cones, given interms of their generators. Applyingideas like Brion’s theorem from polyhedral geometry gives rise to an explicit calculus of rules for intersecting polyhedra with coordinate half-spaces,i.e.,forapplyingtheOmegaoperatoronthelevelofsymboliccones. After the Omega operator has been applied, the final symbolic cone expression can then be con- verted to a rational functionexpressionusing either an explicit formula based on the Smith Normal 6 FELIXBREUERANDZAFEIRAKISZAFEIRAKOPOULOS Formofamatrix,orusingBarvinokdecompositions.TheSmithNormalFormapproachhastheadvan- tageofbeingsimplerandeasiertoimplement,whiletheBarvinokdecompositionapproachismuch fasteronmanyclassesofLDSsandproducesshorterrationalfunctionexpressions.However,formany applicationsitcanbeadvantageoustonotconvertthesymbolicconeexpressiontoarationalfunction expressionatall,sincethisconversioncanincreasethesizeoftheoutputdramatically.Especiallywhen theoutputisintendedforinspectionbyhumans,westronglyrecommendworkingwiththesymbolic coneexpressionasitistypicallymuchmoreintelligible. PolyhedralOmegahasseveraladvantages. Incomparisonwithpartitionanalysismethods,Polyhe- dralOmegaismuchfaster. Infact,oncertainclassesofexamplesitoffersexponentialspeedupsover previouspartitionanalysisalgorithms–evenifBarvinokdecompositionsarenotused! –byvirtueof workingwithsymbolicconesinsteadofrationalfunctionsandthroughanimprovedchoiceoftransfor- mationrules.Incomparisonwithpolyhedralgeometrymethods,PolyhedralOmegaismuchsimpler.It doesnotrequiresophisticatedalgorithmsforexplicitlycomputingallverticesandedge-directionsofa polyhedronorfortriangulatingnon-simplicialcones. Thesecomputationshappenimplicitlythrough thestraightforwardapplicationofafewsimplerulesgivenbyexplicitformulas. Atthesametime, if Barvinokdecompositionsareusedfortheconversionofsymbolicconesintorationalfunctions,then PolyhedralOmegarunsinpolynomialtimeinfixeddimension,i.e.,itliesinthesamecomplexityclass asthefastest-knownalgorithmsfortherfsLDSproblem. PolyhedralOmegaisthefirstalgorithmfrom thepartitionanalysisfamilythatachievespolynomialruntimeinfixeddimension. Moreover,Polyhe- dralOmegaisthefirstpartitionanalysisalgorithmforwhichadetailedcomplexityanalysisisavailable. Thisanalysisismadepossiblebyourgeometricpointofview. On the whole, Polyhedral Omega is the overall simplest algorithm for solving linear Diophantine systems,whileatthesametimeitsperformanceiscompetitivewiththecurrentstate-of-the-art. Inthisarticle,wegiveatheoreticaldescriptionandanalysisofPolyhedralOmega,togetherwitha detailedmotivationandillustrationofthecentralideasthatmakestheexpositionaccessibletoreaders withoutpriorexperienceineitherpartitionanalysisorpolyhedralgeometry. Thedefinitionoftheal- gorithmisexplicitenoughtomakeitpossibleforinterestedreaderstoimplementthealgorithmeven withoutbeingexpertsinthefield. Moreover, wehaveimplementedPolyhedralOmegaontopofthe Sagesystem[57]andthecodeispubliclyavailable[23]. OrganizationofthisArticle. InSection2,wegiveageometricinterpretationofpreviouspartitionanal- ysiscalculibyElliott,MacMahon,Andrews-Paule-RieseandXin.InSection3,wepresentthePolyhedral Omegaalgorithmandproveitscorrectness. Wetakeparticularcaretomotivateandillustratethekey ideasinanaccessiblemanner. InSection4,wegiveadetailedcomplexityanalysisofouralgorithm. InSection5,wethencomparePolyhedralOmegatootheralgorithmsfrombothpartitionanalysisand polyhedralgeometry. Finally,weconcludeinSection6,wherewealsopointoutdirectionsforfuture research. 2. PARTITIONANALYSISANDITSGEOMETRICINTERPRETATION Inthissectionwepresentsomeofthemajorlandmarksofpartitionanalysisandinterprettheresults, whicharetypicallyphrasedinthelanguageofanalysis,fromtheperspectiveofpolyhedralgeometry. Inthefollowing,wewillintroducetherequiredconceptsandterminologyfrompolyhedralgeometry aswegoalong. However,acomprehensiveintroductiontopolyhedralgeometryisoutofscopeofthis article.Wethereforereferthereadertotheexcellenttextbooks[20,31,54,66]forfurtherdetails. 2.1. MacMahon. TheΩ(cid:202)operatorandthecrudegeneratingfunction. MacMahonpresentedhisinvestigationsconcern- ingintegerpartitiontheoryinaseriesof(seven)memoirs,publishedbetween1895and1916. Inthe secondmemoir[49],MacMahonobservesthatthetheoryofpartitionsofnumbersdevelopsinparallel tothatoflinearDiophantineequations.Inparticular,hedefinestheΩ(cid:202)operatorinArticle66: “SupposewehaveafunctionF(x,a)whichcanbeexpandedinascendingpowersofx. Suchex- pansionbeingeitherfiniteorinfinite,thecoefficientsofthevariouspowersofxarefunctionsofa whichingeneralinvolvebothpositiveandnegativepowersofa.Wemayrejectalltermscontain- ingnegativepowersofaandsubsequentlyputaequaltounity. Wethusarriveatafunctionofx POLYHEDRALOMEGA:SOLVINGLINEARDIOPHANTINESYSTEMS 7 only,whichmayberepresentedafterCayley(modifiedbytheassociationwiththesymbol(cid:202))by Ω(cid:202)F(x,a)thesymbol(cid:202)denotingthatthetermsretainedarethoseinwhichthepowerofais(cid:202)0.” Amodernwordingofthedefinitionisgivenin[4]. TheΩ(cid:202)operatorisdefinedonfunctionswithabso- lutelyconvergentmultisumexpansions ∞ ∞ ∞ (cid:88) (cid:88) ··· (cid:88) A λs1λs2···λsr s1,s2,...,sr 1 2 r s1=−∞s2=−∞ sr=−∞ inanopenneighborhoodofthecomplexcircles|λi|=1.TheactionofΩ(cid:202)isgivenby ∞ ∞ ∞ ∞ ∞ ∞ Ω(cid:202) (cid:88) (cid:88) ··· (cid:88) As1,s2,...,srλ1s1λ2s2···λrsr := (cid:88) (cid:88) ··· (cid:88) As1,s2,...,sr. s1=−∞s2=−∞ sr=−∞ s1=0s2=0 sr=0 MacMahonproceedswithdefiningwhathecallsthecrudegeneratingfunction.Thisisa(multivariate) generalizationofanideaappearinginCayleyaswellasinElliott’swork(seenextsection).Givenalinear systemofDiophantineinequalities,MacMahonintroducesanextravariableforeachinequality.These extravariablesaredenotedbyλ.GivenA∈(cid:90)m×d andb∈(cid:90)mweconsiderthegeneratingfunction m Φλ (z,λ)= (cid:88) zx(cid:89)λAix−bi A,b i x∈(cid:90)d(cid:202)0 i=1 andbasedonthegeometricseriesexpansionformula(1−z)−1=(cid:80)x(cid:202)0zx wecantransformtheseries intoarationalfunction.TherationalformofΦλ (z)isdenotedbyρλ (z),i.e., A,b A,b ρλ (z,λ)=λ−b(cid:89)m 1 . A,b i=1(cid:179)1−ziλ(AT)i(cid:180) Wecallρλ (z,λ)theλ-generatingfunctionof Aandb. Thecrudegeneratingfunctionisthesyntactic A,b expression   Ω λ−b(cid:89)m 1 (cid:202) (cid:179) (cid:180) i=1 1−ziλ(AT)i obtainedbyprependingΩ(cid:202)totheλ-generatingfunction. Wepresentanexamplethatwealsousefor thegeometricinterpretation.LetA=[ 2 3 ]andb=5andconsidertheDiophantinesystemAx(cid:202)b. Then λ−5 Φλ (z ,z ,λ)= (cid:88) λ2x1+3x2−5zx1zx2 and ρλ (z ,z ,λ)= . A,b 1 2 x1,x2∈(cid:90)(cid:202)0 1 2 A,b 1 2 (1−z1λ2)(1−z2λ3) Thus,thecrudegeneratingfunctionforthissystemis λ−5 Ω (cid:202) (1−z λ2)(1−z λ3) 1 2 Afterdefiningthecrudegeneratingfunction, MacMahonproceedstoevaluatethisexpressionbythe useofformalrules,suchas 1 1 Ω(cid:202)(1−λx)(cid:161)1− y (cid:162)=(1−x)(1−xsy). λs Inwhatfollowswewillexplainwhatthismeansgeometrically. Conesandfundamentalparallelepipeds. Theconnectionbetweenpartitionanalysisandgeometryhinges onthefactthatcertainrationalfunctionexpressionscorresponddirectlytogeneratingfunctionsoflat- ticepointsinpolyhedralcones. Alatticepoint isanyintegerlinearcombinationofafixedsetofbasis vectors.Withoutlossofgenerality,inthisarticle,wewillalwaysworkwiththethestandardunitvectors asourbasis,sothat“latticepoint”issynonymouswith“integerpoint”. WealreadyusedΦtodenotetheformalpowerseriesgeneratingfunctionforlinearDiophantinesys- tems.Moregenerally,forasetS⊆(cid:90)d oflatticepointswewilluseΦ todenotethegeneratingfunction S ofS,i.e.,ΦS=(cid:80)s∈Szs.Forconvenience,ifS⊆(cid:82)d isasetofrealvectors,wewillalsowriteΦS todenote thegeneratingfunctionofalllatticepointsinS,namelyΦS =(cid:80)s∈S∩(cid:90)dzs. Ofcourse,ifandwhereΦS 8 FELIXBREUERANDZAFEIRAKISZAFEIRAKOPOULOS (A) C(cid:90)((1,0),(1,3)). p(Bo)sed(cid:90)2i∩ntCo(cid:82)(t(h1r,e0e),(t1r,a3n))sladteecsomo-f C(C(cid:82))((1T,h0e),(c1o,3rr)e)sbpyotnrdaninsglatteilsinogfiotsf C(cid:90)((1,0),(1,3)). fundamentalparallelepiped. FIGURE2. Generatingfunctionsofcones. convergesdependsonthesetS,butthesequestionswillbeofnoparticularconcerntous,aswewillsee below.WeareprimarilyinterestedingeneratingfunctionsΦ whereC isasimplicialpolyhedralcone. C WedefinetheconeC generatedbyv ,v ,...,v as 1 2 d (cid:40) (cid:175) (cid:41) C(cid:82)(v1,...,vd)= (cid:88)d λivi (cid:175)(cid:175)(cid:175)0(cid:201)λi ∈(cid:82) . i=1 (cid:175) Thev arethegeneratorsofC. Ifthev arelinearlyindependent,theconeC issimplicial. Allconesin i i thisarticlearerationalpolyhedra,whichmeansthattheycanalsobedefinedasthesetofrealvectors satisfyingafinitesystemoflinearinequalitieswithintegercoefficients.AconeC ispointedorline-free ifitdoesnotcontainaline. Allsimplicialconesarepointed. ForpointedconesC thepowerseriesΦ C alwayshasanon-trivialradiusofconvergence,seealsoSection2.4. IfagivenconeC ispointedand thechosensetofgeneratorsvi isinclusion-minimal,thenthesets(cid:82)(cid:202)0vi arecalledtheextremeraysof thecone. ThesetofextremeraysR ofapointedconeisuniquelydetermined,butwithineachraythe i choiceofgenerators0(cid:54)=v ∈R isarbitrary.Avectorv∈(cid:90)n isprimitive,ifthegreatestcommondivisor i i ofitsentriesis1. Foreveryvector0(cid:54)=v ∈(cid:81)n thereexistsauniqueprimitivevectorprim(v)thatisa positivemultipleofv. Equivalently,prim(v)istheuniqueshortestnon-zerointegervectorontheray (cid:82)(cid:202)0v. Withthesenormalizationswefindthateachpointedrationalconehasauniqueminimalsetof primitivegenerators. Later, wearealsogoingtoneedthenotionofaffinecones, whichjustmeansconeslikethosede- finedabovetranslatedbyarationalvectorq.Fortranslatesofpointedcones,thetranslationvectorqis uniquelydetermined,inwhichcaseq iscalledtheapexoftheaffinecone. Wewillfrequentlydropthe adjective“affine”andjustrefertoconesingeneral. Itwillbeclearfromcontextwhetheragivencone hasitsapexattheoriginornot.Foragivenq∈(cid:81)d weadoptthenotation C(cid:82)(cid:161)v1,...,vd;q(cid:162)=q+C(cid:82)(v1,...,vd). Here,+denotestheMinkowskisum,i.e.,fortwosetsAandBwehaveA+B={a+b|a∈A,b∈B}. Ifinsteadofconsideringallconiccombinationsofthegeneratorswithrealcoefficients,weconsider onlynon-negativeintegercombinations,thenweobtainthesemigroup (cid:40) (cid:175) (cid:41) C(cid:90)(v1,...,vd)= (cid:88)d λivi (cid:175)(cid:175)(cid:175)0(cid:201)λi ∈(cid:90) . i=1 (cid:175) which we also call the discrete cone generated by the v . As above we write C (cid:161)v ,...,v ;q(cid:162):= q+ i (cid:90) 1 d C (v ,...,v ). (cid:90) 1 d Asanexample,considerthesemigroupgeneratedby(1,0)and(1,3)asshowninFigure2a.Ifwetake thegeometricseries 1 andexpandit,i.e.,1+z +z2+...,weobservethatthisseriesisthegenerating (1−z1) 1 1 functionofthelatticepointsonthenon-negativex-axis. Inordertotakethegeneratingfunctionfor thelatticepointsatheight3,weneedtomultiplytheseriesby(cid:161)z z3(cid:162). Accordingly,totakethelattice 1 2 POLYHEDRALOMEGA:SOLVINGLINEARDIOPHANTINESYSTEMS 9 pointsatheight6wemultiplyby(cid:161)z z3(cid:162)2andsoon.Itiseasytosee,thatthegeneratingfunctionofthe 1 2 semigroupistheproductoftwogeometricseries,namely (1−1z1) and (cid:161)1−z11z23(cid:162),sothatweobtain 1 Φ = . C(cid:90)((1,0),(1,3)) (1−z )(cid:161)1−z z3(cid:162) 1 1 2 Followingthesameidea,wecantranslatethewholesemigroupbymultiplyingbyamonomial.Ifwe multiply(1−z1)(cid:161)11−z1z23(cid:162)byz1z2weobtainthegeneratingfunctionoftheredpointsinFigure2b.Similarly, we multiply by z z2 to obtain the generating function of the green points. Now, since there are no 1 2 overlaps,wecansumthegeneratingfunctionsinordertoobtainthegeneratingfunctionforalllattice pointsintheconegeneratedby(1,0)and(1,3)withrealcoefficients, 1+z z +z z2 Φ =(cid:161)1+z z +z z2(cid:162)·Φ = 1 2 1 2 . C(cid:82)((1,0),(1,3)) 1 2 1 2 C(cid:90)((1,0),(1,3)) (1−z )(cid:161)1−z z3(cid:162) 1 1 2 Thenumeratorofthisrationalfunctionencodesthelatticepoints(0,0),(1,1)and(1,2)whicharepre- ciselythelatticepointsintheparallelogramΠ={λ (1,0)+λ(1,3)|0(cid:201)λ <1}asshowninblueinFig- 1 i ure2c.Moreover,theconeC ((1,0),(1,3))istiledbytranslatesofΠ.Thisobservationgeneralizes. (cid:82) GivenasimplicialconeC=C (v ,v ,...,v )∈(cid:82)d,wedefineitsfundamentalparallelepipedas (cid:82) 1 2 d (cid:40) (cid:175) (cid:41) Π(cid:82)(C)=Π(cid:82)(v1,...,vd)= (cid:88)d λivi (cid:175)(cid:175)(cid:175)0(cid:201)λi <1,λi ∈(cid:82) . i=1 (cid:175) andthesetoflatticepointsinthefundamentalparallelepipedasΠ (C)=Π (C)∩(cid:90)d.Giventhisdefini- (cid:90)d (cid:82) tion,thegeneralformofthegeneratingfunctionofasimplicialconeC=C (v ,...,v )⊂(cid:82)d generated (cid:82) 1 m byv ,...,v ⊂(cid:90)d is 1 m Φ (z ,...,z ) Π (C) 1 d (2) Φ (z ,...,z ) = (cid:90)d C 1 d (cid:81)mi=1(1−zvi) where Φ (z ,...,z ) is the generating function for the lattice points in the fundamental paral- Π (C) 1 d (cid:90)d lelepipedofC. ThisfundamentalobservationgoesbacktoEhrhart[36,37,38]. Fromthisidentityit isimmediatelyobviousthatconesC withfewlatticepointsintheirfundamentalparallelepipedΠ (C) (cid:90)d have particularly short representations as rational functions. ConesC where Π (C) contains just a (cid:90)d singlelatticepointarecalledunimodular. Animportantpropertyofthefundamentalparallelepipedisthatitslatticepointsarecosetrepresen- tativesforthelatticepointsintheconemodulothesemigroupofthegenerators(asexemplifiedbythe constructionofthegeneratingfunction). Toputitinanotherway: thesetoflatticepointsinthecone is the set of lattice points in the fundamental parallelepiped, translated by all non-negative integral combinationsofthegeneratorsofthecone,i.e., (cid:195) (cid:33) n (3) C∩(cid:90)d = (cid:71) (cid:88)ijvj+Π(cid:90)d(C) = C(cid:90)(v1,...,vn)+Π(cid:90)d(C) i1,i2,...,in∈(cid:90)(cid:202)0 j=1 (cid:70) Here denotesanordinaryuniontogetherwiththeassertionthattheoperandsaredisjoint. GeometryoftheΩ(cid:202)operatorandtheMacMahonlifting. NowwecanseetheconnectionofMacMahon’s constructiontopolyhedralgeometry. InFigure3,weseethegeometricstepsforcreatingtheλ-cone (equivalentoftheλ-generatingfunction)bymeansoftheMacMahonliftingandthenapplyingtheΩ(cid:202) operator. Given an LDS Ax (cid:202)b with m inequalities and d variables, we first lift the standard generators of thepositiveorthantof(cid:82)d byappendingtheexponentsoftheλvariables(thusliftingthegeneratorsin (cid:82)d+m). ThisliftingideafirstappearsinananalyticcontextinCayley’swork[27]. Considerthecone (cid:179) (cid:180) generatedbytheMacMahonlifting,i.e.,theconegeneratedbythecolumnsofmatrixV = Idd . The A generatingfunctionofthiscone,dueto(2),isexactly 1 . ThisistruebecausethematrixV is (cid:179)1−ziλ(AT)i(cid:180) unimodularbyconstruction,whencethefundamentalparallelepipedoftheconecontainsexactlyone 10 FELIXBREUERANDZAFEIRAKISZAFEIRAKOPOULOS λ λ y λ λ x y y x y x x (A) Consider the first (D) The projected quadrant and lift the (B) Shift according to (C) Take the part of polyhedron contains standard generators the inhomogeneous the cone (a polyhe- the solutions to the accordingtotheinput. part in he negative dron)livingabovethe inequality. λ direction the cone xy-plane and project generatedbythelifted itslatticepoints. generators. Intersect withthexy-plane. FIGURE3. Ω(cid:202)asliftingandprojection. latticepoint,namelytheorigin. Next,wetranslatetheconebyq =(cid:161) 0 (cid:162)anddenotethenewconeby −b C=C(V)+(cid:161) 0 (cid:162).ThenthegeneratingfunctionofconeC isexactlytheλ-generatingfunction,i.e., −b Φ (z)=ρλ (z)=λ−b(cid:89)m 1 . C A,b i=1(cid:179)1−ziλ(AT)i(cid:180) ApplyingtheOmegaoperatortoΦ nowcorrespondsgeometricallytointersectingC withallhalf- C spaceswheretheλ-variablesarenon-negativeandthenapplyingaprojectionthatforgetstheλ-variables. Tomakethisprecise,letusdenotebyHλtheintersectionofcoordinatehalfspaceswhereallλvariables arenon-negative,i.e., Hλ=(cid:110)(cid:161)x1,x2,...,xd,xλ1,xλ2,...,xλm(cid:162)∈(cid:82)d+m(cid:175)(cid:175)(cid:175)xλi (cid:202)0for1(cid:201)i(cid:201)m(cid:111). Next,letπdenotetheprojectionfrom(cid:82)d+mto(cid:82)d thatforgetsthelastmcoordinateswhichcorrespond toλ-variables. Giventhisnotation,theactionofΩ(cid:202)ontheλ-generatingfunctioncorrespondstofirst intersectingC withtheHλandthenprojectingtheresultingpolyhedronusingπ,i.e., Ω(cid:202)ΦC(z)=Ω(cid:202)ρλA,b(z)=Φπ(C∩Hλ)(z). Continuingwiththeexamplewehadchosenpreviously,assumeA=[23]andb=5.InFigure3a,we constructthevectors(1,0,2)and(0,1,3),liftingthetwostandardbasisvectors(1,0)and(0,1)according totheentriesofmatrix A. InFigure3b,weshifttheconegeneratedby(1,0,2)and(0,1,3)by−b=−5, and consider the halfspace defined by the xy-plane, where λ is non-negative. The two intersection pointsoftheconewiththehyperplanearemarked. Now,inFigure3c,weconsidertheintersectionof thehalfspacedefinedabovewiththeshiftedcone. Thelatticepointsintheintersectionareshownin bluemarked,andtheirprojectionsontotheλ=0planeareshowninred. Theseprojectionsarethe non-negativesolutionstotheoriginalinequality,asshowninFigure3d. MacMahon’s rules. MacMahon’s method was general and in principle algorithmic, based on Elliott’s work (see next section). As we shall see, Elliott’s method is simple but computationally very ineffi- cient. MacMahon introduced a set of rules in order to deal with particular combinatorial problems computationally–inparticular,ashehadtocarryouthiscomputationsbyhand. Wepresentsomeof MacMahon’srulesfromageometricpointofview.