Polyhedral Omega: A New Algorithm for Solving Linear Diophantine Systems

0.67 MB·English

by  Felix Breuer

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 poly-
hedral 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.