Internals

Debugging

ToricGB`Debug`HemmeckeToricGroebnerBasis=1;
ToricGB`Debug`geometricBuchberger=2;
ToricGB`Debug`binomialReduce=4;

Flags to enable assertions, profiling and tune debug output.

ToricGB`Debug`assert=False;
ToricGB`Debug`logging=False;
ToricGB`Debug`profile=False;
ToricGB`Debug`output=BitOr[0(*,ToricGB`Debug`HemmeckeToricGroebnerBasis,ToricGB`Debug`geometricBuchberger,ToricGB`Debug`binomialReduce*)];

Main Debugging Functions

assert is  always enabled, dbgAssert depends on ToricGB`Debug`assert. dbgPrint receives as parameter from which function it is called and prints depending on the bits set in ToricGB`Debug`output. dbgLog writes messages to a file depending on ToricGB`Debug`logging.

SyntaxInformation[assert]={"ArgumentsPattern"->{_,_, __}};
SetAttributes[assert,HoldAll];
assert[condition_,function_,error__]:=If[!condition,Print["Error in ",Unevaluated[function],": ",error];Abort[]]

SyntaxInformation[dbgAssert]={"ArgumentsPattern"->{_,_, __}};
SetAttributes[dbgAssert,HoldAll];
dbgAssert[condition_,function_,error__]:=If[ToricGB`Debug`assert,If[!condition,Print["Error in ",Unevaluated[function],": ",error];Abort[]]]

SyntaxInformation[dbgPrint]={"ArgumentsPattern"->{_,__}};
SetAttributes[dbgPrint,HoldAll];
dbgPrint[switch_,message__]:=If[BitAnd[ToricGB`Debug`output,switch]≠0,Print[message]]

If[ValueQ[Internal`toricGBlogStream],Close[Internal`toricGBlogStream]];
Internal`toricGBlogStream=OpenWrite["log.txt"];
SyntaxInformation[dbgLog]={"ArgumentsPattern"->{__}};
SetAttributes[dbgLog,HoldAll];
dbgLog[message__]:=If[ToricGB`Debug`logging,WriteString[Internal`toricGBlogStream,message];WriteString[Internal`toricGBlogStream,"\n"]]

Verification Functions

The macros LM and S compute the leading monomial and the S-polynomial of a polynomial.

SyntaxInformation[LM]={"ArgumentsPattern"->{_,_, _}};
LM[f_,order_,vars_]:=First[MonomialList[f,vars,order]]

SyntaxInformation[S]={"ArgumentsPattern"->{_,_, _,_}};
S[f_,g_,order_,vars_]:=Expand[(f*LM[g,order,vars]-g*LM[f,order,vars])/PolynomialGCD[LM[f,order,vars],LM[g,order,vars]]]

dbgGroebnerBasisUnique brings a reduced Gröbner basis into a unique representation by sorting the polynomials and assigning unique signs. The leading coefficients are assumed to be 1.

SyntaxInformation[dbgGroebnerBasisUnique]={"ArgumentsPattern"->{_,_}};
dbgGroebnerBasisUnique[G_List,variables_List]:=Sort[Map[(If[(LM[#,Lexicographic,variables]/.Thread[variables->1])>0,#,-#])&,G]]

dbgGroebnerBasisQ checks whether the G (given in vector form) is a Gröbner basis with respect to monomialOrder and projection.

SyntaxInformation[dbgGroebnerBasisQ]={"ArgumentsPattern"->{_,_, _}};dbgGroebnerBasisQ[G_List,monomialOrder_List,projection_List]:=Module[{GP,mop,n=Length[projection[[ToricGB`Private`vectorList]]],GPP,i,j},
GP=Map[(#[[projection[[ToricGB`Private`vectorList]]]])&,binomialStrip[G]];
mop=Map[(#[[projection[[ToricGB`Private`vectorList]]]])&,monomialOrder];
mop=mop[[latticeZeroKernelProjection[Transpose[mop]][[ToricGB`Private`vectorList]]]];
GPP=VectorToBinomial[GP,vars[n]];
For[i=1,i≤Length[GPP],i++,
For[j=i+1,j≤Length[GPP],j++,
If[PolynomialReduce[S[GPP[[i]],GPP[[j]],mop,vars[n]],GPP,vars[n],MonomialOrder->mop][[2]]=!=0,
Print["groebnerBasisQ: G=",G,", monomialOrder=",monomialOrder,", projection=",projection," fails at ",{i,j}];
Return[False]];
];
];
True
]

Profiling

resetProfiling initializes profiling variables, printProfiling outputs them. addTime measures the execution time of a command, adds it to a specified variable and returns the value of the function. profile initializes the profiling, then executes the specified function, measures its time and prints the profiling variables.

resetProfiling:=(timeGroebner=0;timeRedGB=0;timePrereduce=0;timeCM=0;timeInterRed=0;timeMinimize=0;timeRebuild=0;
cntCS=0;cntCT=0;cntCM=0;
cntReduce0=0;cntReduce1=0;cntReduce2=0;cntReduce2Q=0;
cntReduceQDivQ=0;cntReduceDivQ=0;cntReduceDiv=0;cntMSDivQ=0;cntDividesIter=0;
cntInsert=0;cntChange=0;cntDelete=0;cntQuery=0;cntQueryFail=0;cntRebuild=0;cntKDCallback=0;cntKDRecursions=0;cntCompare=0;cntBinaryIter=0;sizeG=0;
timeDiv=0;cntDiv=0;maxExponent={})
resetProfiling;

printProfiling:=Print["timePrereduce=",timePrereduce,", timeGroebner=",timeGroebner,", timeRedGB=",timeRedGB,", timeCM=",timeCM,", timeInterRed=",timeInterRed,", timeMinimize=",timeMinimize,", timeRebuild=", timeRebuild,", cntCS=",cntCS,", cntCT=",cntCT,", cntCM=",cntCM,", cntReduce0=",cntReduce0,", cntReduce1=",cntReduce1,", cntReduce2=",cntReduce2,", cntReduce2Q=",cntReduce2Q,", cntReduceQDivQ=",cntReduceQDivQ,", cntReduceDivQ=",cntReduceDivQ,", cntReduceDiv=",cntReduceDiv,", cntMSDivQ=",cntMSDivQ,", cntDividesIter=",cntDividesIter,", cntInsert=",cntInsert,", cntChange=",cntChange,", cntDelete=",cntDelete,", cntQuery=",cntQuery,", cntQueryFail=",cntQueryFail,", cntRebuild=",cntRebuild,", cntKDCallback=",cntKDCallback,", cntKDRecursions=",cntKDRecursions,", cntCompare=",cntCompare,", cntBinaryIter=",cntBinaryIter];

SetAttributes[addTime,HoldAll];
addTime[func_,var_]:=Module[{return,time},
{time,return}=AbsoluteTiming[func];
var+=time;
return
];

SetAttributes[profile,HoldAll];
profile[func_]:=(resetProfiling;Print[AbsoluteTiming[func][[1]]];printProfiling;)

Common Parameters for Examples

Functions that give the lexicographic order for n variables, the n-th character in the alphabet as symbol, the first n characters and the projection from n variables to the first k variables.

lexOrder[n_]:=Array[(UnitVector[n,#])&,n];
var[n_]:=Symbol[FromCharacterCode[n-1+ToCharacterCode["a"]]];
vars[n_]:=Array[(var[#])&,n];
proj[k_,n_]:=positionsToBitMask[Table[i,{i,k}],n];

Data Structures

Queue (from http://library.wolfram.com/conferences/devconf99/lichtblau)

An implementation of an queue using DownValues.

ToricGB`Private`queuePointer=1;
ToricGB`Private`queueElement=2;
ToricGB`Private`queueFront=3;
ToricGB`Private`queueBack=4;
ClearAll[newQueue,clearQueue,queueLength,queueEmpty,enQueue,deQueue];

SyntaxInformation[newQueue]={"ArgumentsPattern"->{}};
newQueue[]:={Unique[],Null,0,0};

SyntaxInformation[clearQueue]={"ArgumentsPattern"->{_}};
SetAttributes[clearQueue,HoldFirst];
clearQueue[Q_Symbol]:=(ClearAll[Evaluate[Q[[ToricGB`Private`queuePointer]]]];ClearAll[Q])

SyntaxInformation[queueLength]={"ArgumentsPattern"->{_}};
SetAttributes[queueLength,HoldFirst];
queueLength[Q_Symbol]:=(Q[[ToricGB`Private`queueBack]]-Q[[ToricGB`Private`queueFront]])

SyntaxInformation[queueEmpty]={"ArgumentsPattern"->{_}};
SetAttributes[queueEmpty,HoldFirst];
queueEmpty[Q_Symbol]:=(Q[[ToricGB`Private`queueBack]]===Q[[ToricGB`Private`queueFront]])

SyntaxInformation[enQueue]={"ArgumentsPattern"->{_,_}};
SetAttributes[enQueue,HoldFirst];
enQueue[Q_Symbol,elem_]:=(Q[[ToricGB`Private`queueFront]]++;Q[[ToricGB`Private`queuePointer]][Q[[ToricGB`Private`queueFront]]]=elem;)

SyntaxInformation[deQueue]={"ArgumentsPattern"->{_}};
SetAttributes[deQueue,HoldFirst];
deQueue[Q_Symbol]:=(Q[[ToricGB`Private`queueBack]]++;Q[[ToricGB`Private`queueElement]]=Q[[ToricGB`Private`queuePointer]][Q[[ToricGB`Private`queueBack]]];Q[[ToricGB`Private`queuePointer]][Q[[ToricGB`Private`queueBack]]]=.;Q[[ToricGB`Private`queueElement]])

Binary Search

Search the (first) position in which x should be inserted into L such that it remains ordered. If x is contained in L, the position of x in L is returned.

SyntaxInformation[binarySearch]={"ArgumentsPattern"->{_,_}};
binarySearch[L_List,x_,start_Integer:1]:=Module[{a=start,b=Length[L]+1,m},
While[a<b,
(*cntBinaryIter++;*)
m=Quotient[a+b,2];
If[x>L[[m]],a=m+1,b=m];
];
Return[a];
];

binarySearch[{1,4,6,9},5]

binarySearch[{1,4,6,9},9]

binarySearch[{1,4,6,9},12]

kd-List

A kd-list is a k-dimensional (linked) list, similar to kd-trees, which is used to store k-dimensional integer vectors. The first level is a list containing the first component of each vector (each value is only stored once). Then each entry x points to a list of level two in which we store the second component of the vectors whose first component has value x. This is continued up to level k.

This implementation maintains the lists at each level in a sorted order, resulting in slow insertion times but very fast queries of the kind "all entries with components smaller x". Behaviour of this kind is desirable since the list is rarely updated and many queries are run. It operates on binomials with extra information, v[[ToricGB`Private`binomialExponent]] being the integer vector used as "key" and v[[ToricGB`Private`binomialIndex]] being the index of the vector in the plain list that is maintained for direct access. Additionally, only the components listed in projection are used for branching. Those vectors whose coordinate is negative or zero are stored in one bucket with value 0 (this is possible since all queries are nonnegative).

It is possiple to insert, change and (kind of) delete entries. The functions with the parameter hint are helper functions that only perform a partial job and not meant to be called externally.

The function queryKDList iterates over all entries of the kd-list that are componentwise less or equal than the query. Depending on the parameter largeFirst and whether the query value is positive or zero, it starts with the vectors with large respectively small entries. For each entry the callback function is called. As soon as the callback function returns a value different from Null, the query is aborted and the value of the callback function is returned. The vector query must be nonnegative. If inverseQuery is True, instead all elements that are componentwise greater or equal than the query are looked for.

There are functions to remove and change elements in the kd-list. rebuildKDList is only provided to conform with the implementations of kd-trees and Mathematica lists and does nothing.

The first entry of a KDList is a plain list of all elements. The second stores the projection which is passed to newKDList as parameter. The head symbol of the KDList is stored in the third entry of the list returned by newKDList. Each element has an entry ToricGB`Private`binomialIndex which corresponds to its position in the plain list. An element with index i is considered as disabled (temporarily removed) if KDL[[ToricGB`Private`basisKDList,i,ToricGB`Private`binomialMarker]]=!=1.

Attention

Inserting an element twice breaks the data structure!

ToricGB`Private`basisList=1;
ToricGB`Private`basisProjection=2;
ToricGB`Private`basisKDList=3;
ToricGB`Private`basisChanges=4;
ClearAll[newKDList,clearKDList,insertKDList,deleteKDList,changeKDList,queryKDList];

SyntaxInformation[newKDList]={"ArgumentsPattern"->{_}};
newKDList[projection_List]:=Module[{KDL},
(*dbgLog["New kd list with projection=",projection];*)
KDL={{},projection,Unique[]};
KDL[[ToricGB`Private`basisKDList]][{}]={};
KDL
]

SyntaxInformation[clearKDList]={"ArgumentsPattern"->{_}};
SetAttributes[clearKDList,HoldFirst];
clearKDList[KDL_Symbol]:=(ClearAll[Evaluate[KDL[[ToricGB`Private`basisKDList]]]];ClearAll[KD])

SyntaxInformation[rebuildKDList]={"ArgumentsPattern"->{_}};
SetAttributes[rebuildKDList,HoldFirst];
rebuildKDList[KDT_Symbol,idx_Integer]:=idx;

SyntaxInformation[insertKDList]={"ArgumentsPattern"->{_,_,_.}};
SetAttributes[insertKDList,HoldFirst];
insertKDList[KDL_Symbol,v_List,hint_List]:=Module[{it,pos,level=Length[hint]+1,coord,value},
If[level>Length[KDL[[ToricGB`Private`basisProjection,ToricGB`Private`vectorList]]],
(* end of projection -> simple list *)
KDL[[ToricGB`Private`basisKDList]][hint]=Append[KDL[[ToricGB`Private`basisKDList]][hint],v];
,
(* search v[[ToricGB`Private`binomialExponent,coord]] in level 'hint' *)
coord=KDL[[ToricGB`Private`basisProjection,ToricGB`Private`vectorList,level]];
value=Max[v[[ToricGB`Private`binomialExponent,coord]],0];
it=Append[hint,value];
pos=binarySearch[KDL[[ToricGB`Private`basisKDList]][hint],value];
If[pos>Length[KDL[[ToricGB`Private`basisKDList]][hint]]||KDL[[ToricGB`Private`basisKDList]][hint][[pos]]=!=value,
(* not found -> insert *)
KDL[[ToricGB`Private`basisKDList]][hint]=Insert[KDL[[ToricGB`Private`basisKDList]][hint],value,pos];
KDL[[ToricGB`Private`basisKDList]][it]={};
];
(* go to next level *)
insertKDList[KDL,v,it];
];
]

insertKDList[KDL_Symbol,v_List]:=Module[{it,pos,coord,u=v,value},
(*cntInsert++;
dbgLog["insertKDList[",v,"]"];*)
u[[ToricGB`Private`binomialIndex]]=Length[KDL[[ToricGB`Private`basisList]]]+1;
u[[ToricGB`Private`binomialMarker]]=1;
(* append to list of all elements *)
KDL[[ToricGB`Private`basisList]]=Append[KDL[[ToricGB`Private`basisList]],u];
(* insert in real kd-list *)
insertKDList[KDL,u,{}];
(*dbgLog["G=",KDL[[ToricGB`Private`basisList]]];*)
]

SyntaxInformation[deleteKDList]={"ArgumentsPattern"->{_,_,_}};
SetAttributes[deleteKDList,HoldFirst];
(* returns True if list KDL[[ToricGB`Private`basisKDList]][hint] === {} after call *)
deleteKDList[KDL_Symbol,v_List,hint_]:=Module[{it,pos,level=Length[hint]+1,coord,value},
If[level>Length[KDL[[ToricGB`Private`basisProjection,ToricGB`Private`vectorList]]],
(* end of projection -> simple list *)
KDL[[ToricGB`Private`basisKDList]][hint]=Select[KDL[[ToricGB`Private`basisKDList]][hint],(#[[ToricGB`Private`binomialIndex]]=!=v[[ToricGB`Private`binomialIndex]])&];
,
(* search v[[ToricGB`Private`binomialExponent,coord]] in level 'hint' *)
coord=KDL[[ToricGB`Private`basisProjection,ToricGB`Private`vectorList,level]];
value=Max[v[[ToricGB`Private`binomialExponent,coord]],0];
it=Append[hint,value];
pos=binarySearch[KDL[[ToricGB`Private`basisKDList]][hint],value];
(*assert[pos<=Length[KDL[[ToricGB`Private`basisKDList]][hint]]&&KDL[[ToricGB`Private`basisKDList]][hint][[pos]]===value,deleteKDList,"Could not find ",v," in kd-list."];*)
(* go to next level *)
If[deleteKDList[KDL,v,it],
KDL[[ToricGB`Private`basisKDList]][it]=.;
KDL[[ToricGB`Private`basisKDList]][hint]=Complement[KDL[[ToricGB`Private`basisKDList]][hint],{value}];
];
];
KDL[[ToricGB`Private`basisKDList]][hint]==={}
];

deleteKDList[KDL_Symbol,index_Integer]:=Module[{u=v},
(*cntDelete++;
dbgLog["deleteKDList[",index,"]"];*)
KDL[[ToricGB`Private`basisList,index,ToricGB`Private`binomialMarker]]=0;
deleteKDList[KDL,KDL[[ToricGB`Private`basisList,index]],{}];
(*dbgLog["G=",KDL[[ToricGB`Private`basisList]]];*)
];

SyntaxInformation[changeKDList]={"ArgumentsPattern"->{_,_,_}};
SetAttributes[changeKDList,HoldFirst];
changeKDList[KDL_Symbol,index_Integer,v_List]:=Module[{it,pos,coord,u=v,value},
(*cntChange++;
dbgLog["changeKDList[",index,", ",v,"]"];*)
deleteKDList[KDL,KDL[[ToricGB`Private`basisList,index]],{}];
u[[ToricGB`Private`binomialIndex]]=index;
u[[ToricGB`Private`binomialMarker]]=1;
KDL[[ToricGB`Private`basisList,index]]=u;
insertKDList[KDL,u,{}];
(*dbgLog["G=",KDL[[ToricGB`Private`basisList]]];*)
]

SyntaxInformation[queryKDList]={"ArgumentsPattern"->{_,_,_,_,_.,_.,_.}};
SetAttributes[queryKDList,HoldAll];
queryKDList[KDL_Symbol,query_,cnt_,callback_,inverseQuery_,largeFirst_,hint_]:=Module[{it,level=Length[hint]+1,coord,pos,return,value},
(*cntKDRecursions++;*)
If[level>Length[KDL[[ToricGB`Private`basisProjection,ToricGB`Private`vectorList]]],
(* end of projection -> simple list *)
For[pos=1,pos≤Length[KDL[[ToricGB`Private`basisKDList]][hint]],pos++,
(*cntKDCallback++;*)
return=callback[KDL,query,KDL[[ToricGB`Private`basisKDList]][hint][[pos]]];
If[return=!=Null,Sow[return];If[--cnt==0,Return[]]];
],
(* traverse list in level hint up to value v[[coord]] *)
coord=KDL[[ToricGB`Private`basisProjection,ToricGB`Private`vectorList,level]];
value=query[[ToricGB`Private`binomialExponent,coord]];
it=Append[hint,0];
If[inverseQuery,
(* look for greater values than query *)
For[pos=Length[KDL[[ToricGB`Private`basisKDList]][hint]],pos≥1&&KDL[[ToricGB`Private`basisKDList]][hint][[pos]]>=value,pos--,
it[[level]]=KDL[[ToricGB`Private`basisKDList]][hint][[pos]];
queryKDList[KDL,query,cnt,callback,inverseQuery,largeFirst,it];
If[cnt==0,Return[]];
];
,
(* look for smaller values than query *)
If[!largeFirst||value≤0,
(* start search with small values *)
For[pos=1,pos<=Length[KDL[[ToricGB`Private`basisKDList]][hint]]&&KDL[[ToricGB`Private`basisKDList]][hint][[pos]]≤value,pos++,
it[[level]]=KDL[[ToricGB`Private`basisKDList]][hint][[pos]];
queryKDList[KDL,query,cnt,callback,inverseQuery,largeFirst,it];
If[cnt==0,Return[]];
];
,
(* start search with large values -> binary search *)
pos=binarySearch[KDL[[ToricGB`Private`basisKDList]][hint],value];
If[pos>Length[KDL[[ToricGB`Private`basisKDList]][hint]]||KDL[[ToricGB`Private`basisKDList]][hint][[pos]]>value,pos--];
While[pos>0,
it[[level]]=KDL[[ToricGB`Private`basisKDList]][hint][[pos]];
queryKDList[KDL,query,cnt,callback,inverseQuery,largeFirst,it];
If[cnt==0,Return[]];
pos--;
];
];
];
];
]

queryKDList[KD_Symbol,query_,count_,callback_,inverseQuery_:False,largeFirst_:False]:=Module[{cnt=count,results},
(*cntQuery++;*)
results=Reap[queryKDList[KD,query,cnt,callback,inverseQuery,largeFirst,{}]][[2]];(*dbgLog["queryKDList[",query,", ",callback,", ",count,", ",inverseQuery,", ",largeFirst,"] = ",results];
If[results==={},cntQueryFail++];*)
results
]

kdl=newKDList[proj[3,3]]
insertKDList[kdl,binomialSupport[{4,-4,10}]];
insertKDList[kdl,binomialSupport[{4,-1,2}]];
insertKDList[kdl,binomialSupport[{4,3,1}]];
insertKDList[kdl,binomialSupport[{4,5,-1}]];

Information[Evaluate[kdl[[ToricGB`Private`basisKDList]]]]

queryKDList[kdl,binomialSupport[{4,4,9}],Infinity,(#3)&]

changeKDList[kdl,2,binomialSupport[{5,4,2}]];

Information[Evaluate[kdl[[ToricGB`Private`basisKDList]]]]

queryKDList[kdl,binomialSupport[{4,0,4}],1,(#3)&,False]

queryKDList[kdl,binomialSupport[{4,0,4}],1,(#3)&,True]

kd-Tree

This is an implementation of the classical kd-trees. The elements are only stored in the leaves of the tree. At each level one coordinate is used to split the tree into two subtrees with nonoverlapping value intervals. Usually the coordinates are used in a circular fashion. If splitting with respect to a coordinate is not possible, it is skipped. If less than ToricGB`Private`kdtCutOff elements are remaining in the subtree, they are stored as a plain list.

Each inner node consist of a list containing the splitting index, the value v used for splitting (if this node represented the interval [a,b) in this coordinate, the left subtree represents the interval [a,v), the right [v,b) ). The following four values represent the minimum resp. maximum actual value in both subtrees. They might be changed by insertKDTree and deleteKDTree (as opposed to v) (the latter is not implemented). Those are also the values used by queryKDTree.

The function rebuildKDTree removes all disabled elements and recreates the tree from scratch. It returns the index of the element, whose index was passed to the function - after the changes made.

queryTree maintains a list of coordinates which still have to be checked. In theory, computing the Complement each time an entry is removed could be avoided by using a doubly linked list (or vector) of the entries and an array with pointers to the entries of the list which can be accessed by the coordinate value. I decided against the implementation of this technique due to the overhead introduced by Mathematica.

The first three entry of a KDTree are like in KDList. The fourth counts the changes (insertions, changes, deletions) to the data structure. This value is used to decide in which case the data structure is rebuild.

The rest of the implementation is analogous to KDList.

ClearAll[newKDTree,clearKDTree,rebuildKDTree,insertKDTree,deleteKDTree,changeKDTree,queryKDTree];
ToricGB`Private`kdtRangeStart=1;
ToricGB`Private`kdtRangeEnd=2;
ToricGB`Private`kdtIndex=1;
ToricGB`Private`kdtSeparator=2;
ToricGB`Private`kdtLowMin=3;
ToricGB`Private`kdtLowMax=4;
ToricGB`Private`kdtHighMin=5;
ToricGB`Private`kdtHighMax=6;
ToricGB`Private`kdtCutOff=1;

SyntaxInformation[newKDTree]={"ArgumentsPattern"->{_}};
newKDTree[projection_List]:=Module[{KDT},
(*dbgLog["New kd tree with projection=",projection];*)
KDT={{},projection,Unique[],0};
KDT[[ToricGB`Private`basisKDList]][Table[{0,Infinity},{Length[projection[[ToricGB`Private`vectorList]]]}]]={};
KDT
]

SyntaxInformation[clearKDTree]={"ArgumentsPattern"->{_}};
SetAttributes[clearKDTree,HoldFirst];
clearKDTree[KDT_Symbol]:=(ClearAll[Evaluate[KDT[[ToricGB`Private`basisKDList]]]];ClearAll[KDT])

SyntaxInformation[rebuildKDTree]={"ArgumentsPattern"->{_,_.,_.,_.}};
SetAttributes[rebuildKDTree,HoldFirst];
rebuildKDTree[KDT_Symbol,F_List,hint_List,levelP_Integer]:=Module[{coord,G,mid,value,it,level=levelP},
(*dbgPrint["rebuildKDTree[",F,",",hint,",",level,"]"];*)
If[Length[F]≤ToricGB`Private`kdtCutOff,
(* rest of elements as list *)
KDT[[ToricGB`Private`basisKDList]][hint]=F;
,
(* build tree recursively *)
While[True,
(* search feasible coordinate for splitting *)
coord=KDT[[ToricGB`Private`basisProjection,ToricGB`Private`vectorList,level]];
G=Sort[F,(#1[[ToricGB`Private`binomialExponent,coord]]<#2[[ToricGB`Private`binomialExponent,coord]])&];
If[Max[G[[1,ToricGB`Private`binomialExponent,level]],0]≠Max[G[[Length[G],ToricGB`Private`binomialExponent,level]],0],Break[]];
level=Mod[level,Length[KDT[[ToricGB`Private`basisProjection,ToricGB`Private`vectorList]]]]+1;
If[level==levelP,
(* no coordinate found *)
KDT[[ToricGB`Private`basisKDList]][hint]=F;
Return[];
];
];
mid=Quotient[Length[G],2]+1;
value=Max[G[[mid,ToricGB`Private`binomialExponent,coord]],0];
While[mid>1&&value==Max[G[[mid-1,ToricGB`Private`binomialExponent,coord]],0],mid--];
If[mid==1&&value==Max[G[[1,ToricGB`Private`binomialExponent,coord]],0],
While[value==Max[G[[mid,ToricGB`Private`binomialExponent,coord]],0],mid++];
value=Max[G[[mid,ToricGB`Private`binomialExponent,coord]],0];
];
KDT[[ToricGB`Private`basisKDList]][hint]={level,value,Max[G[[1,ToricGB`Private`binomialExponent,coord]],0],Max[G[[mid-1,ToricGB`Private`binomialExponent,coord]],0],value,Max[G[[Length[G],ToricGB`Private`binomialExponent,coord]],0]};
(* lower values *)
it=hint;it[[level,ToricGB`Private`kdtRangeEnd]]=value;
rebuildKDTree[KDT,G[[;;mid-1]],it,Mod[level,Length[KDT[[ToricGB`Private`basisProjection,ToricGB`Private`vectorList]]]]+1];
(* higher values *)
it=hint;it[[level,ToricGB`Private`kdtRangeStart]]=value;
rebuildKDTree[KDT,G[[mid;;]],it,Mod[level,Length[KDT[[ToricGB`Private`basisProjection,ToricGB`Private`vectorList]]]]+1];
];
];

rebuildKDTree[KDT_Symbol,idx_Integer]:=Module[{u,newidx,i},
If[KDT[[ToricGB`Private`basisChanges]]<0.2*Length[KDT[[ToricGB`Private`basisList]]],Return[idx]];
KDT[[ToricGB`Private`basisChanges]]=0;
(*cntRebuild++;
dbgLog["rebuildKDTree[",idx,"]"];
dbgPrint["rebuildKDTree[",idx,"]"];*)
u=KDT[[ToricGB`Private`basisList,idx]];
(*dbgAssert[u[[ToricGB`Private`binomialMarker]]===1,rebuildKDTree,"Binomial with index ",idx," is not active."];*)
KDT[[ToricGB`Private`basisList]]=Select[KDT[[ToricGB`Private`basisList]],(#[[ToricGB`Private`binomialMarker]]===1)&];
newidx=Position[KDT[[ToricGB`Private`basisList]],u,{1},1,Heads->False][[1,1]];
For[i=1,i≤Length[KDT[[ToricGB`Private`basisList]]],i++,KDT[[ToricGB`Private`basisList,i,ToricGB`Private`binomialIndex]]=i];
ClearAll[Evaluate[KDT[[ToricGB`Private`basisKDList]]]];
rebuildKDTree[KDT,KDT[[ToricGB`Private`basisList]],Table[{0,Infinity},{Length[KDT[[ToricGB`Private`basisProjection,ToricGB`Private`vectorList]]]}],1];
newidx
];

SyntaxInformation[insertKDTree]={"ArgumentsPattern"->{_,_,_.}};
SetAttributes[insertKDTree,HoldFirst];
insertKDTree[KDT_Symbol,v_List,hint_List]:=Module[{it,coord,value,level},
(*Print["insertKDTree[",v,",",hint];
Print[KDT[[ToricGB`Private`basisKDList]][hint]];*)
If[KDT[[ToricGB`Private`basisKDList]][hint]==={}||Head[KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`binomialExponent]]]===List,
(* rest of elements as list *)
KDT[[ToricGB`Private`basisKDList]][hint]=Append[KDT[[ToricGB`Private`basisKDList]][hint],v];
,
level=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtIndex]];
coord=KDT[[ToricGB`Private`basisProjection,ToricGB`Private`vectorList,level]];
value=Max[v[[ToricGB`Private`binomialExponent,coord]],0];
it=hint;
If[value>=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtSeparator]],
KDT[[ToricGB`Private`basisKDList]][hint]={KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtIndex]],KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtSeparator]],KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtLowMin]],KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtLowMax]],Min[value,KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtHighMin]]],
Max[value,KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtHighMax]]]};
it[[level,ToricGB`Private`kdtRangeStart]]=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtSeparator]];
,
KDT[[ToricGB`Private`basisKDList]][hint]={KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtIndex]],KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtSeparator]],Min[value,KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtLowMin]]],Max[value,KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtLowMax]]],KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtHighMin]],
KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtHighMax]]};
it[[level,ToricGB`Private`kdtRangeEnd]]=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtSeparator]];
];
insertKDTree[KDT,v,it];
];
]

insertKDTree[KDT_Symbol,v_List]:=Module[{it,pos,coord,u=v,value},
(*cntInsert++;
KDT[[ToricGB`Private`basisChanges]]++;
dbgLog["insertKDTree[",v,"]"];*)
u[[ToricGB`Private`binomialIndex]]=Length[KDT[[ToricGB`Private`basisList]]]+1;
u[[ToricGB`Private`binomialMarker]]=1;
(* append to list of all elements *)
KDT[[ToricGB`Private`basisList]]=Append[KDT[[ToricGB`Private`basisList]],u];
(* insert in real kd-list *)
insertKDTree[KDT,u,Table[{0,Infinity},{Length[KDT[[ToricGB`Private`basisProjection,ToricGB`Private`vectorList]]]}]];
(*Information[Evaluate[KDT[[ToricGB`Private`basisKDList]]]];*)
(*dbgLog["G=",KDT[[ToricGB`Private`basisList]]];*)
]

SyntaxInformation[deleteKDTree]={"ArgumentsPattern"->{_,_,_.}};
SetAttributes[deleteKDTree,HoldFirst];
deleteKDTree[KDT_Symbol,v_List,hint_]:=Module[{it,coord,value,level},
(*dbgAssert[KDT[[ToricGB`Private`basisKDList]][hint]=!={},deleteKDTree,"Could not find ",v," in kd-tree."];*)
If[Head[KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`binomialExponent]]]===List,
(* rest of elements as list *)
KDT[[ToricGB`Private`basisKDList]][hint]=Select[KDT[[ToricGB`Private`basisKDList]][hint],(#[[ToricGB`Private`binomialIndex]]=!=v[[ToricGB`Private`binomialIndex]])&];
,
level=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtIndex]];
coord=KDT[[ToricGB`Private`basisProjection,ToricGB`Private`vectorList,level]];
value=Max[v[[ToricGB`Private`binomialExponent,coord]],0];
it=hint;
If[value>=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtSeparator]],
it[[level,ToricGB`Private`kdtRangeStart]]=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtSeparator]];
,
it[[level,ToricGB`Private`kdtRangeEnd]]=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtSeparator]];
];
deleteKDTree[KDT,v,it];
];
];

deleteKDTree[KDT_Symbol,index_Integer]:=Module[{},
(*cntDelete++;
dbgLog["deleteKDTree[",index,"]"];*)
KDT[[ToricGB`Private`basisChanges]]++;
KDT[[ToricGB`Private`basisList,index,ToricGB`Private`binomialMarker]]=0;
deleteKDTree[KDT,KDT[[ToricGB`Private`basisList,index]],Table[{0,Infinity},{Length[KDT[[ToricGB`Private`basisProjection,ToricGB`Private`vectorList]]]}]];
(*dbgLog["G=",KDT[[ToricGB`Private`basisList]]];*)
];

SyntaxInformation[changeKDTree]={"ArgumentsPattern"->{_,_,_}};
SetAttributes[changeKDTree,HoldFirst];
changeKDTree[KDT_Symbol,index_Integer,v_List]:=Module[{it,pos,coord,u=v,value},
KDT[[ToricGB`Private`basisChanges]]++;
(*cntChange++;
dbgLog["changeKDTree[",index,", ",v,"]"];*)
deleteKDTree[KDT,KDT[[ToricGB`Private`basisList,index]],Table[{0,Infinity},{Length[KDT[[ToricGB`Private`basisProjection,ToricGB`Private`vectorList]]]}]];
u[[ToricGB`Private`binomialIndex]]=index;
u[[ToricGB`Private`binomialMarker]]=1;
KDT[[ToricGB`Private`basisList,index]]=u;
insertKDTree[KDT,u,Table[{0,Infinity},{Length[KDT[[ToricGB`Private`basisProjection,ToricGB`Private`vectorList]]]}]];
(*Information[Evaluate[KDT[[ToricGB`Private`basisKDList]]]];*)
(*dbgLog["G=",KDT[[ToricGB`Private`basisList]]];*)
]

SyntaxInformation[queryKDTree]={"ArgumentsPattern"->{_,_,_,_,_.,_.,_.,_.}};
SetAttributes[queryKDTree,HoldAll];
queryKDTree[KDT_Symbol,query_,cnt_,callback_,inverseQuery_,largeFirst_,hint_,check_]:=Module[{it,coord,pos,return,value,level},
(*cntKDRecursions++;*)
If[KDT[[ToricGB`Private`basisKDList]][hint]==={},Return[]];
If[Head[KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`binomialExponent]]]===List,
(* rest of elements as list *)
For[pos=1,pos≤Length[KDT[[ToricGB`Private`basisKDList]][hint]],pos++,
If[inverseQuery,
return=binomialDividesQ[query,KDT[[ToricGB`Private`basisKDList]][hint][[pos]],{check}],
return=binomialDividesQ[KDT[[ToricGB`Private`basisKDList]][hint][[pos]],query,{check}]
];
If[return,
(*cntKDCallback++;*)
return=callback[KDT,query,KDT[[ToricGB`Private`basisKDList]][hint][[pos]]];
If[return=!=Null,Sow[return];If[--cnt==0,Return[]]];
];
];
,
level=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtIndex]];
coord=KDT[[ToricGB`Private`basisProjection,ToricGB`Private`vectorList,level]];
value=Max[query[[ToricGB`Private`binomialExponent,coord]],0];
If[inverseQuery,
(*cntCompare++;*)
 If[value<=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtHighMax]],
it=hint;it[[level,ToricGB`Private`kdtRangeStart]]=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtSeparator]];
(*cntCompare++;*)
If[value<=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtHighMin]],
queryKDTree[KDT,query,cnt,callback,inverseQuery,largeFirst,it,Complement[check,{coord}]];
If[cnt==0,Return[]];
(*cntCompare++;*)
If[value<=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtLowMax]],
it=hint;it[[level,ToricGB`Private`kdtRangeEnd]]=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtSeparator]];
(*cntCompare++;*)
If[value<=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtLowMin]],
queryKDTree[KDT,query,cnt,callback,inverseQuery,largeFirst,it,Complement[check,{coord}]],
queryKDTree[KDT,query,cnt,callback,inverseQuery,largeFirst,it,check];
];
];
,
queryKDTree[KDT,query,cnt,callback,inverseQuery,largeFirst,it,check];
];
];
,
If[largeFirst,
(*cntCompare++;*)
If[value>=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtHighMin]],
it=hint;it[[level,ToricGB`Private`kdtRangeStart]]=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtSeparator]];
(*cntCompare++;*)
If[value>=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtHighMax]],
queryKDTree[KDT,query,cnt,callback,inverseQuery,largeFirst,it,Complement[check,{coord}]],
queryKDTree[KDT,query,cnt,callback,inverseQuery,largeFirst,it,check];
];
If[cnt==0,Return[]];
it=hint;it[[level,ToricGB`Private`kdtRangeEnd]]=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtSeparator]];
queryKDTree[KDT,query,cnt,callback,inverseQuery,largeFirst,it,Complement[check,{coord}]];
,
(*cntCompare++;*)
If[value>=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtLowMin]],
it=hint;it[[level,ToricGB`Private`kdtRangeEnd]]=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtSeparator]];
(*cntCompare++;*)
If[value>=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtLowMax]],
queryKDTree[KDT,query,cnt,callback,inverseQuery,largeFirst,it,Complement[check,{coord}]],
queryKDTree[KDT,query,cnt,callback,inverseQuery,largeFirst,it,check];
];
];
];
,
(*cntCompare++;*)
If[value>=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtLowMin]],
it=hint;it[[level,ToricGB`Private`kdtRangeEnd]]=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtSeparator]];
(*cntCompare++;*)
If[value>=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtLowMax]],
queryKDTree[KDT,query,cnt,callback,inverseQuery,largeFirst,it,Complement[check,{coord}]];
If[cnt==0,Return[]];
(*cntCompare++;*)
If[value>=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtHighMin]],
it=hint;it[[level,ToricGB`Private`kdtRangeStart]]=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtSeparator]];
(*cntCompare++;*)
If[value>=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtHighMax]],
queryKDTree[KDT,query,cnt,callback,inverseQuery,largeFirst,it,Complement[check,{coord}]],
queryKDTree[KDT,query,cnt,callback,inverseQuery,largeFirst,it,check];
];
];
,
queryKDTree[KDT,query,cnt,callback,inverseQuery,largeFirst,it,check];
];
];
];
];
];
]

queryKDTree[KDT_Symbol,query_,count_,callback_,inverseQuery_:False,largeFirst_:False]:=Module[{cnt=count,results,k},
(*cntQuery++;*)
k=Length[KDT[[ToricGB`Private`basisProjection,ToricGB`Private`vectorList]]];
results=Reap[queryKDTree[KDT,query,cnt,callback,inverseQuery,largeFirst,Table[{0,Infinity},{k}],KDT[[ToricGB`Private`basisProjection,ToricGB`Private`vectorList]]]][[2]];(*dbgLog["queryKDTree[",query,", ",callback,", ",count,", ",inverseQuery,", ",largeFirst,"] = ",results];
If[results==={},cntQueryFail++];*)
results
]

kdt=newKDTree[proj[3,3]];
kdt[[ToricGB`Private`basisList]]=binomialSupport[{{1,2,3},{1,4,6},{4,5,2},{5,3,5},{4,3,2},{6,6,1}}];
rebuildKDTree[kdt,1];
Information[Evaluate[kdt[[ToricGB`Private`basisKDList]]]]

insertKDTree[kdt,binomialSupport[{4,4,4}]];
deleteKDTree[kdt,2];
Information[Evaluate[kdt[[ToricGB`Private`basisKDList]]]]

queryKDTree[kdt,binomialSupport[{4,4,4}],Infinity,#3&]

queryKDTree[kdt,binomialSupport[{1,4,-2}],Infinity,#3&,True]

Mathematica Lists

These are mapping functions which should provide an interface for Mathematica lists exchangeable with kd-lists and kd-trees.

ClearAll[newMathList,clearMathList,rebuildMathList,deleteMathList,changeMathList,queryMathList];

SyntaxInformation[newMathList]={"ArgumentsPattern"->{_}};
newMathList[projection_List]:={{},projection};

SyntaxInformation[clearMathList]={"ArgumentsPattern"->{_}};
SetAttributes[clearMathList,HoldFirst];
clearMathList[ML_Symbol]:=ClearAll[ML];

SyntaxInformation[rebuildMathList]={"ArgumentsPattern"->{_}};
SetAttributes[rebuildMathList,HoldFirst];
rebuildMathList[ML_Symbol,idx_Integer]:=Module[{u,i,newidx},
(*cntRebuild++;
dbgLog["rebuildMathList[",idx,"]"];*)
u=ML[[ToricGB`Private`basisList,idx]];
(*dbgAssert[u[[ToricGB`Private`binomialMarker]]===1,rebuildMathList,"Binomial with index ",idx," is not active."];*)
ML[[ToricGB`Private`basisList]]=Select[ML[[ToricGB`Private`basisList]],(#[[ToricGB`Private`binomialMarker]]===1)&];
newidx=Position[ML[[ToricGB`Private`basisList]],u,{1},1,Heads->False][[1,1]];
For[i=1,i≤Length[ML[[ToricGB`Private`basisList]]],i++,ML[[ToricGB`Private`basisList,i,ToricGB`Private`binomialIndex]]=i];
newidx
];

SyntaxInformation[insertMathList]={"ArgumentsPattern"->{_}};
SetAttributes[insertMathList,HoldFirst];
insertMathList[ML_Symbol,v_List]:=Module[{u=v},
(*cntInsert++;
dbgLog["insertMathList[",v,"]"];*)
u[[ToricGB`Private`binomialIndex]]=Length[ML[[ToricGB`Private`basisList]]]+1;
u[[ToricGB`Private`binomialMarker]]=1;
ML[[ToricGB`Private`basisList]]=Append[ML[[ToricGB`Private`basisList]],u];
(*dbgLog["G=",ML[[ToricGB`Private`basisList]]];*)
];

SyntaxInformation[changeMathList]={"ArgumentsPattern"->{_}};
SetAttributes[changeMathList,HoldFirst];
changeMathList[ML_Symbol,index_Integer,v_List]:=Module[{u=v},
(*cntChange++;
dbgLog["changeMathList[",index,", ",v,"]"];*)
u[[ToricGB`Private`binomialIndex]]=index;
ML[[ToricGB`Private`basisList,index]]=u;
(*dbgLog["G=",ML[[ToricGB`Private`basisList]]];*)
];

SyntaxInformation[deleteMathList]={"ArgumentsPattern"->{_}};
SetAttributes[deleteMathList,HoldFirst];
deleteMathList[ML_Symbol,index_Integer]:=Module[{u=v},
(*cntDelete++;
dbgLog["deleteMathList[",index,"]"];*)
ML[[ToricGB`Private`basisList,index,ToricGB`Private`binomialMarker]]=0;
(*dbgLog["G=",ML[[ToricGB`Private`basisList]]];*)
];

SyntaxInformation[queryMathList]={"ArgumentsPattern"->{_,_,_,_,_.,_.,_.}};
SetAttributes[queryMathList,HoldAll];
queryMathList[ML_Symbol,query_,count_,callback_,inverseQuery_:False,largeFirst_:False]:=Module[{i,results,return,cnt=count},
(*cntQuery++;*)
results=Reap[For[i=1,i≤Length[ML[[ToricGB`Private`basisList]]],i++,
If[ML[[ToricGB`Private`basisList,i,ToricGB`Private`binomialMarker]]=!=1,Continue[]];
If[inverseQuery,
return=binomialDividesQ[query,ML[[ToricGB`Private`basisList,i]],ML[[ToricGB`Private`basisProjection]]],
return=binomialDividesQ[ML[[ToricGB`Private`basisList,i]],query,ML[[ToricGB`Private`basisProjection]]]
];
If[return,
(*cntKDCallback++;*)
return=callback[ML,query,ML[[ToricGB`Private`basisList,i]]];
If[return=!=Null,Sow[return];If[--cnt==0,Break[]]];
];
];
][[2]];
(*dbgLog["queryMathList[",query,", ",callback,", ",count,", ",inverseQuery,", ",largeFirst,"] = ",results];
If[results==={},cntQueryFail++];*)
results
];

Lattices

LatticeByNullSpace

The integer matrix A defines a lattice L = { v integer vector: A.v=0 }.
LatticeByNullSpace computes a basis B (given as list of basis vectors) of L, i.e. A . B[[i]] = 0 for all i.

SyntaxInformation[LatticeByNullSpace]={"ArgumentsPattern"->{_}};
LatticeByNullSpace[A_List]:=Module[{U,H,N},
{U,H} =HermiteDecomposition[Transpose[A]];
N=NullSpace[Transpose[H],Method->DivisionFreeRowReduction];
If[N==={},{},N=N.U;dbgAssert[A.Transpose[N]===Array[0&,{Length[A],Length[N]}],LatticeByNullSpace,"Lattice basis is not orthogonal to A."];N]
];

A={{2,1,1}};
L=LatticeByNullSpace[A]

Attention

NullSpace[A,Method->DivisionFreeRowReduction]

NullSpace doesn't work here. It give a Q-basis, but not a Z-basis. {-1,1,1} = 0.5 * {-1,0,2} + 0.5 *  {-1,2,0}, but there is no integral representation. Thus {{-1, 0, 2}, {-1, 2, 0}} is a basis of a proper sublattice of  L=span{{-1, 0, 2}, {-1, 1, 1}}.

zeroKernelProjection

Find projection wrt. which lattice L (given as lattice basis) is full-dimensional (MatrixRank[L] = Length[projection]) and which has no kernel.

SyntaxInformation[latticeZeroKernelProjection]={"ArgumentsPattern"->{_}};
latticeZeroKernelProjection[L_List]:=Module[{projection={},i,M=Transpose[L],A={},B},
positionsToBitMask[Flatten[Map[(Position[#,a_/;a=!=0,{1},1,Heads->False])&,RowReduce[L]]],Length[M]]
]

latticeZeroKernelProjection[{{1,0,1},{3,0,1}}]

latticeZeroKernelProjection[{{1,1,1},{3,0,1}}]

latticePositiveElement

Finds an small integral vector in the lattice L = { v integer vector: A.v=0 } whose projection is non-negative and whose i-th coordinate is positive.

SyntaxInformation[latticePositiveElement]={"ArgumentsPattern"->{_,_,_}};latticePositiveElement[L_List,projection_List,i_Integer]:=Module[{x,c,b,lu,var,equ},
(* OPT: Use ray algorithms instead of LinearProgramming -> smaller elements, faster? *)
ClearAll[z];
var=Array[z,Length[L]];
equ=var.L;
x=(*addTime[*)FindInstance[Flatten[{equ[[i]]≥1,Thread[equ[[projection[[ToricGB`Private`vectorList]]]]≥0]}],var,Integers](*,timeMinimize]*);
If[x==={},{},binomialSupport[equ/.x[[1]]]]
(*Off[Minimize::infeas];
x=addTime[Minimize[{Table[1,{n}].Abs[var],Apply[And,Thread[A.var==0]]&&var[[i]]≥1&&Apply[And,Thread[var[[Complement[projection[[ToricGB`Private`vectorList]],{i}]]]≥0]]},var,Integers],timeMinimize];
x=addTime[Minimize[{Table[1,{Length[projection[[ToricGB`Private`vectorList]]]}].var[[projection[[ToricGB`Private`vectorList]]]],Apply[And,Thread[A.var==0]]&&var[[i]]≥1&&Apply[And,Thread[var[[Complement[projection[[ToricGB`Private`vectorList]],{i}]]]≥0]]},var,Integers],timeMinimize];
On[Minimize::infeas];
If[x[[1]]===Infinity,Null,binomialSupport[Map[(#[[2]])&,x[[2]]]]]*)
(*c=Table[1,{n}];
b=Table[{0,0},{Length[A]}];
lu=Table[{-Infinity,Infinity},{n}];
lu[[projection[[ToricGB`Private`vectorList]],1]]=0;
lu[[i,1]]=1;
Print[c,A,b,lu,Integers];
x=LinearProgramming[c,A,b,lu,Integers];
Print[x];
If[x[[i]]>0,
x=x*Apply[LCM,Denominator[x]];
dbgAssert[Apply[And,Map[(IntegerQ[#])&,{1,0,2}]]&&A.x===Table[0,{Length[A]}],latticePositiveElement,"x is not an integral element of the lattice."];
binomialSupport[x],
{}]*)
];

latticePositiveElement[LatticeByNullSpace[{{2,1,-1}}],proj[2,3],3]

latticePositiveElement[LatticeByNullSpace[{{2,1,1}}],proj[2,3],3]

latticeMaximalPositiveElements

Finds a set of integral vectors in the lattice L = { v integer vector: A.v=0 } whose sum is non-negative (wrt. projection) with as many positive coordinates as possible.

SyntaxInformation[latticeMaximalPositiveElements]={"ArgumentsPattern"->{_,_}};latticeMaximalPositiveElements[L_List,projection_List]:=Module[{list,y,x,c,b,lu,n=Length[L[[1]]],i,coord},
y=binomialSupport[Table[0,{n}]];
list=Reap[For[i=1,i≤Length[projection[[ToricGB`Private`vectorList]]],i++,
coord=projection[[ToricGB`Private`vectorList,i]];
If[y[[ToricGB`Private`binomialExponent,coord]]>0,Continue[]];
x=latticePositiveElement[L,projection,coord];
If[x=!={},
Sow[x];
y=binomialAdd[y,x];
];
]][[2]];
If[list=!={},list=list[[1]]];
{y,list}
];

latticeMaximalPositiveElements[{{2,1,-1}},proj[3,3]]

maximalGrading

Computer a grading w.r.t. which the vectors in lattice L (given as lattice basis) are homogeneous and which has maximal support.
Therefore find gradings (nonnegative vectors) that are orthogonal to all vectors in L. To obtain maximal support, sum up gradings with different support.

SyntaxInformation[maximalGrading]={"ArgumentsPattern"->{_,_}};
maximalGrading[L_List,projection_List]:=Module[{grading,i,lu,n=Length[L[[1]]]},
grading=Table[0,{n}];
lu=Table[{0,0},{n}];
Map[(lu[[#]]={0,1})&,projection[[ToricGB`Private`vectorList]]];
For[i=1,i<=n,i++,
If[grading[[i]]>0||!MemberQ[projection[[ToricGB`Private`vectorList]],i],Continue[]];
grading+=LinearProgramming[-UnitVector[n,i],L,Table[{0,0},{Length[L]}],lu];
];
grading=grading*Apply[LCM,Denominator[grading]];
dbgAssert[L.grading===Table[0,{Length[L]}],maximalGrading,"L not homogeneous wrt. grading."];{grading,binomialSupportB[grading][[ToricGB`Private`binomialPosSupp]]}
];

maximalGrading[{{-1,0,2},{-1,1,1}},proj[3,3]]

Binomials

ToricGB`KDList=1;
ToricGB`KDTree=2;
ToricGB`MathList=1;

ToricGB`bitMasks=True;
ToricGB`dataStructure=ToricGB`MathList;
ToricGB`useToricGroebnerBasis=True;
ToricGB`intermediateReduce=True;
ToricGB`cascadeIntermediateReduce=True;
ToricGB`alwaysCriterionMinimalSet=True;

In most functions, binomials are stored with some extra information. The format is {v, marker, index, [posSupp, negSupp]}. v ist the exponent vector of the binomial, marker indicates whether it is active (1) or inactive (0). The index is only valid if the binomial is stored in G in geometricBuchberger and refers to the index of the binomial (needed for criterionMinimalSetQ). Otherwise index=0. The last two entries are only created/used in functions ending in B/BQ (stands for bitmask) which are used if ToricGB`bitMasks is set.
A projection is usualy given as list with two elements. The first entry contains the list with the coordinate indices and the second list contains a bitmask of the indices.

Depending on ToricGB`dataStructure, lists of binomials are stored as kd-list, kd-trees or as list of a list of binomials (in all cases, the first two entries are the Mathematica list of binomials and the projection).

ToricGB`Private`binomialExponent=1;
ToricGB`Private`binomialMarker=2;
ToricGB`Private`binomialIndex=3;
ToricGB`Private`binomialPosSupp=4;
ToricGB`Private`binomialNegSupp=5;
ToricGB`Private`vectorList=1;
ToricGB`Private`vectorSupport=2;

Set functions corresponding to flags

If[ToricGB`bitMasks,
binomialSupport:=binomialSupportB;
binomialPositivePart:=binomialPositivePartB;
binomialNegativePart:=binomialNegativePartB;
binomialNegate:=binomialNegateB;
binomialLCM:=binomialLCMB;
binomialAdd:=binomialAddB;
criterionSupportQ:=criterionSupportBQ;
criterionTailQ:=criterionTailBQ;
,
binomialSupport:=binomialSupportN;
binomialPositivePart:=binomialPositivePartN;
binomialNegativePart:=binomialNegativePartN;
binomialNegate:=binomialNegateN;
binomialLCM:=binomialLCMN;
binomialAdd:=binomialAddN;
criterionSupportQ:=criterionSupportNQ;
criterionTailQ:=criterionTailNQ;
];
Switch[ToricGB`dataStructure,ToricGB`KDList,
newList:=newKDList;
clearList:=clearKDList;
rebuildList:=rebuildKDList;
insertList:=insertKDList;
changeList:=changeKDList;
deleteList:=deleteKDList;
queryList:=queryKDList;
binomialReducesBy:=binomialReducesByKD;
binomialReduce:=binomialReduceKD;
criterionMinimalSetQ:=criterionMinimalSetKDQ;

ToricGB`KDTree,
newList:=newKDTree;
clearList:=clearKDTree;
rebuildList:=rebuildKDTree;
insertList:=insertKDTree;
changeList:=changeKDTree;
deleteList:=deleteKDTree;
queryList:=queryKDTree;
binomialReducesBy:=binomialReducesByKD;
binomialReduce:=binomialReduceKD;
criterionMinimalSetQ:=criterionMinimalSetKDQ;

ToricGB`MathList,
newList:=newMathList;
clearList:=clearMathList;
rebuildList:=rebuildMathList;
insertList:=insertMathList;
changeList:=changeMathList;
deleteList:=deleteMathList;
queryList:=queryMathList;
binomialReducesBy:=binomialReducesByML;
binomialReduce:=binomialReduceML;
criterionMinimalSetQ:=criterionMinimalSetMLQ;
];

positionsToBitMask

Given a list of indexes, create a bitmask representing the indexes and return the tuple of the list and the bitmask. The order of the bits has to agree with binomialSupport (see criterionSupportB), i.e. the first digit is the most significant bit.

SyntaxInformation[positionsToBitMask]={"ArgumentsPattern"->{_,_}};
positionsToBitMask[positions_List,n_Integer]:={positions,Total[2^(n-positions)]}

positionsToBitMask[{1,3,4},4]

binomialSupport

Given a vector v, return the triple of v, the marker (set to 1), the product of monomialOrder and v. In the version ending in B append the bitmasks of the support of v+ and the bitmask of the support of v- (N stands for No bitmasks). Also works for lists of vectors.

SyntaxInformation[binomialSupport]={"ArgumentsPattern"->{_}};
binomialSupportN[v_List]:=If[Head[v[[ToricGB`Private`binomialExponent]]]===List,Map[(binomialSupportN[#])&,v],{v,1,0}]
binomialSupportB[v_List]:= If[Head[v[[ToricGB`Private`binomialExponent]]]===List,Map[(binomialSupportB[#])&,v],{v,1,0,FromDigits[Clip[v,{0,1}],2],FromDigits[Clip[-v,{0,1}],2]}]

binomialSupportN[{2,0,-2,1}]

binomialSupportN[{{2,0},{0,-3}}]

binomialSupportB[{2,0,-2,1}]

binomialSupportB[{{2,0},{0,-3}}]

binomialStrip

Given a vector {v,marker,monomialOrder.v,[posSupp,negSupp]}, return v. Also works for lists of vectors.

SyntaxInformation[binomialStrip] = {"ArgumentsPattern" -> {_}};
binomialStrip[v_List] := If[Head[v[[ToricGB`Private`binomialExponent,ToricGB`Private`binomialExponent]]] === List, Map[(binomialStrip[#]) &, v],
  v[[ToricGB`Private`binomialExponent]]]

binomialStrip[binomialSupport[{2,0,-2,1}]]

binomialStrip[binomialSupportB[{{2,0},{0,-3}}]]

binomialPositivePart/binomialNegativePart

Given a vector v, return the vector v+ (resp. v-).

SyntaxInformation[binomialPositivePart]={"ArgumentsPattern"->{_}};
binomialPositivePartN[v_List]:={Clip[v[[ToricGB`Private`binomialExponent]],{0,Infinity}],v[[ToricGB`Private`binomialMarker]],0};binomialPositivePartB[v_List]:={Clip[v[[ToricGB`Private`binomialExponent]],{0,Infinity}],v[[ToricGB`Private`binomialMarker]],0,v[[ToricGB`Private`binomialPosSupp]],0};

SyntaxInformation[binomialNegativePart]={"ArgumentsPattern"->{_}};
binomialNegativePartN[v_List]:={Clip[-v[[ToricGB`Private`binomialExponent]],{0,Infinity}],v[[ToricGB`Private`binomialMarker]],0};binomialNegativePartB[v_List]:={Clip[-v[[ToricGB`Private`binomialExponent]],{0,Infinity}],v[[ToricGB`Private`binomialMarker]],0,v[[ToricGB`Private`binomialNegSupp]],0};

binomialPositivePartN[binomialSupportN[{-3,2,2}]]

binomialNegativePartN[binomialSupportN[{-3,2,2}]]

binomialPositivePartB[binomialSupportB[{-3,2,2}]]

binomialNegativePartB[binomialSupportB[{-3,2,2}]]

binomialNegate

Given v, return -v with corrected extra information.

SyntaxInformation[binomialNegate]={"ArgumentsPattern"->{_}};
binomialNegateN[v_List]:={-v[[ToricGB`Private`binomialExponent]],v[[ToricGB`Private`binomialMarker]],0}
binomialNegateB[v_List]:={-v[[ToricGB`Private`binomialExponent]],v[[ToricGB`Private`binomialMarker]],0,v[[ToricGB`Private`binomialNegSupp]],v[[ToricGB`Private`binomialPosSupp]]}

binomialNegateN[binomialSupportN[{3,-1,0}]]

binomialNegateB[binomialSupportB[{3,-1,0}]]

binomialAdd

Compute v+c*u with corrected extra information.

SyntaxInformation[binomialAdd]={"ArgumentsPattern"->{_,_,_.}};
binomialAddN[v_List,u_List,c_Integer:1]:={v[[ToricGB`Private`binomialExponent]]+c*u[[ToricGB`Private`binomialExponent]],1,0};
binomialAddB[v_List,u_List,c_Integer:1]:=binomialSupportB[v[[ToricGB`Private`binomialExponent]]+c*u[[ToricGB`Private`binomialExponent]]];

binomialAddN[binomialSupportN[{1,2,0}],binomialSupportN[{-1,2,1}]]

binomialAddN[binomialSupportN[{1,2,0}],binomialSupportN[{-1,2,1}],-1]

binomialAddB[binomialSupportB[{1,2,0}],binomialSupportB[{-1,2,1}]]

binomialAddB[binomialSupportB[{1,2,0}],binomialSupportB[{-1,2,1}],-1]

binomialLCM

Computes the least common multiple of (x^v+) and (x^u+) as vector where v+ contains the positive entries of v and v- = (-v)+, i.e. the entrywise maximum of v, u, and 0.

SyntaxInformation[binomialLCM]={"ArgumentsPattern"->{_,_}};
binomialLCMN[v_List,u_List]:=(dbgAssert[Length[v[[ToricGB`Private`binomialExponent]]]===Length[u[[ToricGB`Private`binomialExponent]]],binomialLCM,"Lengths of v and u differ."];binomialPositivePartN[{Array[(Max[v[[ToricGB`Private`binomialExponent,#]],u[[ToricGB`Private`binomialExponent,#]]])&,Length[v[[ToricGB`Private`binomialExponent]]]],1,0}])
binomialLCMB[v_List,u_List]:=(dbgAssert[Length[v[[ToricGB`Private`binomialExponent]]]===Length[u[[ToricGB`Private`binomialExponent]]],binomialLCM,"Lengths of v and u differ."];binomialPositivePartB[{Array[(Max[v[[ToricGB`Private`binomialExponent,#]],u[[ToricGB`Private`binomialExponent,#]]])&,Length[v[[ToricGB`Private`binomialExponent]]]],1,0,BitOr[v[[ToricGB`Private`binomialPosSupp]],u[[ToricGB`Private`binomialPosSupp]]],0}])

binomialLCMN[binomialSupportN[{3,-1,-4}],binomialSupportN[{2,4,-1}]]

binomialLCMB[binomialSupportB[{3,-1,-4}],binomialSupportB[{0,4,-1}]]

binomialZeroQ

Check whether projection (given as list of remaining indeces) of vector v is zero. If v==={} this indicates a zero vector, too.

SyntaxInformation[binomialZeroQ]={"ArgumentsPattern"->{_,_}};
binomialZeroQ[v_,projection_List]:=(v==={}||v[[ToricGB`Private`binomialExponent,projection[[ToricGB`Private`vectorList]]]]===Table[0,{Length[projection[[ToricGB`Private`vectorList]]]}]);

binomialZeroQ[binomialSupport[{0,3,0}],positionsToBitMask[{1,3},3]]

binomialZeroQ[binomialSupport[{0,3,0}],positionsToBitMask[{2,3},3]]

binomialEqualQ

Check whether the projections of the vectors v and u are the same.

SyntaxInformation[binomialEqualQ]={"ArgumentsPattern"->{_,_,_}};
binomialEqualQ[v_List,u_List,projection_List]:=(dbgAssert[Length[v[[ToricGB`Private`binomialExponent]]]===Length[u[[ToricGB`Private`binomialExponent]]],binomialEqualQ,"Lengths of v and u differ."];v[[ToricGB`Private`binomialExponent,projection[[ToricGB`Private`vectorList]]]]===u[[ToricGB`Private`binomialExponent,projection[[ToricGB`Private`vectorList]]]])

binomialEqualQ[binomialSupport[{1,2,3}],binomialSupport[{1,2,4}],proj[2,3]]

binomialDividesQ

Checks whether v+ divides u+ wrt. projection.

SyntaxInformation[binomialDividesQ]={"ArgumentsPattern"->{_,_,_,_.}};
binomialDividesQ[v_List,u_List,projection_List,startidx_Integer:1]:=((*dbgAssert[Length[v[[ToricGB`Private`binomialExponent]]]==Length[u[[ToricGB`Private`binomialExponent]]],binomialDividesQ,"Lengths of v and u differ."];*)
Apply[And,Map[(v[[ToricGB`Private`binomialExponent,#]]≤Max[u[[ToricGB`Private`binomialExponent,#]],0])&,projection[[ToricGB`Private`vectorList]]]])
(*binomialDividesQ[v_List,u_List,projection_List,startidx_Integer:1]:=(binomialDivides[v,u,projection,startidx]>0)*)

binomialDividesQ[binomialSupport[{1,2,-3}],binomialSupport[{1,2,0}],proj[3,3]]

binomialDividesQ[binomialSupport[{1,3,1}],binomialSupport[{1,2,0}],proj[3,3]]

binomialDividesQ[binomialSupport[{2,3,-3}],binomialSupport[{1,-2,0}],proj[3,3],2]

binomialDivides

Computes which power of v+ divides u+ wrt. projection.

SyntaxInformation[binomialDivides]={"ArgumentsPattern"->{_,_,_,_.}};
binomialDivides[v_List,u_List,projection_List,startidx_Integer:1]:=Module[{i,coord,exponent=Infinity},
(*dbgAssert[Length[v[[ToricGB`Private`binomialExponent]]]==Length[u[[ToricGB`Private`binomialExponent]]],binomialDivides,"Lengths of v and u differ."];*)
For[i=startidx,i≤Length[projection[[ToricGB`Private`vectorList]]],i++,
(*cntDividesIter++;*)
coord=projection[[ToricGB`Private`vectorList,i]];
If[v[[ToricGB`Private`binomialExponent,coord]]≤0,Continue[]];
If[v[[ToricGB`Private`binomialExponent,coord]]>u[[ToricGB`Private`binomialExponent,coord]],Return[0]];
exponent=Min[exponent,Quotient[u[[ToricGB`Private`binomialExponent,coord]],v[[ToricGB`Private`binomialExponent,coord]]]];
];
exponent
]

binomialDivides[binomialSupport[{1,2,-3}],binomialSupport[{3,4,0}],proj[3,3]]

binomialDivides[binomialSupport[{1,3,1}],binomialSupport[{1,2,0}],proj[3,3]]

binomialOrder

Returns 1 if v+ > u+ wrt. monomialOrder, given as weight matrix, and projection, -1 if v+ < u+ and 0 otherwise.

SyntaxInformation[binomialOrder]={"ArgumentsPattern"->{_,_,_,_}};
binomialOrder[v_List,u_List,monomialOrder_List,projection_List]:=Module[{i,j,coord,x},
dbgAssert[Length[v[[ToricGB`Private`binomialExponent]]]===Length[u[[ToricGB`Private`binomialExponent]]],binomialOrder,"Lengths of v and u differ."];
For[i=1,i<=Length[monomialOrder],i++,
dbgAssert[Length[v[[ToricGB`Private`binomialExponent]]]===Length[monomialOrder[[i]]],binomialOrder,"monomialOrder has wrong number of columns."];
For[x=0;j=1,j<=Length[projection[[ToricGB`Private`vectorList]]],j++,
coord=projection[[ToricGB`Private`vectorList,j]];
If[v[[ToricGB`Private`binomialExponent,coord]]>0,x+=v[[ToricGB`Private`binomialExponent,coord]]*monomialOrder[[i,coord]]];
If[u[[ToricGB`Private`binomialExponent,coord]]>0,x-=u[[ToricGB`Private`binomialExponent,coord]]*monomialOrder[[i,coord]]];
];
If[x≠0,Return[If[x>0,1,-1]]];
];
0
];

binomialOrder[binomialSupport[{3,5,0}],binomialSupport[{3,4,0}],lexOrder[3],proj[3,3]]

binomialOrder[binomialSupport[{-3,5,-1}],binomialSupport[{-4,5,1}],lexOrder[3],proj[3,3]]

binomialOrder[binomialSupport[{3,5,-5}],binomialSupport[{3,4,-3}],lexOrder[3],positionsToBitMask[{1,3},3]]

binomialStandardFormQ

Returns True if v+ >= v- wrt. monomialOrder and projection.

SyntaxInformation[binomialStandardFormQ]={"ArgumentsPattern"->{_,_,_}};
binomialStandardFormQ[v_List,monomialOrder_List,projection_List]:=(binomialOrder[v,binomialNegate[v],monomialOrder,projection]>=0)

binomialStandardFormQ[binomialSupport[{3,-4,0}],lexOrder[3],proj[3,3]]

binomialStandardFormQ[binomialSupport[{3,-4,0}],lexOrder[3],positionsToBitMask[{2,3},3]]

binomialStandardForm

Returns u in {v, -v} such that isBinomialStandardForm[u], i.e. u+ > u- wrt. monomialOrder and projection.

SyntaxInformation[binomialStandardForm]={"ArgumentsPattern"->{_,_,_}};
binomialStandardForm[v_List,monomialOrder_List,projection_List]:=If[binomialStandardFormQ[v,monomialOrder,projection],v,binomialNegate[v]]

binomialStandardForm[binomialSupport[{3,-4,0}],lexOrder[3],positionsToBitMask[{2,3},3]]

binomialReducesBy

Internal`binomialHead=1;
Internal`binomialTail=2;

Check whether binomial v reduces with respect to binomials in F which have to be in standard form wrt. the desired monomialOrder and projection. monomialSelection controls, whether the head of the binomial, the tail or both are checked. startidx indicates that only the part of F starting with this index has to be checked. The binomials in F contain markers, i.e. are of the form {u,1} if active and {u,0} if inactive. The return value is the index of the first reducer binomial in F or 0 if none is reducable. The sign of the return value indicates whether it reduces v+ (positive) or v- (negative).

SyntaxInformation[binomialReducesBy]={"ArgumentsPattern"->{_,_,_,_,_.,_.}};

binomialReducesByML[F_List,v_List, monomialOrder_List,monomialSelection_Integer,startidx_Integer:1,endidx_Integer:-1]:=Module[{u=binomialNegate[v],i,end=Length[F[[ToricGB`Private`basisList]]]},
If[endidx≥0,end=Min[end,endidx]];
(*dbgLog["FindLessQ[",v,"]"];
cntQuery++;*)
For[i=startidx,i≤end,i++,
If[F[[ToricGB`Private`basisList,i,ToricGB`Private`binomialMarker]]=!=1,Continue[]];
If[BitAnd[monomialSelection,Internal`binomialHead]≠0,
(*cntReduceQDivQ++;*)
If[binomialDividesQ[F[[ToricGB`Private`basisList,i]],v,F[[ToricGB`Private`basisProjection]]],Return[i]]];
If[BitAnd[monomialSelection,Internal`binomialTail]≠0,
(*cntReduceQDivQ++;*)
If[binomialDividesQ[F[[ToricGB`Private`basisList,i]],u,F[[ToricGB`Private`basisProjection]]],Return[-i]]];
];
0
]

SyntaxInformation[callbackReducesByKD]={"ArgumentsPattern"->{_,_,_}};
SetAttributes[callbackReducesByKD,HoldFirst];
callbackReducesByKDL[F_Symbol,query_List,v_List]:=If[F[[ToricGB`Private`basisList,v[[ToricGB`Private`binomialIndex]],ToricGB`Private`binomialMarker]]=!=0,v[[ToricGB`Private`binomialIndex]]]

SetAttributes[binomialReducesByKD,HoldFirst];
binomialReducesByKD[F_Symbol,v_List, monomialOrder_List,monomialSelection_Integer]:=Module[{i},
If[BitAnd[monomialSelection,Internal`binomialHead]≠0,
i=queryList[F,binomialPositivePart[v],1,callbackReducesByKD];
If[i=!={},Return[i[[1,1]]]];
];
If[BitAnd[monomialSelection,Internal`binomialTail]≠0,
i=queryList[F,binomialNegativePart[v],1,callbackReducesByKD];
If[i=!={},Return[-i[[1,1]]]];
];
0
];

binomialReduce

Reduces binomial v with respect to binomials in F which have to be in standard form wrt. the desired monomialOrder and projection. If grading is not Null, apply criterionTailQ. monomialSelection controls, whether the head of the binomial, the tail or both are reduced. startidx indicates that only the part of F starting with this index has to be checked the first time. The binomials in F contain markers, i.e. are of the form {u,1} if active and {u,0} if inactive.

SyntaxInformation[binomialReduce]={"ArgumentsPattern"->{_,_,_,_,_,_.}};

binomialReduceML[F_List,v_List,monomialOrder_List,grading_List,monomialSelection_Integer,startidx_Integer:1]:=Module[{u=v,w=binomialNegate[v],i=startidx,j,c},
(*dbgPrint[ToricGB`Debug`binomialReduce,"binomialReduceML[",F,", ",v,", ",monomialOrder,", ",grading,", ",monomialSelection,", ",startidx];
dbgPrint[ToricGB`Debug`binomialReduce,"FindLess[",u,"]"];*)
For[j=startidx,j≤Length[F[[ToricGB`Private`basisList]]],j++;i=Mod[i,Length[F[[ToricGB`Private`basisList]]]]+1,
If[BitAnd[monomialSelection,Internal`binomialHead]≠0,
(*cntReduceQDivQ++;*)
If[binomialDividesQ[F[[ToricGB`Private`basisList,i]],u,F[[ToricGB`Private`basisProjection]]],
(*cntReduceDiv++;*)
If[grading[[ToricGB`Private`vectorSupport]]≠0&&criterionTailQ[u,F[[ToricGB`Private`basisList,i]],grading,F[[ToricGB`Private`basisProjection]]],(*dbgPrint[ToricGB`Debug`geometricBuchberger,"Drop ",v," during binomialReduce because of criterionTail."];*)Return[{}]];
c=binomialDivides[F[[ToricGB`Private`basisList,i]],u,F[[ToricGB`Private`basisProjection]]];
u=binomialAdd[u,F[[ToricGB`Private`basisList,i]],-c];
If[binomialZeroQ[u,F[[ToricGB`Private`basisProjection]]],Return[{}]];
w=binomialNegate[u];
j=0;
(*dbgPrint[ToricGB`Debug`binomialReduce,"FindLess[",u,"]"];*)
];
];
If[BitAnd[monomialSelection,Internal`binomialTail]≠0,
(*cntReduceQDivQ++;*)
If[binomialDividesQ[F[[ToricGB`Private`basisList,i]],w,F[[ToricGB`Private`basisProjection]]],
(*cntReduceDiv++;*)
If[grading[[ToricGB`Private`vectorSupport]]≠0&&criterionTailQ[w,F[[ToricGB`Private`basisList,i]],grading,F[[ToricGB`Private`basisProjection]]],(*dbgPrint[ToricGB`Debug`geometricBuchberger,"Drop ",v," during binomialReduce because of criterionTail."];*)Return[{}]];
c=binomialDivides[F[[ToricGB`Private`basisList,i]],w,F[[ToricGB`Private`basisProjection]]];
u=binomialAdd[u,F[[ToricGB`Private`basisList,i]],c];
If[binomialZeroQ[u,F[[ToricGB`Private`basisProjection]]],Return[{}]];
w=binomialNegate[u];
j=0;
(*dbgPrint[ToricGB`Debug`binomialReduce,"FindLess[",u,"]"];*)
];
];
];
binomialStandardForm[u,monomialOrder,F[[ToricGB`Private`basisProjection]]]
];

SyntaxInformation[callbackReduceKD]={"ArgumentsPattern"->{_,_,_}};
SetAttributes[callbackReduceKD,HoldFirst];
callbackReduceKD[F_Symbol,query_List,v_List]:=(If[F[[ToricGB`Private`basisList,v[[ToricGB`Private`binomialIndex]],ToricGB`Private`binomialMarker]]=!=0,v])

SyntaxInformation[binomialReduceKD]={"ArgumentsPattern"->{_,_,_,_,_}};
SetAttributes[binomialReduceKD,HoldFirst];
binomialReduceKD[F_Symbol,v_List,monomialOrder_List,grading_List,monomialSelection_Integer]:=Module[{u=v,w,changed=True},
While[changed,
changed=False;
If[BitAnd[monomialSelection,Internal`binomialHead]≠0,
(*dbgPrint[ToricGB`Debug`binomialReduce,"FindLess[",binomialPositivePart[u],"]"];*)
w=queryList[F,binomialPositivePart[u],1,callbackReduceKD,False,True];
If[w=!={},
w=w[[1,1]];
(*dbgPrint[ToricGB`Debug`binomialReduce,"-> ",w];*)
If[grading[[ToricGB`Private`vectorSupport]]≠0&&criterionTailQ[u,w,grading,F[[ToricGB`Private`basisProjection]]],(*dbgPrint[ToricGB`Debug`geometricBuchberger,"Drop ",v," during binomialReduce because of criterionTail."];*)Return[{}]];
u=binomialAdd[u,w,-binomialDivides[w,u,F[[ToricGB`Private`basisProjection]]]];
If[binomialZeroQ[u,F[[ToricGB`Private`basisProjection]]],Return[{}]];
changed=True];
];
If[BitAnd[monomialSelection,Internal`binomialTail]≠0,
(*dbgPrint[ToricGB`Debug`binomialReduce,"FindLess[",binomialNegativePart[u],"]"];*)
w=queryList[F,binomialNegativePart[u],1,callbackReduceKD,False,True];
If[w=!={},
w=w[[1,1]];
(*dbgPrint[ToricGB`Debug`binomialReduce,"-> ",w];*)
If[grading[[ToricGB`Private`vectorSupport]]!=0&&criterionTailQ[binomialNegate[u],w,grading,F[[ToricGB`Private`basisProjection]]],(*dbgPrint[ToricGB`Debug`geometricBuchberger,"Drop ",v," during binomialReduce because of criterionTail."];*)Return[{}]];
u=binomialAdd[u,w,binomialDivides[w,-u,F[[ToricGB`Private`basisProjection]]]];
If[binomialZeroQ[u,F[[ToricGB`Private`basisProjection]]],Return[{}]];
changed=True];
];
];
binomialStandardForm[u,monomialOrder,F[[ToricGB`Private`basisProjection]]]
]

binomialReduceML[{binomialSupport[{{1,-1,-1},{0,1,-1}}],proj[3,3]},binomialSupport[{1,1,0}],lexOrder[3],{{0,0,0},0},3]

F=newKDList[proj[3,3]];
insertKDList[F,binomialSupport[{1,-1,-1}]];
insertKDList[F,binomialSupport[{0,1,-1}]];
binomialReduceKD[F,binomialSupport[{1,1,0}],lexOrder[3],{{0,0,0},0},3]
binomialReducesByKD[F,binomialSupport[{1,1,0}],lexOrder[3],3]

F={binomialSupport[{{1,-1,-1},-{0,1,-1}}],positionsToBitMask[{1,3},3]};
binomialReduceML[F,binomialSupport[{1,1,0}],lexOrder[3],{{0,0,0},0},3]

F=newKDList[proj[1,3]];
insertKDList[F,binomialSupport[{1,-1,-1}]];
insertKDList[F,binomialSupport[-{0,1,-1}]];binomialReduceKD[F,binomialSupport[{1,1,0}],lexOrder[3],{{0,0,0},0},3]

VectorToBinomial

VectorToBinomial converts the vector v to the binomial (x^v+) - (x^v-) where v+ contains the positive entries of v and v- = (-v)+. Also applicable to list of vectors.

SyntaxInformation[VectorToBinomial]={"ArgumentsPattern"->{_,_}};VectorToBinomial[v_List,variables_List]:=If[Head[v[[1]]]===List,Map[(VectorToBinomial[#,variables])&,v],
assert[Length[v]===Length[variables],VectorToBinomial,"Lengths of v and variables differ."];Numerator[Together[Apply[Times,variables^v]-1]]];

VectorToBinomial[{1,2,0},vars[3]]

VectorToBinomial[{{1,2,-4},{3,-6,1},{1,0,-1}},vars[3]]

BinomialToVector

Provides the inverse function to VectorToBinomial. Might inverse the sign, though.

SyntaxInformation[BinomialToVector]={"ArgumentsPattern"->{_,_}};BinomialToVector[F_List,variables_List]:=Map[(BinomialToVector[#,variables])&,F];
BinomialToVector[p_,variables_List]:=Module[{v=Table[0,{Length[variables]}],q,r,c,e},
q=Collect[p,variables];
assert[Head[q]===Plus,BinomialToVector,"Only binomials may be given as input."];
r=First[q]; q=Rest[q];
assert[Head[q]=!=Plus,BinomialToVector,"Only binomials may be given as input."];
e=Exponent[r,variables];
c=Extract[CoefficientList[r,variables],e+1];
assert[c===1||c===-1,BinomialToVector,"All coefficients must be +-1."];
v+=e;
q*=c;
e=Exponent[q,variables];
c=Extract[CoefficientList[q,variables],e+1];
assert[c===-1,BinomialToVector,"The second coefficient of the binomial must have the opposite sign of the first coefficient."];
v-e
];

BinomialToVector[VectorToBinomial[{1,2,0},vars[3]],vars[3]]

BinomialToVector[VectorToBinomial[{{1,2,-4},{3,-6,1},{1,0,-1}},vars[3]],vars[3]]

Criteria for Buchberger's Algorithm

ToricGB`Private`criterionSupport=1;
ToricGB`Private`criterionMinimalSet=2;
ToricGB`Private`criterionTail=4;

criterionSupportQ

A pair {v,u} can be dropped if the intersection of the supports of u+ and v+ (w.r.t. projection) is empty. This is Buchbergers well-known criterion reformulated for vectors (see Hemmecke-Malkin).

SyntaxInformation[criterionSupportQ]={"ArgumentsPattern"->{_,_,_}};criterionSupportNQ[v_List, u_List, projection_List]:=Apply[And,Map[((v[[ToricGB`Private`binomialExponent,#]] <= 0||u[[ToricGB`Private`binomialExponent,#]] <= 0))&,projection[[ToricGB`Private`vectorList]]]]
criterionSupportBQ[v_List, u_List, projection_List]:=(BitAnd[v[[ToricGB`Private`binomialPosSupp]],u[[ToricGB`Private`binomialPosSupp]],projection[[ToricGB`Private`vectorSupport]]]==0)

criterionSupportNQ[binomialSupportN[{0,1,0,1}],binomialSupportN[{4,0,-4,1}],proj[4,4]]

criterionSupportNQ[binomialSupportN[{1,1,0,1}],binomialSupportN[{1,0,0,0}],positionsToBitMask[{2,3,4},4]]

criterionSupportBQ[binomialSupportB[{0,1,0,1}],binomialSupportB[{4,0,-4,1}],proj[4,4]]

criterionSupportBQ[binomialSupportB[{1,1,0,1}],binomialSupportB[{1,0,0,0}],positionsToBitMask[{2,3,4},4]]

criterionMinimalSetQ

It suffices to consider pairs that generate the module of syzygies. criterionMinimalSet drops pairs to get a pretty good approximation. It was originally given by Gebauer and Möller and reformulated to this form by Hemmecke and Malkin. The binomials in F contain markers, i.e. are of the form {u,1} if active and {u,0} if inactive.

SyntaxInformation[criterionMinimalSetQ]={"ArgumentsPattern"->{_,_,_}};

criterionMinimalSetMLQ[F_List,i_Integer,j_Integer]:=Module[{w,k},
(*dbgAssert[i<j,criterionMinimalSetQ,"i must be smaller than j."];*)
w=binomialLCM[F[[ToricGB`Private`basisList,i]],F[[ToricGB`Private`basisList,j]]];
For[k=1,k<i,k++,(*cntMSDivQ++;*)If[F[[ToricGB`Private`basisList,k,ToricGB`Private`binomialMarker]]===1&&binomialDividesQ[F[[ToricGB`Private`basisList,k]],w,F[[ToricGB`Private`basisProjection]]],Return[True]]];
For[k=i+1,k<j,k++,(*cntMSDivQ++;*)If[F[[ToricGB`Private`basisList,k,ToricGB`Private`binomialMarker]]===1&&binomialDividesQ[F[[ToricGB`Private`basisList,k]],w,F[[ToricGB`Private`basisProjection]]]&&!binomialEqualQ[binomialLCM[F[[ToricGB`Private`basisList,k]],F[[ToricGB`Private`basisList,j]]],w,F[[ToricGB`Private`basisProjection]]],Return[True]]];
For[k=j+1,k≤Length[F[[ToricGB`Private`basisList]]],k++,(*cntMSDivQ++;*)If[F[[ToricGB`Private`basisList,k,ToricGB`Private`binomialMarker]]===1&&binomialDividesQ[F[[ToricGB`Private`basisList,k]],w,F[[ToricGB`Private`basisProjection]]]&&!binomialEqualQ[binomialLCM[F[[ToricGB`Private`basisList,k]],F[[ToricGB`Private`basisList,i]]],w,F[[ToricGB`Private`basisProjection]]]&&!binomialEqualQ[binomialLCM[F[[ToricGB`Private`basisList,k]],F[[ToricGB`Private`basisList,j]]],w,F[[ToricGB`Private`basisProjection]]],Return[True]]];
False
];

SetAttributes[callbackMinimalSetKD,HoldFirst];
callbackMinimalSetKD[F_Symbol,query_List,v_List,i_Integer,j_Integer]:=
(If[F[[ToricGB`Private`basisList,v[[ToricGB`Private`binomialIndex]],ToricGB`Private`binomialMarker]]===1,
If[v[[ToricGB`Private`binomialIndex]]<i,True,
If[i<v[[ToricGB`Private`binomialIndex]]&&v[[ToricGB`Private`binomialIndex]]<j,If[!binomialEqualQ[binomialLCM[v,F[[ToricGB`Private`basisList]][[j]]],query,F[[ToricGB`Private`basisProjection]]],True],
If[j<v[[ToricGB`Private`binomialIndex]],If[!binomialEqualQ[binomialLCM[v,F[[ToricGB`Private`basisList]][[i]]],query,F[[ToricGB`Private`basisProjection]]]&&!binomialEqualQ[binomialLCM[v,F[[ToricGB`Private`basisList]][[j]]],query,F[[ToricGB`Private`basisProjection]]],True]
]
]
]
])

SetAttributes[criterionMinimalSetKDQ,HoldFirst];
criterionMinimalSetKDQ[F_Symbol,i_Integer,j_Integer]:=Module[{w},
(*dbgAssert[i<j,criterionMinimalSetQ,"i must be smaller than j."];*)
w=binomialLCM[F[[ToricGB`Private`basisList]][[i]],F[[ToricGB`Private`basisList]][[j]]];
queryList[F,w,1,Function[{KD,query,it},callbackMinimalSetKD[KD,query,it,i,j],{HoldAll}]]=!={}
];

criterionTailQ

This criterion is specific to homogeneous ideals (w.r.t. some grading) that are saturated w.r.t. the product of the variables. If the S-polynomial spol reduces to some polynomial recuct of lower grading, the S-polynomial can be dropped. When called by geometricBuchberger, it checks, whether the reduction of the actual S-polynomial to v-u already is a decrease of the grading. This function is also called by monomialReduce. Here it checks, whether the reduction of u by v to u-v (or vice versa) decreases the grading. This criterion is special to toric ideals and was given by Hemmecke and Malkin.

SyntaxInformation[criterionTailQ]={"ArgumentsPattern"->{_,_,_,_}};
criterionTailNQ[u_List,v_List,grading_List,projection_List]:=(dbgAssert[Length[v[[ToricGB`Private`binomialExponent]]]===Length[u[[ToricGB`Private`binomialExponent]]],criterionTailQ,"Lengths of v and u differ."];
Apply[Or,Map[((v[[ToricGB`Private`binomialExponent,#]]<0&&u[[ToricGB`Private`binomialExponent,#]]<0&&grading[[ToricGB`Private`vectorList,#]]>0))&,projection[[ToricGB`Private`vectorList]]]])
criterionTailBQ[u_List,v_List,grading_List,projection_List]:=(dbgAssert[Length[v[[ToricGB`Private`binomialExponent]]]===Length[u[[ToricGB`Private`binomialExponent]]],criterionTailQ,"Lengths of v and u differ."];
BitAnd[v[[ToricGB`Private`binomialNegSupp]],u[[ToricGB`Private`binomialNegSupp]],grading[[ToricGB`Private`vectorSupport]],projection[[ToricGB`Private`vectorSupport]]]≠0)

grad1={{1,2,0},binomialSupportB[{1,2,0}][[ToricGB`Private`binomialPosSupp]]};
criterionTailNQ[binomialSupportN[{-1,0,-1}],binomialSupportN[{-2,0,0}],grad1,proj[3,3]]

grad2={{0,1,1},binomialSupportB[{0,1,1}][[ToricGB`Private`binomialPosSupp]]};criterionTailNQ[binomialSupportN[{-1,-1,0}],binomialSupportN[{-1,0,-1}],grad2,proj[3,3]]

criterionTailBQ[binomialSupportB[{-1,0,-1}],binomialSupportB[{-2,0,0}],grad1,proj[3,3]]

criterionTailBQ[binomialSupportB[{-1,-1,0}],binomialSupportB[{-1,0,-1}],grad2,proj[3,3]]

Algorithms

intermediateReduce

During the Buchberger algorithm sometimes one gets cascading sequences of polynomials where the later ones reduce the previous ones. In order to keep the intermediate bases small, an interreduction is performed. In this implementation, intermediateReduce is called for each new polynomial. Instead of reducing each polynomial wrt. all other polynomials, a quick scan is performed for polynomials that are reducible by the new polynomial. For those a full reduction is done. Then (if ToricGB`cascadeIntermediateReduce is set) those are treated as new and therefore polynomials which are reducible by them are searched.

SetAttributes[callbackIntermediateReduce,HoldFirst];
callbackIntermediateReduce[F_Symbol,query_List,v_List]:=(If[F[[ToricGB`Private`basisList,v[[ToricGB`Private`binomialIndex]],ToricGB`Private`binomialMarker]]=!=0,v[[ToricGB`Private`binomialIndex]]]);

SyntaxInformation[intermediateReduce]={"ArgumentsPattern"->{_,_,_,_}};
SetAttributes[intermediateReduce,HoldFirst];
intermediateReduce[F_Symbol,min_Integer,new_Integer,monomialOrder_List]:=Module[{changed=newQueue[],results,i,j,u,v,e,cur,nullGrading={Table[0,{Length[F[[ToricGB`Private`basisList,1]]]}],0}},
(*dbgPrint[ToricGB`Debug`geometricBuchberger,"intermediateReduce: F=",F[[ToricGB`Private`basisList]],", min=",min,", new=",new,", monomialOrder=",monomialOrder,", projection=",F[[ToricGB`Private`basisProjection]]];*)
enQueue[changed,{new,F[[ToricGB`Private`basisList,new]]}];
While[!queueEmpty[changed],
{i,u}=deQueue[changed];
If[F[[ToricGB`Private`basisList,i]]=!=u,Continue[]];
results=queryList[F,u,Infinity,callbackIntermediateReduce,True];
If[results=!={},results=results[[1]]];
For[j=1,j≤Length[results],j++,
If[results[[j]]==i,Continue[]];
(*cntReduceDiv++;*)
e=binomialDivides[u,F[[ToricGB`Private`basisList,results[[j]]]],F[[ToricGB`Private`basisProjection]]];
(*dbgAssert[e>0,intermediateReduce,"Query returned element ",F[[ToricGB`Private`basisList,results[[j]]]]," not divivisible by ",u,". F=",F[[ToricGB`Private`basisList]],", results=",results];*)
v=binomialAdd[F[[ToricGB`Private`basisList,results[[j]]]],u,-e];
F[[ToricGB`Private`basisList,results[[j]],ToricGB`Private`binomialMarker]]=0;
v=binomialReduce[F,v,monomialOrder,nullGrading,BitOr[Internal`binomialHead,Internal`binomialTail]];
(*dbgPrint[ToricGB`Debug`geometricBuchberger,F[[ToricGB`Private`basisList,results[[j]]]]," reduces to ",v];*)
If[!binomialZeroQ[v,F[[ToricGB`Private`basisProjection]]],
If[results[[j]]≤min,
deleteList[F,results[j]];
insertList[F,v];
If[ToricGB`cascadeIntermediateReduce,enQueue[changed,{Length[F[[ToricGB`Private`basisList]]],F[[ToricGB`Private`basisList,Length[F[[ToricGB`Private`basisList]]]]]}]];
,
changeList[F,results[[j]],v];
If[ToricGB`cascadeIntermediateReduce,enQueue[changed,{results[[j]],F[[ToricGB`Private`basisList,results[[j]]]]}]];
];
];
];
];
clearQueue[changed];
];

geometricBuchberger

Performs the Buchberger algorithm for the toric ideal I=VectorToBinomial[F] doing all computations with vectors (of the exponents). The result is a Gröbner basis (in vector form) of an ideal containing I and contained in the radical of I.
In the main loop, for all pairs {v,u} in the approximation to the Gröbner basis G the S-polynomial S(v,u)=v-u is calculated and reduced w.r.t. G. If the remainder is non-zero, it is added to G. Different criteria tell which pairs {v,u} may be omitted.
There are several criteria for dropping pairs which can be controlled by the bitmask criteria.
Finally the Gröbner basis is reduced. Head-reducible binomials are superfluous. Intermediate reductions are applied at the very beginning and on insertion of new elements if ToricGB`intermediateReduce is set.

SyntaxInformation[geometricBuchberger]={"ArgumentsPattern"->{_,_,_,_,_}};geometricBuchberger[F_List, monomialOrder_List, gradingP_List,projection_List,criteria_Integer] :=Module[{G=newList[projection],H,spol,i, j,k,l,v,u,grading=gradingP,n=Length[F[[1]]]},
If[BitAnd[criteria,ToricGB`Private`criterionTail]==0,grading={Table[0,{n}],0}];
(*dbgPrint[ToricGB`Debug`geometricBuchberger,"geometricBuchberger[",F,", ",monomialOrder,", ",grading,", ",projection,", ",criteria,"]"];
dbgLog["geometricBuchberger[",F,", ",monomialOrder,", ",grading,", ",projection,", ",criteria,"]"];*)
For[i=1,i<=Length[F],i++,
(*cntInsert++;*)
insertList[G,binomialStandardForm[F[[i]],monomialOrder,projection]];
];
(*timePrereduce+=AbsoluteTiming[*)For[i=1,i<=Length[G[[ToricGB`Private`basisList]]],i++,
If[G[[ToricGB`Private`basisList,i,ToricGB`Private`binomialMarker]]===0,Continue[]];
(*cntReduce0++;*)
intermediateReduce[G,0,i,monomialOrder];
](*][[1]]*);
(*dbgPrint[ToricGB`Debug`geometricBuchberger,"Prereduced: G=",G[[ToricGB`Private`basisList]]];*)
(*timeGroebner+=AbsoluteTiming[*)
For[j=1,j≤Length[G[[ToricGB`Private`basisList]]],j++,
If[G[[ToricGB`Private`basisList,j,ToricGB`Private`binomialMarker]]===0,Continue[]];
j=(*addTime[*)rebuildList[G,j](*,timeRebuild]*);
For[i=1,i<j,i++,
If[G[[ToricGB`Private`basisList,i,ToricGB`Private`binomialMarker]]===0,Continue[]];
v=G[[ToricGB`Private`basisList,i]]; u=G[[ToricGB`Private`basisList,j]];
(*cntCS++;*)
If[BitAnd[criteria,ToricGB`Private`criterionSupport]≠0&&criterionSupportQ[v,u,projection],(*dbgPrint[ToricGB`Debug`geometricBuchberger,"Drop ",{i,j}," because of criterionSupportQ."];*)Continue[]];
(*cntCT++;*)
If[BitAnd[criteria,ToricGB`Private`criterionTail]≠0&&criterionTailQ[v,u,grading,projection],(*dbgPrint[ToricGB`Debug`geometricBuchberger,"Drop ",{i,j}," because of criterionTailQ."];*)Continue[]];
(*cntCM++;*)
If[BitAnd[criteria,ToricGB`Private`criterionMinimalSet]≠0&&(*addTime[*)criterionMinimalSetQ[G,i,j](*,timeCM]*),(*dbgPrint[ToricGB`Debug`geometricBuchberger,"Drop ",{i,j}," because of criterionMinimalSetQ."];*)Continue[]];
spol=binomialAdd[v,u,-1];
(*cntReduce1++;*)
spol=binomialReduce[G,spol,monomialOrder,grading,BitOr[Internal`binomialHead,Internal`binomialTail]];
If[!binomialZeroQ[spol,projection],
(*cntInsert++;*)
insertList[G,spol];
If[ToricGB`intermediateReduce,(*addTime[*)intermediateReduce[G,j+1,Length[G[[ToricGB`Private`basisList]]],monomialOrder](*,timeInterRed]*)];
];
(*dbgPrint[ToricGB`Debug`geometricBuchberger,"Pair ",{i,j},": ",binomialAdd[v,u,-1]," reduces to ",spol];
dbgPrint[ToricGB`Debug`geometricBuchberger,"G=",G[[ToricGB`Private`basisList]]];*)
];
](*][[1]]*);
(* reduce Gröbner basis *)
(*dbgPrint[ToricGB`Debug`geometricBuchberger,"Reduce... G=",G[[ToricGB`Private`basisList]]];*)
(*timeRedGB+=AbsoluteTiming[*)H=Reap [For[i=1,i≤Length[G[[ToricGB`Private`basisList]]],i++,
If[G[[ToricGB`Private`basisList,i,ToricGB`Private`binomialMarker]]===0,Continue[]];
v=G[[ToricGB`Private`basisList,i]];
G[[ToricGB`Private`basisList,i,ToricGB`Private`binomialMarker]]=0;
If[!ToricGB`intermediateReduce||!ToricGB`cascadeIntermediateReduce,
(*cntReduce2Q++;*)
If[binomialReducesBy[G,v,monomialOrder,Internal`binomialHead]≠0,
(*dbgPrint[ToricGB`Debug`geometricBuchberger,"Drop G[[",i,"]]=",v," because of binomialReduceQ"];*)
Continue[]];
];
(*cntReduce2++;*)
v=binomialReduce[G,v,monomialOrder,{Table[0,{n}],0},Internal`binomialTail];
G[[ToricGB`Private`basisList,i,ToricGB`Private`binomialMarker]]=1;
(*dbgPrint[ToricGB`Debug`geometricBuchberger,G[[ToricGB`Private`basisList,i]]," reduces to ",v];*)
If[binomialZeroQ[v,projection],
deleteList[G,i,v];
,
If[v[[ToricGB`Private`binomialExponent]]!=G[[ToricGB`Private`basisList,i,ToricGB`Private`binomialExponent]],
v[[ToricGB`Private`binomialMarker]]=1;
changeList[G,i,v];
];
Sow[v];
];
]][[2,1]](*][[1]]*);
clearList[G];
(*dbgAssert[dbgGroebnerBasisQ[H,monomialOrder,projection],geometricBuchberger,"G should be a Gröbner basis now. F=",F, ", monomialOrder=",monomialOrder,", grading=",gradingP,", projection=",projection,", criteria=",criteria];
dbgPrint[ToricGB`Debug`geometricBuchberger,"G=",H];*)
H
];

F = {{1, 2, -3}, {0, -2, 4}, {-2, 3, 1}}; dbgGroebnerBasisUnique[VectorToBinomial[binomialStrip[geometricBuchberger[binomialSupport[F], lexOrder[3],{{0,0,0},0},proj[3,3],3]], vars[3]],vars[3]]

dbgGroebnerBasisUnique[GroebnerBasis[VectorToBinomial[F, vars[3]], vars[3]],vars[3]]

Attention

The is not result equal since geometricBuchberger constructs a Gröbner basis of an ideal I < J < Sqrt[I].

F={{13,-64,0,29},{-5,24,1,-10}};
order={{8,(13/8),0,0},{1,0,0,0},{0,1,0,0},{0,0,1,0},{0,0,0,1}};
dbgGroebnerBasisUnique[VectorToBinomial[binomialStrip[geometricBuchberger[binomialSupport[F], order,{{0,0,0,0},0},proj[2,4],7]],vars[4]],vars[4]]

profile[geometricBuchberger[binomialSupport[{{1,0,38,9,-19,13},{0,1,366,78,-184,126},{0,0,1123,239,-565,386},{239,0,98,239,-21,19},{1003,1025,0,1025,263,141},{471,544,263,544,0,169}}],lexOrder[6],{{0,0,0,0,0,0},0},proj[6,6],3]]

VEC, bitmasks, geometricBuchberger, profiling, LatticeReduce, part intermediateReduce

VEC, bitmasks, geometricBuchberger, profiling, LatticeReduce, intermediateReduce

profile[geometricBuchberger[binomialSupport[{{1,1,-1,1,-2},{0,10,1,21,-46},{0,0,3,0,2},{2,2,1,2,-2},{42,22,7,0,42}}],lexOrder[5],{{0,0,0,0,0},0},proj[5,5],3]]

VEC, bitmasks, geometricBuchberger, profiling, LatticeReduce, part intermediateReduce

VEC, bitmasks, geometricBuchberger, profiling, LatticeReduce, intermediateReduce

projectionOrdering

Computes an ordering such that the Gröbner basis G with respect to this ordering and projection is also a Gröbner basis with respect to this ordering and Sort[Append[projection,i]].

SyntaxInformation[projectionOrdering]={"ArgumentsPattern"->{_,_,_,_}};
projectionOrdering[L_List,projection_List,i_Integer,monomialOrder_List]:=Module[{x,c,A,b,n=Length[L[[1]]],projectionComplement},
c=Table[1,{n}];
projectionComplement=Complement[Range[n],projection[[ToricGB`Private`vectorList]]];
A=FlattenAt[Append[L,Map[(UnitVector[n,#])&,projectionComplement]],Length[L]+1];
b=-Array[({L[[#,i]],0})&,Length[L]];
b=FlattenAt[Append[b,Table[{0,0},{Length[projectionComplement]}]],Length[b]+1];
x=LinearProgramming[c,A,b];
x=x*Apply[LCM,Denominator[x]];
Prepend[monomialOrder,x]
];

L={{1,-5,-3,3},{0,8,4,-5}};
order=projectionOrdering[L,positionsToBitMask[{1,2,3},4],4,{}]

L.order[[1]]

L={{1,6,-7,-4,2,-15},{0,7,2,-8,9,-19},{0,0,23,-10,19,-7}};
order=projectionOrdering[L,positionsToBitMask[{1,2,3,4,5},6],6,{}]

L.order[[1]]

HemmeckeToricGroebnerBasis

The integer matrix A defines a lattice L = { v integer vector: A.v=0 }. (see also LatticeNullSpace)
The corresponding monomial ideal is defined by I = < (x^v+) - (x^v-) : v in L >.
HemmeckeToricGroebnerBasis computes a Gröbner basis G (given as list of polynomials) of I .

LatticeReduce is used as preprocessing such that the vectors have small norms. Then a projection is found such that the lift of this projection is unique and we easily can compute a generating set (aka Markov basis). For this projection there always exists a totally positive element u. Adding this to the basis G ensures that it is a generating set. Now we compute the unbounded variables and lift the generating set wrt. them. Finally we lift wrt. the bounded variables and compute a Gröbner basis.

SyntaxInformation[HemmeckeToricGroebnerBasis]={"ArgumentsPattern"->{_,_,_}};
HemmeckeToricGroebnerBasis[B_List,variables_List,monomialOrder_List]:=Module[{L,G,grading,projection,i,u,v,list,time,H,criteria,order,n=Length[B[[1]]],unbounded,A},
assert[Length[variables]==n,HemmeckeToricGroebnerBasis,"The number of variables provided is wrong."];
A=LatticeReduce[B];
L=LatticeByNullSpace[A];
If[L==={},Return[{}]];
L=LatticeReduce[L];
(*dbgPrint[ToricGB`Debug`HemmeckeToricGroebnerBasis,"HemmeckeToricGroebnerBasis: A=",A,", L=",L, "monomialOrder=",monomialOrder];*)
G=binomialSupport[L];
(* find projection and suitable generating set *)
projection=latticeZeroKernelProjection[L];
{u,list}=latticeMaximalPositiveElements[L,projection]; (* Could be done faster using Integer-Gauß since projection of L is full rank *)
(*assert[Apply[And,Thread[u[[ToricGB`Private`binomialExponent,projection[[ToricGB`Private`vectorList]]]]>0]],HemmeckeToricGroebnerBasis,"No totally positive element in projection of lattice found."];
dbgPrint[ToricGB`Debug`HemmeckeToricGroebnerBasis,"Initial projection=",projection,", G=",G,", u=",u];*)
G=Join[G,list];
(* find globally bounded / unbounded variables *)
{u,list}=latticeMaximalPositiveElements[L,proj[n,n]];
unbounded=Select[Complement[Table[i,{i,1,n}],projection[[ToricGB`Private`vectorList]]],(u[[ToricGB`Private`binomialExponent,#]]≤0)&];
(*dbgPrint[ToricGB`Debug`HemmeckeToricGroebnerBasis,"Unbounded variables=",unbounded,", u=",u];*)
(* lift generating set from projection to all unbounded variables *)
While[Length[unbounded]>0,
i=First[unbounded];unbounded=Rest[unbounded];
v=latticePositiveElement[L,projection,i];
(*dbgPrint[ToricGB`Debug`HemmeckeToricGroebnerBasis,"G=",G,", projection=", projection,", coordinate=",i,": latticePositiveElement=",u];*)
If[v=!={},
If[!MemberQ[binomialStrip[G],binomialStrip[v]],G=Append[G,v]],
grading=maximalGrading[L,projection];
criteria=ToricGB`Private`criterionSupport;
If[ToricGB`alwaysCriterionMinimalSet||Count[grading[[ToricGB`Private`vectorList]],x_/;x>0]<=Length[variables]/2,criteria=BitOr[criteria,ToricGB`Private`criterionMinimalSet]];
If[grading[[ToricGB`Private`vectorSupport]]≠0,criteria=BitOr[criteria,ToricGB`Private`criterionTail]];
{time,G}=AbsoluteTiming[geometricBuchberger[G,projectionOrdering[L,projection,i,monomialOrder],grading,
projection,criteria]];
(*dbgPrint[ToricGB`Debug`HemmeckeToricGroebnerBasis,"geometricBuchberger: G=",G,", order=",projectionOrdering[L,projection,i,monomialOrder],", grading=",grading,", projection=",projection,", criteria=",criteria," : time = ",time];
dbgAssert[dbgGroebnerBasisQ[G,projectionOrdering[L,projection,i,monomialOrder],positionsToBitMask[Sort[Append[projection[[ToricGB`Private`vectorList]],i]],n]],HemmeckeToricGroebnerBasis,"G should be a Gröbner basis now."];*)
];
projection=positionsToBitMask[Sort[Append[projection[[ToricGB`Private`vectorList]],i]],n];
];
(* lift generating set to bounded variables *)
If[!binomialZeroQ[u,projection],
G=Join[G,list];
(*dbgPrint[ToricGB`Debug`HemmeckeToricGroebnerBasis,"Final lift with u=",u];*)
];
(* compute Gröbner basis *)
projection=proj[n,n];
grading=maximalGrading[L,projection];
criteria=ToricGB`Private`criterionSupport;
If[ToricGB`alwaysCriterionMinimalSet||Count[grading[[ToricGB`Private`vectorList]],x_/;x>0]<=Length[variables]/2,criteria=BitOr[criteria,ToricGB`Private`criterionMinimalSet]];
If[grading[[ToricGB`Private`vectorSupport]]≠0,criteria=BitOr[criteria,ToricGB`Private`criterionTail]];
(*dbgAssert[dbgSaturatedQ[G],HemmeckeToricGroebnerBasis,"Ideal should be saturated by now. (A=",A,")"];*)
If[ToricGB`useToricGroebnerBasis,
{time,H}=AbsoluteTiming[GroebnerBasis`ToricGroebnerBasis[VectorToBinomial[binomialStrip[G],variables],variables,MonomialOrder->monomialOrder]];
Print[time];
(*dbgPrint[ToricGB`Debug`HemmeckeToricGroebnerBasis,"GroebnerBasis`ToricGroebnerBasis: time=",time];*)
H,
(*dbgPrint[ToricGB`Debug`HemmeckeToricGroebnerBasis,"geometricBuchberger: G=",G,", order=",monomialOrder,", grading=", grading,", projection=",projection,", criteria=",criteria," : time = ",time];*)
(*{time,H}=AbsoluteTiming[geometricBuchberger[G,monomialOrder,grading,projection,criteria]];*)
H=geometricBuchberger[G,monomialOrder,grading,projection,criteria];
(*dbgPrint[ToricGB`Debug`HemmeckeToricGroebnerBasis,"geometricBuchberger: G=",G,", order=",monomialOrder,", grading=", grading,", projection=",projection,", criteria=",criteria," : time = ",time];*)
VectorToBinomial[binomialStrip[H],variables]]
];

Example 1

B[1]={{2,1,1}};
dbgGroebnerBasisUnique[HemmeckeToricGroebnerBasis[B[1],vars[3],lexOrder[3]],vars[3]]

dbgGroebnerBasisUnique[GroebnerBasis`ToricGroebnerBasis[Append[VectorToBinomial[LatticeByNullSpace[B[1]],vars[3]],1-t*Apply[Times,vars[3]]],vars[3],t],vars[3]]

Example 2

B[2]={{-6,-2,9,2},{-1,3,-10,1}};
dbgGroebnerBasisUnique[HemmeckeToricGroebnerBasis[B[2],vars[4],lexOrder[4]],vars[4]]

dbgGroebnerBasisUnique[
GroebnerBasis`ToricGroebnerBasis[Append[VectorToBinomial[LatticeByNullSpace[B[2]],vars[4]],1-t*Apply[Times,vars[4]]],vars[4],t],vars[4]]

Example 3

B[3]={{-3,8,-5,-8}};
dbgGroebnerBasisUnique[HemmeckeToricGroebnerBasis[B[3],vars[4],lexOrder[4]],vars[4]]

dbgGroebnerBasisUnique[GroebnerBasis`ToricGroebnerBasis[Append[VectorToBinomial[LatticeByNullSpace[B[3]],vars[4]],1-t*Apply[Times,vars[4]]],vars[4],t],vars[4]]

Example 4

B[4]={{-3,-10,-8}};
dbgGroebnerBasisUnique[HemmeckeToricGroebnerBasis[B[4],vars[3],lexOrder[3]],vars[3]]

dbgGroebnerBasisUnique[GroebnerBasis`ToricGroebnerBasis[Append[VectorToBinomial[LatticeByNullSpace[B[4]],vars[3]],1-t*Apply[Times,vars[3]]],vars[3],t],vars[3]]

Example 5

B[5]={{-9,-1,3,-5,-5,10},{5,0,-2,3,5,1}};
dbgGroebnerBasisUnique[HemmeckeToricGroebnerBasis[B[5],vars[6],lexOrder[6]],vars[6]]

dbgGroebnerBasisUnique[GroebnerBasis`ToricGroebnerBasis[Append[VectorToBinomial[LatticeByNullSpace[B[5]],vars[6]],1-t*Apply[Times,vars[6]]],vars[6],t],vars[6]]

Example 6

B[6]={{2,-9,-4}};
dbgGroebnerBasisUnique[HemmeckeToricGroebnerBasis[B[6],vars[3],lexOrder[3]],vars[3]]

dbgGroebnerBasisUnique[GroebnerBasis`ToricGroebnerBasis[Append[VectorToBinomial[LatticeByNullSpace[B[6]],vars[3]],1-t*Apply[Times,vars[3]]],vars[3],t],vars[3]]

Example 7

B[7]={{-1,3,-5,-5,10}};
dbgGroebnerBasisUnique[HemmeckeToricGroebnerBasis[B[7],vars[5],lexOrder[5]],vars[5]]

dbgGroebnerBasisUnique[GroebnerBasis`ToricGroebnerBasis[Append[VectorToBinomial[LatticeByNullSpace[B[7]],vars[5]],1-t*Apply[Times,vars[5]]],vars[5],t],vars[5]]

Example 8

B[8]={{3,9,7,-4},{6,2,7,7}};
dbgGroebnerBasisUnique[HemmeckeToricGroebnerBasis[B[8],vars[4],lexOrder[4]],vars[4]]

dbgGroebnerBasisUnique[
GroebnerBasis`ToricGroebnerBasis[Append[VectorToBinomial[LatticeByNullSpace[B[8]],vars[4]],1-t*Apply[Times,vars[4]]],vars[4],t],vars[4]]

Example 9

B[9]={{-10,5,-9,-8,10,9},{-6,-2,-4,-3,4,2},{-3,3,-3,5,7,2}};
dbgGroebnerBasisUnique[HemmeckeToricGroebnerBasis[B[9],vars[6],lexOrder[6]],vars[6]]

dbgGroebnerBasisUnique[GroebnerBasis`ToricGroebnerBasis[Append[VectorToBinomial[LatticeByNullSpace[B[9]],vars[6]],1-t*Apply[Times,vars[6]]],vars[6],t],vars[6]]

Example 10

B[10]={{1,6,-5,-5},{0,-6,-10,-9}};
dbgGroebnerBasisUnique[HemmeckeToricGroebnerBasis[B[10],vars[4],lexOrder[4]],vars[4]]

dbgGroebnerBasisUnique[GroebnerBasis`ToricGroebnerBasis[Append[VectorToBinomial[LatticeByNullSpace[B[10]],vars[4]],1-t*Apply[Times,vars[4]]],vars[4],t],vars[4]]

Example 11

B[11]={{-10,9,9}};
dbgGroebnerBasisUnique[HemmeckeToricGroebnerBasis[B[11],vars[3],lexOrder[3]],vars[3]]

dbgGroebnerBasisUnique[GroebnerBasis`ToricGroebnerBasis[Append[VectorToBinomial[LatticeByNullSpace[B[11]],vars[3]],1-t*Apply[Times,vars[3]]],vars[3],t],vars[3]]

Example 12

B[12]={{-8,0,-7,10,-9,1},{9,2,3,-10,-1,-4},{9,0,-2,-9,-3,7}};
dbgGroebnerBasisUnique[HemmeckeToricGroebnerBasis[B[12],vars[6],lexOrder[6]],vars[6]]

dbgGroebnerBasisUnique[GroebnerBasis`ToricGroebnerBasis[Append[VectorToBinomial[LatticeByNullSpace[B[12]],vars[6]],1-t*Apply[Times,vars[6]]],vars[6],t],vars[6]]

Example 13

B[13]={{1,-7,-2,10,3},{-9,7,-4,10,6}};
dbgGroebnerBasisUnique[HemmeckeToricGroebnerBasis[B[13],vars[5],lexOrder[5]],vars[5]]

dbgGroebnerBasisUnique[GroebnerBasis`ToricGroebnerBasis[Append[VectorToBinomial[LatticeByNullSpace[B[13]],vars[5]],1-t*Apply[Times,vars[5]]],vars[5],t],vars[5]]

Example 14

B[14]={{-5,4,-1,-3,-3,10},{0,3,1,-7,7,-1},{7,-2,10,8,8,4}};
profile[dbgGroebnerBasisUnique[HemmeckeToricGroebnerBasis[B[14],vars[6],lexOrder[6]],vars[6]]]

KDL semibin, bitmasks, geometricBuchberger, profiling, LatticeReduce, intermediateReduce, FindInstance

VEC, bitmasks, geometricBuchberger, profiling, LatticeReduce, intermediateReduce, FindInstance

VEC, bitmasks, geometricBuchberger, profiling, LatticeReduce, part intermediateReduce, FindInstance

VEC, bitmasks, geometricBuchberger, profiling, LatticeReduce, part intermediateReduce, Minimize

kdl semibin, bitmasks, geometricBuchberger, profiling, LatticeReduce, part intermediateReduce, Minimize

VEC, bitmasks, geometricBuchberger, profiling, LatticeReduce, part intermediateReduce

kdl semibin, LatticeReduce+bounded variables

vec

kdl fullbin + vec

kdl nobin + vec (=kdl semibin + vec)

kdl fullbin

kdl semibin

kdl nobin

{31.8-10.3,132.0-106.3,51.0-22.6,124.4-105.5,43.8-26.3,43.2-26.3}

Example 15

B[15]={{1,-9,2,-6,-7,2},{-8,-1,7,7,-9,8},{10,-3,9,5,5,10}};
profile[dbgGroebnerBasisUnique[HemmeckeToricGroebnerBasis[B[15],vars[6],lexOrder[6]],vars[6]]]

KDL semibin, bitmasks, geometricBuchberger, profiling, LatticeReduce, intermediateReduce, FindInstance

VEC, bitmasks, geometricBuchberger, profiling, LatticeReduce, intermediateReduce, FindInstance

VEC, bitmasks, geometricBuchberger, profiling, LatticeReduce, part intermediateReduce, FindInstance

VEC, bitmasks, geometricBuchberger, profiling, LatticeReduce, part intermediateReduce, Minimize

kdl semibin, bitmasks, geometricBuchberger, profiling, LatticeReduce, part intermediateReduce, Minimize

VEC, bitmasks, geometricBuchberger, profiling, LatticeReduce, part intermediateReduce

kdl semibin, LatticeReduce+bounded variables

vec

kdl fullbin + vec

kdl nobin + vec (=kdl semibin + vec)

kdl fullbin

kdl semibin

kdl nobin

{55.3-7.1,64.3-2.2,61.8-4.0,59.9-2.2,58.0-4.0,53.1-4.0}

Example 16

B[16]={{-5,6,8,3,5,3},{7,9,-9,8,-3,10},{-2,4,-9,-3,-8,-5}};
profile[dbgGroebnerBasisUnique[HemmeckeToricGroebnerBasis[B[16],vars[6],lexOrder[6]],vars[6]]]

KDL semibin, bitmasks, geometricBuchberger, profiling, LatticeReduce, intermediateReduce, FindInstance

VEC, bitmasks, geometricBuchberger, profiling, LatticeReduce, intermediateReduce, FindInstance

VEC, bitmasks, geometricBuchberger, profiling, LatticeReduce, part intermediateReduce, FindInstance

VEC, bitmasks, geometricBuchberger, profiling, LatticeReduce, part intermediateReduce, Minimize

kdl semibin, bitmasks, geometricBuchberger, profiling, LatticeReduce, part intermediateReduce, Minimize

VEC, bitmasks, geometricBuchberger, profiling, LatticeReduce, part intermediateReduce

kdl semibin, LatticeReduce+bounded variables

vec

kdl fullbin + vec

kdl nobin + vec (=kdl semibin + vec)

kdl fullbin

kdl semibin

kdl nobin

{56.4-10.3,130.8-76.6,58.6-8.5,108.4-75.9,36.9-8.4,35.3-8.4}

Example 17

B[17]={{-2,-3,-7,6,-5,0},{-3,2,-4,9,-1,9}};
profile[dbgGroebnerBasisUnique[HemmeckeToricGroebnerBasis[B[17],vars[6],lexOrder[6]],vars[6]]]

KDL semibin, bitmasks, geometricBuchberger, profiling, LatticeReduce, intermediateReduce, FindInstance

VEC, bitmasks, geometricBuchberger, profiling, LatticeReduce, intermediateReduce, FindInstance

VEC, bitmasks, geometricBuchberger, profiling, LatticeReduce, part intermediateReduce, FindInstance

VEC, bitmasks, geometricBuchberger, profiling, LatticeReduce, part intermediateReduce, Minimize

kdl semibin, bitmasks, geometricBuchberger, profiling, LatticeReduce, part intermediateReduce, Minimize

VEC, bitmasks, geometricBuchberger, profiling, LatticeReduce, part intermediateReduce

kdl semibin, LatticeReduce+bounded variables

kdl fullbin + vec

kdl nobin + vec (=kdl semibin + vec)

kdl fullbin

kdl semibin

kdl nobin

{547.7-525.2,97.1-79.2,80.1-72.7,85.9-78.4,86.1-78.4}

Example 18

B[18]={{912489947,1256837063,1523174407,4535072249,6927997597,8532266279}};
order[18]={{912489947,1256837063,1523174407,4535072249,6927997597,8532266279},{-1,0,0,0,0,0},{0,-1,0,0,0,0},{0,0,-1,0,0,0},{0,0,0,-1,0,0},{0,0,0,0,-1,0}};
profile[HemmeckeToricGroebnerBasis[B[18],vars[6],order[18]]]



KDL semibin, bitmasks, geometricBuchberger, profiling, LatticeReduce, intermediateReduce, FindInstance

VEC, bitmasks, geometricBuchberger, profiling, LatticeReduce, intermediateReduce, FindInstance

VEC, bitmasks, geometricBuchberger, profiling, LatticeReduce, part intermediateReduce, FindInstance

VEC, bitmasks, geometricBuchberger, profiling, LatticeReduce, part intermediateReduce, Minimize

kdl semibin, bitmasks, geometricBuchberger, profiling, LatticeReduce, part intermediateReduce, Minimize

profile[HemmeckeToricGroebnerBasis[{{-5,4,-1,-3,-3,10},{0,3,1,-7,7,-1},{7,-2,10,8,8,4}},vars[6],lexOrder[6]]]

HemmeckeToricGroebnerBasis[{{-1,3,-5}},{x,y,z},lexOrder[3]]

Verification Functions

The macros LM and S compute the leading monomial and the S-polynomial of a polynomial.

SyntaxInformation[LM]={"ArgumentsPattern"->{_,_, _}};
LM[f_,order_,vars_]:=First[MonomialList[f,vars,order]]

SyntaxInformation[S]={"ArgumentsPattern"->{_,_, _,_}};
S[f_,g_,order_,vars_]:=Expand[(f*LM[g,order,vars]-g*LM[f,order,vars])/PolynomialGCD[LM[f,order,vars],LM[g,order,vars]]]

dbgGroebnerBasisUnique brings a reduced Gröbner basis into a unique representation by sorting the polynomials and assigning unique signs. The leading coefficients are assumed to be 1.

SyntaxInformation[dbgGroebnerBasisUnique]={"ArgumentsPattern"->{_,_}};
dbgGroebnerBasisUnique[G_List,variables_List]:=Sort[Map[(If[(LM[#,Lexicographic,variables]/.Thread[variables->1])>0,#,-#])&,G]]

dbgGroebnerBasisQ checks whether the G (given in vector form) is a Gröbner basis with respect to monomialOrder and projection.

SyntaxInformation[dbgGroebnerBasisQ]={"ArgumentsPattern"->{_,_, _}};dbgGroebnerBasisQ[G_List,monomialOrder_List,projection_List]:=Module[{GP,mop,n=Length[projection[[ToricGB`Private`vectorList]]],GPP,i,j},
GP=Map[(#[[projection[[ToricGB`Private`vectorList]]]])&,binomialStrip[G]];
mop=Map[(#[[projection[[ToricGB`Private`vectorList]]]])&,monomialOrder];
mop=mop[[latticeZeroKernelProjection[Transpose[mop]][[ToricGB`Private`vectorList]]]];
GPP=VectorToBinomial[GP,vars[n]];
For[i=1,i≤Length[GPP],i++,
For[j=i+1,j≤Length[GPP],j++,
If[PolynomialReduce[S[GPP[[i]],GPP[[j]],mop,vars[n]],GPP,vars[n],MonomialOrder->mop][[2]]=!=0,
Print["groebnerBasisQ: G=",G,", monomialOrder=",monomialOrder,", projection=",projection," fails at ",{i,j}];
Return[False]];
];
];
True
]

dbgSaturatedQ checks whether the ideal F (given in vector form) is saturated with respect to all variables.

SyntaxInformation[dbgSaturatedQ]={"ArgumentsPattern"->{_}};dbgSaturatedQ[F_List]:=Module[{n=Length[F[[1]][[ToricGB`Private`binomialExponent]]],FP},
FP=VectorToBinomial[binomialStrip[F],vars[n]];
If[dbgGroebnerBasisUnique[GroebnerBasis`ToricGroebnerBasis[Append[FP,1-t*Apply[Times,vars[n]]],vars[n],t],vars[n]]===dbgGroebnerBasisUnique[GroebnerBasis[FP,vars[n]],vars[n]],True,Print[dbgGroebnerBasisUnique[GroebnerBasis`ToricGroebnerBasis[Append[FP,1-t*Apply[Times,vars[n]]],vars[n],t],vars[n]]];Print[dbgGroebnerBasisUnique[GroebnerBasis[FP,vars[n]],vars[n]]]];
]

Profiling

resetProfiling:=(timeGroebner=0;timeRedGB=0;timePrereduce=0;timeCM=0;timeInterRed=0;timeMinimize=0;timeRebuild=0;
cntCS=0;cntCT=0;cntCM=0;
cntReduce0=0;cntReduce1=0;cntReduce2=0;cntReduce2Q=0;
cntReduceQDivQ=0;cntReduceDivQ=0;cntReduceDiv=0;cntMSDivQ=0;cntDividesIter=0;
cntInsert=0;cntChange=0;cntDelete=0;cntQuery=0;cntRebuild=0;cntKDCallback=0;cntCompare=0;cntKDRecursions=0;cntBinaryIter=0;sizeG=0;
timeDiv=0;cntDiv=0;maxExponent={})
resetProfiling;

printProfiling:=Print["timePrereduce=",timePrereduce,", timeGroebner=",timeGroebner,", timeRedGB=",timeRedGB,", timeCM=",timeCM,", timeInterRed=",timeInterRed,", timeMinimize=",timeMinimize,", timeRebuild=", timeRebuild,", cntCS=",cntCS,", cntCT=",cntCT,", cntCM=",cntCM,", cntReduce0=",cntReduce0,", cntReduce1=",cntReduce1,", cntReduce2=",cntReduce2,", cntReduce2Q=",cntReduce2Q,", cntReduceQDivQ=",cntReduceQDivQ,", cntReduceDivQ=",cntReduceDivQ,", cntReduceDiv=",cntReduceDiv,", cntMSDivQ=",cntMSDivQ,", cntDividesIter=",cntDividesIter,", cntInsert=",cntInsert,", cntChange=",cntChange,", cntDelete=",cntDelete,", cntQuery=",cntQuery,", cntRebuild=",cntRebuild,", cntKDCallback=",cntKDCallback,", cntKDRecursions=",cntKDRecursions,", cntCompare=",cntCompare,", cntBinaryIter=",cntBinaryIter,", maxExponent=",maxExponent];

SetAttributes[profile,HoldAll];
profile[func_]:=(resetProfiling;Print[AbsoluteTiming[func][[1]]];printProfiling;)

SetAttributes[addTime,HoldAll];
addTime[func_,var_]:=Module[{return,time},
{time,return}=AbsoluteTiming[func];
var+=time;
return
];

Common Parameters for Examples

lexOrder[n_]:=Array[(UnitVector[n,#])&,n];
var[n_]:=Symbol[FromCharacterCode[n-1+ToCharacterCode["a"]]];
vars[n_]:=Array[(var[#])&,n];
proj[k_,n_]:=positionsToBitMask[Table[i,{i,k}],n];

Data Structures

Queue (from http://library.wolfram.com/conferences/devconf99/lichtblau)

ToricGB`Private`queuePointer=1;
ToricGB`Private`queueElement=2;
ToricGB`Private`queueFront=3;
ToricGB`Private`queueBack=4;
ClearAll[newQueue,clearQueue,queueLength,queueEmpty,enQueue,deQueue];

SyntaxInformation[newQueue]={"ArgumentsPattern"->{}};
newQueue[]:={Unique[],Null,0,0};

SyntaxInformation[clearQueue]={"ArgumentsPattern"->{_}};
SetAttributes[clearQueue,HoldFirst];
clearQueue[Q_Symbol]:=(ClearAll[Evaluate[Q[[ToricGB`Private`queuePointer]]]];ClearAll[Q])

SyntaxInformation[queueLength]={"ArgumentsPattern"->{_}};
SetAttributes[queueLength,HoldFirst];
queueLength[Q_Symbol]:=(Q[[ToricGB`Private`queueBack]]-Q[[ToricGB`Private`queueFront]])

SyntaxInformation[queueEmpty]={"ArgumentsPattern"->{_}};
SetAttributes[queueEmpty,HoldFirst];
queueEmpty[Q_Symbol]:=(Q[[ToricGB`Private`queueBack]]===Q[[ToricGB`Private`queueFront]])

SyntaxInformation[enQueue]={"ArgumentsPattern"->{_,_}};
SetAttributes[enQueue,HoldFirst];
enQueue[Q_Symbol,elem_]:=(Q[[ToricGB`Private`queueFront]]++;Q[[ToricGB`Private`queuePointer]][Q[[ToricGB`Private`queueFront]]]=elem;)

SyntaxInformation[deQueue]={"ArgumentsPattern"->{_}};
SetAttributes[deQueue,HoldFirst];
deQueue[Q_Symbol]:=(Q[[ToricGB`Private`queueBack]]++;Q[[ToricGB`Private`queueElement]]=Q[[ToricGB`Private`queuePointer]][Q[[ToricGB`Private`queueBack]]];Q[[ToricGB`Private`queuePointer]][Q[[ToricGB`Private`queueBack]]]=.;Q[[ToricGB`Private`queueElement]])

Binary Search

Search the position in which x should be inserted into L such that it remains ordered. If x is contained in L, the position of x in L is returned.

SyntaxInformation[binarySearch]={"ArgumentsPattern"->{_,_,_}};
binarySearch[L_List,x_,ordering_]:=Module[{a,b,m,r},
a=1;b=Length[L]+1;
While[a<b,
cntBinaryIter++;
m=Quotient[a+b,2];
r=ordering[x,L[[m]]];
If[r>0,a=m+1,b=m];
];
Return[a];
];

binarySearch[{1,4,6,9},5,(#1-#2)&]

binarySearch[{1,4,6,9},9,(#1-#2)&]

binarySearch[{1,4,6,9},12,(#1-#2)&]

kd-List

A kd-list is a k-dimensional (linked) list, similar to kd-Trees, which is used to store k-dimensional integer vectors. The first level is a list containing the first component of each vector (each value is only stored once). Then each entry x points to a list of level two in which we store the second component of the vectors whose first component has value x. This is continued up to level k.

This implementation maintains the lists at each level in a sorted order, resulting in slow insertion times but very fast queries of the kind "all entries with components smaller x". Behaviour of this kind is desirable since the list is rarely updated and many queries are run. It operates on binomials with extra information, v[[ToricGB`Private`binomialExponent]] being the integer vector used as "key" and v[[ToricGB`Private`binomialIndex]] being the index of the vector in the plain list that is maintained for direct access. Additionally, only the components listed in projection are used for branching. Those vectors whose coordinate is negative or zero are stored in one bucket with value 0 (this is possible since all queries are nonnegative).

It is possiple to insert, change and (kind of) delete entries. The functions with the parameter hint are helper functions that only perform a partial job and not meant to be called externally.

The function queryKDList iterates over all entries of the kd-list that are componentwise less or equal than the query. Depending on the parameter largeFirst and whether the query value is positive or zero, it starts with the vectors with large respectively small entries. For each entry the callback function is called. As soon as the callback function returns a value different from Null, the query is aborted and the value of the callback function is returned. The vector query must be nonnegative. If inverseQuery is True, instead all elements that are componentwise greater or equal than the query are looked for.

The head symbol of the KDList is stored in the first entry of the list returned by newKDList. The second stores the projection which is passed as parameter. The third is the head symbol of a sparse hash table which is used to disable elements in the list (for the use in binomialReduce). Disabled elements are stored with a value of 0, others are not stored in the hash table. The last entry of a KDList is a plain list of all elements.

Attention

Inserting an element twice breaks the data structure!

ToricGB`Private`basisList=1;
ToricGB`Private`basisProjection=2;
ToricGB`Private`basisKDList=3;
ToricGB`Private`basisChanges=4;
ClearAll[newKDList,clearKDList,insertKDList,deleteKDList,changeKDList,queryKDList];

SyntaxInformation[newKDList]={"ArgumentsPattern"->{_}};
newKDList[projection_List]:=Module[{KDL},
dbgLog["New kd list with projection=",projection];
KDL={{},projection,Unique[]};
KDL[[ToricGB`Private`basisKDList]][{}]={};
KDL
]

SyntaxInformation[clearKDList]={"ArgumentsPattern"->{_}};
SetAttributes[clearKDList,HoldFirst];
clearKDList[KDL_Symbol]:=(ClearAll[Evaluate[KDL[[ToricGB`Private`basisKDList]]]];ClearAll[KD])

SyntaxInformation[rebuildKDList]={"ArgumentsPattern"->{_}};
SetAttributes[rebuildKDList,HoldFirst];
rebuildKDList[KDT_Symbol,idx_Integer]:=idx;

SyntaxInformation[insertKDList]={"ArgumentsPattern"->{_,_,_.}};
SetAttributes[insertKDList,HoldFirst];
insertKDList[KDL_Symbol,v_List,hint_List]:=Module[{it,pos,level=Length[hint]+1,coord,value},
If[level>Length[KDL[[ToricGB`Private`basisProjection,ToricGB`Private`vectorList]]],
(* end of projection -> simple list *)
KDL[[ToricGB`Private`basisKDList]][hint]=Append[KDL[[ToricGB`Private`basisKDList]][hint],v];
,
(* search v[[ToricGB`Private`binomialExponent,coord]] in level 'hint' *)
coord=KDL[[ToricGB`Private`basisProjection,ToricGB`Private`vectorList,level]];
value=Max[v[[ToricGB`Private`binomialExponent,coord]],0];
it=Append[hint,value];
pos=binarySearch[KDL[[ToricGB`Private`basisKDList]][hint],value,(#1-#2)&];
If[pos>Length[KDL[[ToricGB`Private`basisKDList]][hint]]||KDL[[ToricGB`Private`basisKDList]][hint][[pos]]=!=value,
(* not found -> insert *)
KDL[[ToricGB`Private`basisKDList]][hint]=Insert[KDL[[ToricGB`Private`basisKDList]][hint],value,pos];
KDL[[ToricGB`Private`basisKDList]][it]={};
];
(* go to next level *)
insertKDList[KDL,v,it];
];
]

insertKDList[KDL_Symbol,v_List]:=Module[{it,pos,coord,u=v,value},
cntInsert++;
dbgLog["insertKDList[",v,"]"];
u[[ToricGB`Private`binomialIndex]]=Length[KDL[[ToricGB`Private`basisList]]]+1;
u[[ToricGB`Private`binomialMarker]]=1;
(* append to list of all elements *)
KDL[[ToricGB`Private`basisList]]=Append[KDL[[ToricGB`Private`basisList]],u];
(* insert in real kd-list *)
insertKDList[KDL,u,{}];
]

SyntaxInformation[deleteKDList]={"ArgumentsPattern"->{_,_,_}};
SetAttributes[deleteKDList,HoldFirst];
(* returns True if list KDL[[ToricGB`Private`basisKDList]][hint] === {} after call *)
deleteKDList[KDL_Symbol,v_List,hint_]:=Module[{it,pos,level=Length[hint]+1,coord,value},
If[level>Length[KDL[[ToricGB`Private`basisProjection,ToricGB`Private`vectorList]]],
(* end of projection -> simple list *)
KDL[[ToricGB`Private`basisKDList]][hint]=Select[KDL[[ToricGB`Private`basisKDList]][hint],(#[[ToricGB`Private`binomialIndex]]=!=v[[ToricGB`Private`binomialIndex]])&];
,
(* search v[[ToricGB`Private`binomialExponent,coord]] in level 'hint' *)
coord=KDL[[ToricGB`Private`basisProjection,ToricGB`Private`vectorList,level]];
value=Max[v[[ToricGB`Private`binomialExponent,coord]],0];
it=Append[hint,value];
pos=binarySearch[KDL[[ToricGB`Private`basisKDList]][hint],value,(#1-#2)&];
assert[pos<=Length[KDL[[ToricGB`Private`basisKDList]][hint]]&&KDL[[ToricGB`Private`basisKDList]][hint][[pos]]===value,deleteKDList,"Could not find ",v," in kd-list."];
(* go to next level *)
If[deleteKDList[KDL,v,it],
KDL[[ToricGB`Private`basisKDList]][it]=.;
KDL[[ToricGB`Private`basisKDList]][hint]=Complement[KDL[[ToricGB`Private`basisKDList]][hint],{value}];
];
];
KDL[[ToricGB`Private`basisKDList]][hint]==={}
];

deleteKDList[KD_Symbol,index_Integer]:=Module[{u=v},
cntDelete++;
dbgLog["deleteKDList[",index,"]"];
KD[[ToricGB`Private`basisList,index,ToricGB`Private`binomialMarker]]=0;
deleteKDList[KD,KD[[ToricGB`Private`basisList,index]],{}];
];

SyntaxInformation[changeKDList]={"ArgumentsPattern"->{_,_,_}};
SetAttributes[changeKDList,HoldFirst];
changeKDList[KD_Symbol,index_Integer,v_List]:=Module[{it,pos,coord,u=v,value},
cntChange++;
dbgLog["changeKDList[",index,", ",v,"]"];
deleteKDList[KD,KD[[ToricGB`Private`basisList,index]],{}];
u[[ToricGB`Private`binomialIndex]]=index;
u[[ToricGB`Private`binomialMarker]]=1;
KD[[ToricGB`Private`basisList,index]]=u;
insertKDList[KD,u,{}];
]

SyntaxInformation[queryKDList]={"ArgumentsPattern"->{_,_,_,_,_.,_.,_.}};
SetAttributes[queryKDList,HoldAll];
queryKDList[KDL_Symbol,query_,cnt_,callback_,inverseQuery_,largeFirst_,hint_]:=Module[{it,level=Length[hint]+1,coord,pos,return,value},
cntKDRecursions++;
If[level>Length[KDL[[ToricGB`Private`basisProjection,ToricGB`Private`vectorList]]],
(* end of projection -> simple list *)
For[pos=1,pos≤Length[KDL[[ToricGB`Private`basisKDList]][hint]],pos++,
cntKDCallback++;
return=callback[KDL,query,KDL[[ToricGB`Private`basisKDList]][hint][[pos]]];
If[return=!=Null,Sow[return];If[--cnt==0,Return[]]];
],
(* traverse list in level hint up to value v[[coord]] *)
coord=KDL[[ToricGB`Private`basisProjection,ToricGB`Private`vectorList,level]];
value=query[[ToricGB`Private`binomialExponent,coord]];
it=Append[hint,0];
If[inverseQuery,
(* look for greater values than query *)
For[pos=Length[KDL[[ToricGB`Private`basisKDList]][hint]],pos≥1&&KDL[[ToricGB`Private`basisKDList]][hint][[pos]]>=value,pos--,
it[[level]]=KDL[[ToricGB`Private`basisKDList]][hint][[pos]];
queryKDList[KDL,query,cnt,callback,inverseQuery,largeFirst,it];
If[cnt==0,Return[]];
];
,
(* look for smaller values than query *)
If[!largeFirst||value≤0,
(* start search with small values *)
For[pos=1,pos<=Length[KDL[[ToricGB`Private`basisKDList]][hint]]&&KDL[[ToricGB`Private`basisKDList]][hint][[pos]]≤value,pos++,
it[[level]]=KDL[[ToricGB`Private`basisKDList]][hint][[pos]];
queryKDList[KDL,query,cnt,callback,inverseQuery,largeFirst,it];
If[cnt==0,Return[]];
];
,
(* start search with large values -> binary search *)
pos=binarySearch[KDL[[ToricGB`Private`basisKDList]][hint],value,(#1-#2)&];
If[pos>Length[KDL[[ToricGB`Private`basisKDList]][hint]]||KDL[[ToricGB`Private`basisKDList]][hint][[pos]]>value,pos--];
While[pos>0,
it[[level]]=KDL[[ToricGB`Private`basisKDList]][hint][[pos]];
queryKDList[KDL,query,cnt,callback,inverseQuery,largeFirst,it];
If[cnt==0,Return[]];
pos--;
];
];
];
];
]

queryKDList[KD_Symbol,query_,count_,callback_,inverseQuery_:False,largeFirst_:False]:=Module[{cnt=count,results},
cntQuery++;
results=Reap[queryKDList[KD,query,cnt,callback,inverseQuery,largeFirst,{}]][[2]];dbgLog["queryKDList[",query,", ",callback,", ",count,", ",inverseQuery,", ",largeFirst,"] = ",results];
results
]

kdl=newKDList[proj[3,3]]
insertKDList[kdl,binomialSupport[{4,-4,10}]];
insertKDList[kdl,binomialSupport[{4,-1,2}]];
insertKDList[kdl,binomialSupport[{4,3,1}]];
insertKDList[kdl,binomialSupport[{4,5,-1}]];

Information[Evaluate[kdl[[ToricGB`Private`basisKDList]]]]

queryKDList[kdl,binomialSupport[{4,4,9}],Infinity,(#3)&]

changeKDList[kdl,2,binomialSupport[{5,4,2}]];

Information[Evaluate[kdl[[ToricGB`Private`basisKDList]]]]

queryKDList[kdl,binomialSupport[{4,0,4}],1,(#3)&,False]

queryKDList[kdl,binomialSupport[{4,0,4}],1,(#3)&,True]

kd-Tree

ClearAll[newKDTree,clearKDTree,rebuildKDTree,insertKDTree,deleteKDTree,changeKDTree,queryKDTree];
ToricGB`Private`kdtRangeStart=1;
ToricGB`Private`kdtRangeEnd=2;
ToricGB`Private`kdtIndex=1;
ToricGB`Private`kdtSeparator=2;
ToricGB`Private`kdtLowMin=3;
ToricGB`Private`kdtLowMax=4;
ToricGB`Private`kdtHighMin=5;
ToricGB`Private`kdtHighMax=6;
ToricGB`Private`kdtCutOff=1;

SyntaxInformation[newKDTree]={"ArgumentsPattern"->{_}};
newKDTree[projection_List]:=Module[{KDT},
dbgLog["New kd tree with projection=",projection];
KDT={{},projection,Unique[],0};
KDT[[ToricGB`Private`basisKDList]][Table[{0,Infinity},{Length[projection[[ToricGB`Private`vectorList]]]}]]={};
KDT
]

SyntaxInformation[clearKDTree]={"ArgumentsPattern"->{_}};
SetAttributes[clearKDTree,HoldFirst];
clearKDTree[KDT_Symbol]:=(ClearAll[Evaluate[KDT[[ToricGB`Private`basisKDList]]]];ClearAll[KDT])

SyntaxInformation[rebuildKDTree]={"ArgumentsPattern"->{_,_.,_.,_.}};
SetAttributes[rebuildKDTree,HoldFirst];
rebuildKDTree[KDT_Symbol,F_List,hint_List,levelP_Integer]:=Module[{coord,G,mid,value,it,level=levelP},
(*Print["rebuildKDTree[",F,",",hint,",",level,"]"];*)
If[Length[F]≤ToricGB`Private`kdtCutOff,
(* rest of elements as list *)
KDT[[ToricGB`Private`basisKDList]][hint]=F;
,
(* build tree recursively *)
While[True,
(* search feasible coordinate for splitting *)
coord=KDT[[ToricGB`Private`basisProjection,ToricGB`Private`vectorList,level]];
G=Sort[F,(#1[[ToricGB`Private`binomialExponent,coord]]<#2[[ToricGB`Private`binomialExponent,coord]])&];
If[Max[G[[1,ToricGB`Private`binomialExponent,level]],0]≠Max[G[[Length[G],ToricGB`Private`binomialExponent,level]],0],Break[]];
level=Mod[level,Length[KDT[[ToricGB`Private`basisProjection,ToricGB`Private`vectorList]]]]+1;
If[level==levelP,
(* no coordinate found *)
KDT[[ToricGB`Private`basisKDList]][hint]=F;
Return[];
];
];
mid=Quotient[Length[G],2]+1;
value=Max[G[[mid,ToricGB`Private`binomialExponent,coord]],0];
While[mid>1&&value==Max[G[[mid-1,ToricGB`Private`binomialExponent,coord]],0],mid--];
If[mid==1&&value==Max[G[[1,ToricGB`Private`binomialExponent,coord]],0],
While[value==Max[G[[mid,ToricGB`Private`binomialExponent,coord]],0],mid++];
value=Max[G[[mid,ToricGB`Private`binomialExponent,coord]],0];
];
KDT[[ToricGB`Private`basisKDList]][hint]={level,value,Max[G[[1,ToricGB`Private`binomialExponent,coord]],0],Max[G[[mid-1,ToricGB`Private`binomialExponent,coord]],0],value,Max[G[[Length[G],ToricGB`Private`binomialExponent,coord]],0]};
(* lower values *)
it=hint;it[[level,ToricGB`Private`kdtRangeEnd]]=value;
(*Print[it];*)
rebuildKDTree[KDT,G[[;;mid-1]],it,Mod[level,Length[KDT[[ToricGB`Private`basisProjection,ToricGB`Private`vectorList]]]]+1];
(* higher values *)
it=hint;it[[level,ToricGB`Private`kdtRangeStart]]=value;
(*Print[it];*)
rebuildKDTree[KDT,G[[mid;;]],it,Mod[level,Length[KDT[[ToricGB`Private`basisProjection,ToricGB`Private`vectorList]]]]+1];
];
];

rebuildKDTree[KDT_Symbol,idx_Integer]:=Module[{u,newidx,i},
If[KDT[[ToricGB`Private`basisChanges]]<0.3*Length[KDT[[ToricGB`Private`basisList]]],Return[idx]];
KDT[[ToricGB`Private`basisChanges]]=0;
cntRebuild++;
dbgLog["rebuildKDTree[",idx,"]"];
(*Print["rebuildKDTree[",idx,"]"];*)
u=KDT[[ToricGB`Private`basisList,idx]];
dbgAssert[u[[ToricGB`Private`binomialMarker]]===1,rebuildKDTree,"Binomial with index ",idx," is not active."];
KDT[[ToricGB`Private`basisList]]=Select[KDT[[ToricGB`Private`basisList]],(#[[ToricGB`Private`binomialMarker]]===1)&];
newidx=Position[KDT[[ToricGB`Private`basisList]],u,{1},1,Heads->False][[1,1]];
For[i=1,i≤Length[KDT[[ToricGB`Private`basisList]]],i++,KDT[[ToricGB`Private`basisList,i,ToricGB`Private`binomialIndex]]=i];
ClearAll[Evaluate[KDT[[ToricGB`Private`basisKDList]]]];
rebuildKDTree[KDT,KDT[[ToricGB`Private`basisList]],Table[{0,Infinity},{Length[KDT[[ToricGB`Private`basisProjection,ToricGB`Private`vectorList]]]}],1];
(*Information[Evaluate[KDT[[ToricGB`Private`basisKDList]]]];*)
newidx
];

SyntaxInformation[insertKDTree]={"ArgumentsPattern"->{_,_,_.}};
SetAttributes[insertKDTree,HoldFirst];
insertKDTree[KDT_Symbol,v_List,hint_List]:=Module[{it,coord,value,level},
(*Print["insertKDTree[",v,",",hint];
Print[KDT[[ToricGB`Private`basisKDList]][hint]];*)
If[KDT[[ToricGB`Private`basisKDList]][hint]==={}||Head[KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`binomialExponent]]]===List,
(* rest of elements as list *)
KDT[[ToricGB`Private`basisKDList]][hint]=Append[KDT[[ToricGB`Private`basisKDList]][hint],v];
,
level=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtIndex]];
coord=KDT[[ToricGB`Private`basisProjection,ToricGB`Private`vectorList,level]];
value=Max[v[[ToricGB`Private`binomialExponent,coord]],0];
it=hint;
If[value>=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtSeparator]],
KDT[[ToricGB`Private`basisKDList]][hint]={KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtIndex]],KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtSeparator]],KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtLowMin]],KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtLowMax]],Min[value,KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtHighMin]]],
Max[value,KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtHighMax]]]};
it[[level,ToricGB`Private`kdtRangeStart]]=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtSeparator]];
,
KDT[[ToricGB`Private`basisKDList]][hint]={KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtIndex]],KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtSeparator]],Min[value,KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtLowMin]]],Max[value,KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtLowMax]]],KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtHighMin]],
KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtHighMax]]};
it[[level,ToricGB`Private`kdtRangeEnd]]=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtSeparator]];
];
insertKDTree[KDT,v,it];
];
]

insertKDTree[KDT_Symbol,v_List]:=Module[{it,pos,coord,u=v,value},
cntInsert++;
KDT[[ToricGB`Private`basisChanges]]++;
dbgLog["insertKDTree[",v,"]"];
u[[ToricGB`Private`binomialIndex]]=Length[KDT[[ToricGB`Private`basisList]]]+1;
u[[ToricGB`Private`binomialMarker]]=1;
(* append to list of all elements *)
KDT[[ToricGB`Private`basisList]]=Append[KDT[[ToricGB`Private`basisList]],u];
(* insert in real kd-list *)
insertKDTree[KDT,u,Table[{0,Infinity},{Length[KDT[[ToricGB`Private`basisProjection,ToricGB`Private`vectorList]]]}]];
(*Information[Evaluate[KDT[[ToricGB`Private`basisKDList]]]];*)
]

SyntaxInformation[deleteKDTree]={"ArgumentsPattern"->{_,_,_.}};
SetAttributes[deleteKDTree,HoldFirst];
deleteKDTree[KDT_Symbol,v_List,hint_]:=Module[{it,coord,value,level},
(*Print["deleteKDTree[",v,",",hint,"]"];*)
dbgAssert[KDT[[ToricGB`Private`basisKDList]][hint]=!={},deleteKDTree,"Could not find ",v," in kd-tree."];
If[Head[KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`binomialExponent]]]===List,
(* rest of elements as list *)
KDT[[ToricGB`Private`basisKDList]][hint]=Select[KDT[[ToricGB`Private`basisKDList]][hint],(#[[ToricGB`Private`binomialIndex]]=!=v[[ToricGB`Private`binomialIndex]])&];
,
level=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtIndex]];
coord=KDT[[ToricGB`Private`basisProjection,ToricGB`Private`vectorList,level]];
value=Max[v[[ToricGB`Private`binomialExponent,coord]],0];
it=hint;
If[value>=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtSeparator]],
it[[level,ToricGB`Private`kdtRangeStart]]=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtSeparator]];
,
it[[level,ToricGB`Private`kdtRangeEnd]]=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtSeparator]];
];
deleteKDTree[KDT,v,it];
];
];

deleteKDTree[KDT_Symbol,index_Integer]:=Module[{},
cntDelete++;
KDT[[ToricGB`Private`basisChanges]]++;
dbgLog["deleteKDTree[",index,"]"];
KDT[[ToricGB`Private`basisList,index,ToricGB`Private`binomialMarker]]=0;
deleteKDTree[KDT,KDT[[ToricGB`Private`basisList,index]],Table[{0,Infinity},{Length[KDT[[ToricGB`Private`basisProjection,ToricGB`Private`vectorList]]]}]];
];

SyntaxInformation[changeKDTree]={"ArgumentsPattern"->{_,_,_}};
SetAttributes[changeKDTree,HoldFirst];
changeKDTree[KDT_Symbol,index_Integer,v_List]:=Module[{it,pos,coord,u=v,value},
cntChange++;
KDT[[ToricGB`Private`basisChanges]]++;
dbgLog["changeKDTree[",index,", ",v,"]"];
deleteKDTree[KDT,KDT[[ToricGB`Private`basisList,index]],Table[{0,Infinity},{Length[KDT[[ToricGB`Private`basisProjection,ToricGB`Private`vectorList]]]}]];
u[[ToricGB`Private`binomialIndex]]=index;
u[[ToricGB`Private`binomialMarker]]=1;
KDT[[ToricGB`Private`basisList,index]]=u;
insertKDTree[KDT,u,Table[{0,Infinity},{Length[KDT[[ToricGB`Private`basisProjection,ToricGB`Private`vectorList]]]}]];
(*Information[Evaluate[KDT[[ToricGB`Private`basisKDList]]]];*)
]

SyntaxInformation[queryKDTree]={"ArgumentsPattern"->{_,_,_,_,_.,_.,_.,_.}};
SetAttributes[queryKDTree,HoldAll];
queryKDTree[KDT_Symbol,query_,cnt_,callback_,inverseQuery_,largeFirst_,hint_,check_]:=Module[{it,coord,pos,return,value,level},
(*Print["queryKDTree[",query,",",hint,"]"];
Print[KDT[[ToricGB`Private`basisKDList]][hint]];*)
cntKDRecursions++;
If[KDT[[ToricGB`Private`basisKDList]][hint]==={},Return[]];
If[Head[KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`binomialExponent]]]===List,
(* rest of elements as list *)
For[pos=1,pos≤Length[KDT[[ToricGB`Private`basisKDList]][hint]],pos++,
If[inverseQuery,
return=binomialDividesQ[query,KDT[[ToricGB`Private`basisKDList]][hint][[pos]],{check}],
return=binomialDividesQ[KDT[[ToricGB`Private`basisKDList]][hint][[pos]],query,{check}]
];
If[return,
cntKDCallback++;
return=callback[KDT,query,KDT[[ToricGB`Private`basisKDList]][hint][[pos]]];
If[return=!=Null,Sow[return];If[--cnt==0,Return[]]];
];
];
,
level=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtIndex]];
coord=KDT[[ToricGB`Private`basisProjection,ToricGB`Private`vectorList,level]];
value=Max[query[[ToricGB`Private`binomialExponent,coord]],0];

If[inverseQuery,
cntCompare++;
If[value<=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtHighMax]],
it=hint;it[[level,ToricGB`Private`kdtRangeStart]]=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtSeparator]];
cntCompare++;
If[value<=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtHighMin]],
queryKDTree[KDT,query,cnt,callback,inverseQuery,largeFirst,it,Complement[check,{coord}]],
queryKDTree[KDT,query,cnt,callback,inverseQuery,largeFirst,it,check];
];
];
cntCompare++;
If[value<=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtLowMax]],
it=hint;it[[level,ToricGB`Private`kdtRangeEnd]]=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtSeparator]];
cntCompare++;
If[value<=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtLowMin]],
queryKDTree[KDT,query,cnt,callback,inverseQuery,largeFirst,it,Complement[check,{coord}]],
queryKDTree[KDT,query,cnt,callback,inverseQuery,largeFirst,it,check];
];
];
,
If[largeFirst,
cntCompare++;
If[value>=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtHighMin]],
it=hint;it[[level,ToricGB`Private`kdtRangeStart]]=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtSeparator]];
cntCompare++;
If[value>=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtHighMax]],
queryKDTree[KDT,query,cnt,callback,inverseQuery,largeFirst,it,Complement[check,{coord}]],
queryKDTree[KDT,query,cnt,callback,inverseQuery,largeFirst,it,check];
];
];
cntCompare++;
If[value>=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtLowMin]],
it=hint;it[[level,ToricGB`Private`kdtRangeEnd]]=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtSeparator]];
cntCompare++;
If[value>=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtLowMax]],
queryKDTree[KDT,query,cnt,callback,inverseQuery,largeFirst,it,Complement[check,{coord}]],
queryKDTree[KDT,query,cnt,callback,inverseQuery,largeFirst,it,check];
];
];
,
cntCompare++;
If[value>=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtLowMin]],
it=hint;it[[level,ToricGB`Private`kdtRangeEnd]]=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtSeparator]];
cntCompare++;
If[value>=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtLowMax]],
queryKDTree[KDT,query,cnt,callback,inverseQuery,largeFirst,it,Complement[check,{coord}]],
queryKDTree[KDT,query,cnt,callback,inverseQuery,largeFirst,it,check];
];
];
cntCompare++;
If[value>=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtHighMin]],
cntCompare++;
it=hint;it[[level,ToricGB`Private`kdtRangeStart]]=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtSeparator]];
If[value>=KDT[[ToricGB`Private`basisKDList]][hint][[ToricGB`Private`kdtHighMax]],
queryKDTree[KDT,query,cnt,callback,inverseQuery,largeFirst,it,Complement[check,{coord}]],
queryKDTree[KDT,query,cnt,callback,inverseQuery,largeFirst,it,check];
];
];
];
];
];
]

queryKDTree[KDT_Symbol,query_,count_,callback_,inverseQuery_:False,largeFirst_:False]:=Module[{cnt=count,results,k},
cntQuery++;
(*Print[KDT[[ToricGB`Private`basisList]]];
Information[Evaluate[KDT[[ToricGB`Private`basisKDList]]]];
Print["queryKDTree[",query,",",count,",",callback,",",inverseQuery,"]"];*)
k=Length[KDT[[ToricGB`Private`basisProjection,ToricGB`Private`vectorList]]];
results=Reap[queryKDTree[KDT,query,cnt,callback,inverseQuery,largeFirst,Table[{0,Infinity},{k}],KDT[[ToricGB`Private`basisProjection,ToricGB`Private`vectorList]]]][[2]];dbgLog["queryKDTree[",query,", ",callback,", ",count,", ",inverseQuery,", ",largeFirst,"] = ",results];
(*Print["-> ",results];*)
results
]

kdt=newKDTree[proj[3,3]];
kdt[[ToricGB`Private`basisList]]=binomialSupport[{{1,2,3},{1,4,6},{4,5,2},{5,3,5},{4,3,2},{6,6,1}}];
rebuildKDTree[kdt,1];
Information[Evaluate[kdt[[ToricGB`Private`basisKDList]]]]

insertKDTree[kdt,binomialSupport[{4,4,4}]];
deleteKDTree[kdt,2];
Information[Evaluate[kdt[[ToricGB`Private`basisKDList]]]]

queryKDTree[kdt,binomialSupport[{4,4,4}],Infinity,#3&]

queryKDTree[kdt,binomialSupport[{1,4,-2}],Infinity,#3&,True]

Mathematica Lists

These are mapping functions which should provide an interface for Mathematica lists exchangeable with kd-lists.

ClearAll[newMathList,clearMathList,rebuildMathList,deleteMathList,changeMathList,queryMathList];

SyntaxInformation[newMathList]={"ArgumentsPattern"->{_}};
newMathList[projection_List]:={{},projection};

SyntaxInformation[clearMathList]={"ArgumentsPattern"->{_}};
SetAttributes[clearMathList,HoldFirst];
clearMathList[ML_Symbol]:=ClearAll[ML];

SyntaxInformation[rebuildMathList]={"ArgumentsPattern"->{_}};
SetAttributes[rebuildMathList,HoldFirst];
rebuildMathList[ML_Symbol,idx_Integer]:=Module[{u,i,newidx},
cntRebuild++;
dbgLog["rebuildMathList[",idx,"]"];
u=ML[[ToricGB`Private`basisList,idx]];
dbgAssert[u[[ToricGB`Private`binomialMarker]]===1,rebuildMathList,"Binomial with index ",idx," is not active."];
ML[[ToricGB`Private`basisList]]=Select[ML[[ToricGB`Private`basisList]],(#[[ToricGB`Private`binomialMarker]]===1)&];
newidx=Position[ML[[ToricGB`Private`basisList]],u,{1},1,Heads->False][[1,1]];
For[i=1,i≤Length[ML[[ToricGB`Private`basisList]]],i++,ML[[ToricGB`Private`basisList,i,ToricGB`Private`binomialIndex]]=i];
newidx
];

SyntaxInformation[insertMathList]={"ArgumentsPattern"->{_}};
SetAttributes[insertMathList,HoldFirst];
insertMathList[ML_Symbol,v_List]:=Module[{u=v},
cntInsert++;
dbgLog["insertMathList[",v,"]"];
u[[ToricGB`Private`binomialIndex]]=Length[ML[[ToricGB`Private`basisList]]]+1;
u[[ToricGB`Private`binomialMarker]]=1;
ML[[ToricGB`Private`basisList]]=Append[ML[[ToricGB`Private`basisList]],u];
];

SyntaxInformation[changeMathList]={"ArgumentsPattern"->{_}};
SetAttributes[changeMathList,HoldFirst];
changeMathList[ML_Symbol,index_Integer,v_List]:=Module[{u=v},
cntChange++;
dbgLog["changeMathList[",index,", ",v,"]"];
u[[ToricGB`Private`binomialIndex]]=index;
ML[[ToricGB`Private`basisList,index]]=u;
];

SyntaxInformation[deleteMathList]={"ArgumentsPattern"->{_}};
SetAttributes[deleteMathList,HoldFirst];
deleteMathList[ML_Symbol,index_Integer]:=Module[{u=v},
cntDelete++;
dbgLog["deleteMathList[",index,"]"];
ML[[ToricGB`Private`basisList,index,ToricGB`Private`binomialMarker]]=0;
];

SyntaxInformation[queryMathList]={"ArgumentsPattern"->{_,_,_,_,_.,_.,_.}};
SetAttributes[queryMathList,HoldAll];
queryMathList[ML_Symbol,query_,count_,callback_,inverseQuery_:False,largeFirst_:False]:=Module[{i,results,return,cnt=count},
cntQuery++;
results=Reap[For[i=1,i≤Length[ML[[ToricGB`Private`basisList]]],i++,
If[ML[[ToricGB`Private`basisList,i,ToricGB`Private`binomialMarker]]=!=1,Continue[]];
If[inverseQuery,
return=binomialDividesQ[query,ML[[ToricGB`Private`basisList,i]],ML[[ToricGB`Private`basisProjection]]],
return=binomialDividesQ[ML[[ToricGB`Private`basisList,i]],query,ML[[ToricGB`Private`basisProjection]]]
];
If[return,
cntKDCallback++;
return=callback[ML,query,ML[[ToricGB`Private`basisList,i]]];
If[return=!=Null,Sow[return];If[--cnt==0,Break[]]];
];
];
][[2]];
dbgLog["queryMathList[",query,", ",callback,", ",count,", ",inverseQuery,", ",largeFirst,"] = ",results];
results
];

Lattices

LatticeByNullSpace

The integer matrix A defines a lattice L = { v integer vector: A.v=0 }.
LatticeByNullSpace computes a basis B (given as list of basis vectors) of L, i.e. A . B[[i]] = 0 for all i.

SyntaxInformation[LatticeByNullSpace]={"ArgumentsPattern"->{_}};
LatticeByNullSpace[A_List]:=Module[{U,H,N},
{U,H} =HermiteDecomposition[Transpose[A]];
N=NullSpace[Transpose[H],Method->DivisionFreeRowReduction];
If[N==={},{},N=N.U;dbgAssert[A.Transpose[N]===Array[0&,{Length[A],Length[N]}],LatticeByNullSpace,"Lattice basis is not orthogonal to A."];N]
];

A={{2,1,1}};
L=LatticeByNullSpace[A]

Attention

NullSpace[A,Method->DivisionFreeRowReduction]

NullSpace doesn't work here. It give a Q-basis, but not a Z-basis. {-1,1,1} = 0.5 * {-1,0,2} + 0.5 *  {-1,2,0}, but there is no integral representation. Thus {{-1, 0, 2}, {-1, 2, 0}} is a basis of a proper sublattice of  L=span{{-1, 0, 2}, {-1, 1, 1}}.

zeroKernelProjection

Find projection wrt. which lattice L (given as lattice basis) is full-dimensional (MatrixRank[L] = Length[projection]) and which has no kernel.

SyntaxInformation[latticeZeroKernelProjection]={"ArgumentsPattern"->{_}};
latticeZeroKernelProjection[L_List]:=Module[{projection={},i,M=Transpose[L],A={},B},
(* OPT: strategy for choosing projection (= order rows of M differently) *)
positionsToBitMask[Flatten[Map[(Position[#,a_/;a=!=0,{1},1,Heads->False])&,RowReduce[L]]],Length[M]]
]

latticeZeroKernelProjection[{{1,0,1},{3,0,1}}]

latticeZeroKernelProjection[{{1,1,1},{3,0,1}}]

latticePositiveElement

Finds an integral vector in the lattice L = { v integer vector: A.v=0 } whose projection is non-negative and whose i-th coordinate is positive.

SyntaxInformation[latticePositiveElement]={"ArgumentsPattern"->{_,_,_}};latticePositiveElement[L_List,projection_List,i_Integer]:=Module[{x,c,b,lu,var,equ},
(* OPT: Use ray algorithms instead of LinearProgramming -> smaller elements, faster? *)
ClearAll[z];
var=Array[z,Length[L]];
equ=var.L;
x=addTime[FindInstance[Flatten[{equ[[i]]≥1,Thread[equ[[projection[[ToricGB`Private`vectorList]]]]≥0]}],var,Integers],timeMinimize];
If[x==={},{},binomialSupport[equ/.x[[1]]]]
(*Off[Minimize::infeas];
x=addTime[Minimize[{Table[1,{n}].Abs[var],Apply[And,Thread[A.var==0]]&&var[[i]]≥1&&Apply[And,Thread[var[[Complement[projection[[ToricGB`Private`vectorList]],{i}]]]≥0]]},var,Integers],timeMinimize];
x=addTime[Minimize[{Table[1,{Length[projection[[ToricGB`Private`vectorList]]]}].var[[projection[[ToricGB`Private`vectorList]]]],Apply[And,Thread[A.var==0]]&&var[[i]]≥1&&Apply[And,Thread[var[[Complement[projection[[ToricGB`Private`vectorList]],{i}]]]≥0]]},var,Integers],timeMinimize];
On[Minimize::infeas];
If[x[[1]]===Infinity,Null,binomialSupport[Map[(#[[2]])&,x[[2]]]]]*)
(*c=Table[1,{n}];
b=Table[{0,0},{Length[A]}];
lu=Table[{-Infinity,Infinity},{n}];
lu[[projection[[ToricGB`Private`vectorList]],1]]=0;
lu[[i,1]]=1;
Print[c,A,b,lu,Integers];
x=LinearProgramming[c,A,b,lu,Integers];
Print[x];
If[x[[i]]>0,
x=x*Apply[LCM,Denominator[x]];
dbgAssert[Apply[And,Map[(IntegerQ[#])&,{1,0,2}]]&&A.x===Table[0,{Length[A]}],latticePositiveElement,"x is not an integral element of the lattice."];
binomialSupport[x],
{}]*)
];

latticePositiveElement[LatticeByNullSpace[{{2,1,-1}}],proj[2,3],3]

latticePositiveElement[LatticeByNullSpace[{{2,1,1}}],proj[2,3],3]

latticeMaximalPositiveElements

Finds an integral vector in the lattice L = { v integer vector: A.v=0 } whose projection is non-negative with as many positive coordinates as possible.

SyntaxInformation[latticeMaximalPositiveElements]={"ArgumentsPattern"->{_,_}};latticeMaximalPositiveElements[L_List,projection_List]:=Module[{list,y,x,c,b,lu,n=Length[L[[1]]],i,coord},
y=binomialSupport[Table[0,{n}]];
list=Reap[For[i=1,i≤Length[projection[[ToricGB`Private`vectorList]]],i++,
coord=projection[[ToricGB`Private`vectorList,i]];
If[y[[ToricGB`Private`binomialExponent,coord]]>0,Continue[]];
x=latticePositiveElement[L,projection,coord];
If[x=!={},
Sow[x];
y=binomialAdd[y,x];
];
]][[2]];
If[list=!={},list=list[[1]]];
{y,list}
];

latticeMaximalPositiveElements[{{2,1,-1}},proj[3,3]]

maximalGrading

Computer a grading w.r.t. which the vectors in lattice L (given as lattice basis) are homogeneous and which has maximal support.
Therefore find gradings (nonnegative vectors) that are orthogonal to all vectors in L. To obtain maximal support, sum up gradings with different support.

SyntaxInformation[maximalGrading]={"ArgumentsPattern"->{_,_}};
maximalGrading[L_List,projection_List]:=Module[{grading,i,lu,n=Length[L[[1]]]},
grading=Table[0,{n}];
lu=Table[{0,0},{n}];
Map[(lu[[#]]={0,1})&,projection[[ToricGB`Private`vectorList]]];
For[i=1,i<=n,i++,
If[grading[[i]]>0||!MemberQ[projection[[ToricGB`Private`vectorList]],i],Continue[]];
grading+=LinearProgramming[-UnitVector[n,i],L,Table[{0,0},{Length[L]}],lu];
];
grading=grading*Apply[LCM,Denominator[grading]];
dbgAssert[L.grading===Table[0,{Length[L]}],maximalGrading,"L not homogeneous wrt. grading."];{grading,binomialSupportB[grading][[ToricGB`Private`binomialPosSupp]]}
];

maximalGrading[{{-1,0,2},{-1,1,1}},proj[3,3]]

Binomials

ToricGB`KDList=1;
ToricGB`KDTree=2;
ToricGB`MathList=1;

ToricGB`bitMasks=True;
ToricGB`dataStructure=ToricGB`KDList;
ToricGB`useToricGroebnerBasis=False;
ToricGB`intermediateReduce=True;

In most functions, binomials are stored with some extra information. The format is {v, marker, index, [posSupp, negSupp]}. v ist the exponent vector of the binomial, marker indicates whether it is active (1) or inactive (0). The index is only valid if the binomial is stored in G in geometricBuchberger and refers to the index of the binomial (needed for criterionMinimalSetQ). Otherwise index=0. The last two entries are only created/used in functions ending in B/BQ (stands for Bitmask).
A projection is usualy given as list with two elements. The first entry contains the list with the coordinate indices and the second list contains a bitmask of the indices.

Depending on ToricGB`dataStructure, lists of binomials are stored as kd-list, kd-trees or as list of a list of binomials (such that in both cases F[[ToricGB`Private`basisList]] is the (Mathematica) list of binomials.

ToricGB`Private`binomialExponent=1;
ToricGB`Private`binomialMarker=2;
ToricGB`Private`binomialIndex=3;
ToricGB`Private`binomialPosSupp=4;
ToricGB`Private`binomialNegSupp=5;
ToricGB`Private`vectorList=1;
ToricGB`Private`vectorSupport=2;

Set functions corresponding to flags

If[ToricGB`bitMasks,
binomialSupport:=binomialSupportB;
binomialPositivePart:=binomialPositivePartB;
binomialNegativePart:=binomialNegativePartB;
binomialNegate:=binomialNegateB;
binomialLCM:=binomialLCMB;
binomialAdd:=binomialAddB;
criterionSupportQ:=criterionSupportBQ;
criterionTailQ:=criterionTailBQ;
,
binomialSupport:=binomialSupportN;
binomialPositivePart:=binomialPositivePartN;
binomialNegativePart:=binomialNegativePartN;
binomialNegate:=binomialNegateN;
binomialLCM:=binomialLCMN;
binomialAdd:=binomialAddN;
criterionSupportQ:=criterionSupportNQ;
criterionTailQ:=criterionTailNQ;
];
Switch[ToricGB`dataStructure,ToricGB`KDList,
newList:=newKDList;
clearList:=clearKDList;
rebuildList:=rebuildKDList;
insertList:=insertKDList;
changeList:=changeKDList;
deleteList:=deleteKDList;
queryList:=queryKDList;
binomialReducesBy:=binomialReducesByKDL;
binomialReduce:=binomialReduceKDL;
criterionMinimalSetQ:=criterionMinimalSetKDLQ;

ToricGB`KDTree,
newList:=newKDTree;
clearList:=clearKDTree;
rebuildList:=rebuildKDTree;
insertList:=insertKDTree;
changeList:=changeKDTree;
deleteList:=deleteKDTree;
queryList:=queryKDTree;
binomialReducesBy:=binomialReducesByKDL;
binomialReduce:=binomialReduceKDL;
criterionMinimalSetQ:=criterionMinimalSetKDLQ;

ToricGB`MathList,
newList:=newMathList;
clearList:=clearMathList;
rebuildList:=rebuildMathList;
insertList:=insertMathList;
changeList:=changeMathList;
deleteList:=deleteMathList;
queryList:=queryMathList;
binomialReducesBy:=binomialReducesByVEC;
binomialReduce:=binomialReduceVEC;
criterionMinimalSetQ:=criterionMinimalSetVECQ;
];

positionsToBitMask

Given a list of indexes, create a bitmask representing the indexes and return the tuple of the list and the bitmask. The order of the bits has to agree with binomialSupport (see criterionSupportB), i.e. the first digit is the most significant bit.

SyntaxInformation[positionsToBitMask]={"ArgumentsPattern"->{_,_}};
positionsToBitMask[positions_List,n_Integer]:={positions,Total[2^(n-positions)]}

positionsToBitMask[{1,3,4},4]

binomialSupport

Given a vector v, return the triple of v, the marker (set to 1), the product of monomialOrder and v. In the version ending in B append the bitmasks of the support of v+ and the bitmask of the support of v- (N stands for No bitmasks). Also works for lists of vectors.

SyntaxInformation[binomialSupport]={"ArgumentsPattern"->{_}};
binomialSupportN[v_List]:=If[Head[v[[ToricGB`Private`binomialExponent]]]===List,Map[(binomialSupportN[#])&,v],{v,1,0}]
binomialSupportB[v_List]:= If[Head[v[[ToricGB`Private`binomialExponent]]]===List,Map[(binomialSupportB[#])&,v],{v,1,0,FromDigits[Clip[v,{0,1}],2],FromDigits[Clip[-v,{0,1}],2]}]

binomialSupportN[{2,0,-2,1}]

binomialSupportN[{{2,0},{0,-3}}]

binomialSupportB[{2,0,-2,1}]

binomialSupportB[{{2,0},{0,-3}}]

binomialStrip

Given a vector {v,marker,monomialOrder.v,[posSupp,negSupp]}, return v. Also works for lists of vectors.

SyntaxInformation[binomialStrip] = {"ArgumentsPattern" -> {_}};
binomialStrip[v_List] := If[Head[v[[ToricGB`Private`binomialExponent,ToricGB`Private`binomialExponent]]] === List, Map[(binomialStrip[#]) &, v],
  v[[ToricGB`Private`binomialExponent]]]

binomialStrip[binomialSupport[{2,0,-2,1}]]

binomialStrip[binomialSupportB[{{2,0},{0,-3}}]]

binomialPositivePart/binomialNegativePart

Given a vector v, return the vector v+ (resp. v-).

SyntaxInformation[binomialPositivePart]={"ArgumentsPattern"->{_}};
binomialPositivePartN[v_List]:={Clip[v[[ToricGB`Private`binomialExponent]],{0,Infinity}],v[[ToricGB`Private`binomialMarker]],0};binomialPositivePartB[v_List]:={Clip[v[[ToricGB`Private`binomialExponent]],{0,Infinity}],v[[ToricGB`Private`binomialMarker]],0,v[[ToricGB`Private`binomialPosSupp]],0};

SyntaxInformation[binomialNegativePart]={"ArgumentsPattern"->{_}};
binomialNegativePartN[v_List]:={Clip[-v[[ToricGB`Private`binomialExponent]],{0,Infinity}],v[[ToricGB`Private`binomialMarker]],0};binomialNegativePartB[v_List]:={Clip[-v[[ToricGB`Private`binomialExponent]],{0,Infinity}],v[[ToricGB`Private`binomialMarker]],0,v[[ToricGB`Private`binomialNegSupp]],0};

binomialPositivePartN[binomialSupportN[{-3,2,2}]]

binomialNegativePartN[binomialSupportN[{-3,2,2}]]

binomialPositivePartB[binomialSupportB[{-3,2,2}]]

binomialNegativePartB[binomialSupportB[{-3,2,2}]]

binomialNegate

Given v, return -v with corrected extra information.

SyntaxInformation[binomialNegate]={"ArgumentsPattern"->{_}};
binomialNegateN[v_List]:={-v[[ToricGB`Private`binomialExponent]],v[[ToricGB`Private`binomialMarker]],0}
binomialNegateB[v_List]:={-v[[ToricGB`Private`binomialExponent]],v[[ToricGB`Private`binomialMarker]],0,v[[ToricGB`Private`binomialNegSupp]],v[[ToricGB`Private`binomialPosSupp]]}

binomialNegateN[binomialSupportN[{3,-1,0}]]

binomialNegateB[binomialSupportB[{3,-1,0}]]

binomialAdd

Compute v+c*u with corrected extra information.

SyntaxInformation[binomialAdd]={"ArgumentsPattern"->{_,_,_.}};
binomialAddN[v_List,u_List,c_Integer:1]:={v[[ToricGB`Private`binomialExponent]]+c*u[[ToricGB`Private`binomialExponent]],1,0};
binomialAddB[v_List,u_List,c_Integer:1]:=binomialSupportB[v[[ToricGB`Private`binomialExponent]]+c*u[[ToricGB`Private`binomialExponent]]];

binomialAddN[binomialSupportN[{1,2,0}],binomialSupportN[{-1,2,1}]]

binomialAddN[binomialSupportN[{1,2,0}],binomialSupportN[{-1,2,1}],-1]

binomialAddB[binomialSupportB[{1,2,0}],binomialSupportB[{-1,2,1}]]

binomialAddB[binomialSupportB[{1,2,0}],binomialSupportB[{-1,2,1}],-1]

binomialLCM

Computes the least common multiple of (x^v+) and (x^u+) as vector, i.e. the entrywise maximum of v, u, and 0.

SyntaxInformation[binomialLCM]={"ArgumentsPattern"->{_,_}};
binomialLCMN[v_List,u_List]:=(dbgAssert[Length[v[[ToricGB`Private`binomialExponent]]]===Length[u[[ToricGB`Private`binomialExponent]]],binomialLCM,"Lengths of v and u differ."];binomialPositivePartN[{Array[(Max[v[[ToricGB`Private`binomialExponent,#]],u[[ToricGB`Private`binomialExponent,#]]])&,Length[v[[ToricGB`Private`binomialExponent]]]],1,0}])
binomialLCMB[v_List,u_List]:=(dbgAssert[Length[v[[ToricGB`Private`binomialExponent]]]===Length[u[[ToricGB`Private`binomialExponent]]],binomialLCM,"Lengths of v and u differ."];binomialPositivePartB[{Array[(Max[v[[ToricGB`Private`binomialExponent,#]],u[[ToricGB`Private`binomialExponent,#]]])&,Length[v[[ToricGB`Private`binomialExponent]]]],1,0,BitOr[v[[ToricGB`Private`binomialPosSupp]],u[[ToricGB`Private`binomialPosSupp]]],0}])

binomialLCMN[binomialSupportN[{3,-1,-4}],binomialSupportN[{2,4,-1}]]

binomialLCMB[binomialSupportB[{3,-1,-4}],binomialSupportB[{0,4,-1}]]

binomialZeroQ

Check whether projection (given as list of remaining indeces) of vector v is zero. If v==={} this indicates a zero vector, too.

SyntaxInformation[binomialZeroQ]={"ArgumentsPattern"->{_,_}};
binomialZeroQ[v_,projection_List]:=(v==={}||v[[ToricGB`Private`binomialExponent,projection[[ToricGB`Private`vectorList]]]]===Table[0,{Length[projection[[ToricGB`Private`vectorList]]]}]);

binomialZeroQ[binomialSupport[{0,3,0}],positionsToBitMask[{1,3},3]]

binomialZeroQ[binomialSupport[{0,3,0}],positionsToBitMask[{2,3},3]]

binomialEqualQ

Check whether the projections of the vectors v and u are the same.

SyntaxInformation[binomialEqualQ]={"ArgumentsPattern"->{_,_,_}};
binomialEqualQ[v_List,u_List,projection_List]:=(dbgAssert[Length[v[[ToricGB`Private`binomialExponent]]]===Length[u[[ToricGB`Private`binomialExponent]]],binomialEqualQ,"Lengths of v and u differ."];v[[ToricGB`Private`binomialExponent,projection[[ToricGB`Private`vectorList]]]]===u[[ToricGB`Private`binomialExponent,projection[[ToricGB`Private`vectorList]]]])

binomialEqualQ[binomialSupport[{1,2,3}],binomialSupport[{1,2,4}],proj[2,3]]

binomialDividesQ

Checks whether v+ divides u+ wrt. projection.

SyntaxInformation[binomialDividesQ]={"ArgumentsPattern"->{_,_,_,_.}};
(*binomialDividesQ[v_List,u_List,projection_List,startidx_Integer:1]:=(dbgAssert[Length[v[[ToricGB`Private`binomialExponent]]]==Length[u[[ToricGB`Private`binomialExponent]]],binomialDividesQ,"Lengths of v and u differ."];
Apply[And,Map[(v[[ToricGB`Private`binomialExponent,#]]≤Max[u[[ToricGB`Private`binomialExponent,#]],0])&,projection[[ToricGB`Private`vectorList]]]])*)
binomialDividesQ[v_List,u_List,projection_List,startidx_Integer:1]:=(binomialDivides[v,u,projection,startidx]>0)

binomialDividesQ[binomialSupport[{1,2,-3}],binomialSupport[{1,2,0}],proj[3,3]]

binomialDividesQ[binomialSupport[{1,3,1}],binomialSupport[{1,2,0}],proj[3,3]]

binomialDividesQ[binomialSupport[{2,3,-3}],binomialSupport[{1,-2,0}],proj[3,3],2]

binomialDivides

Computes which power of v+ divides u+ wrt. projection.

SyntaxInformation[binomialDivides]={"ArgumentsPattern"->{_,_,_,_.}};
binomialDivides[v_List,u_List,projection_List,startidx_Integer:1]:=Module[{i,coord,exponent=Infinity},
dbgAssert[Length[v[[ToricGB`Private`binomialExponent]]]==Length[u[[ToricGB`Private`binomialExponent]]],binomialDivides,"Lengths of v and u differ."];
For[i=startidx,i≤Length[projection[[ToricGB`Private`vectorList]]],i++,
cntDividesIter++;
coord=projection[[ToricGB`Private`vectorList,i]];
If[v[[ToricGB`Private`binomialExponent,coord]]≤0,Continue[]];
If[v[[ToricGB`Private`binomialExponent,coord]]>u[[ToricGB`Private`binomialExponent,coord]],Return[0]];
exponent=Min[exponent,Quotient[u[[ToricGB`Private`binomialExponent,coord]],v[[ToricGB`Private`binomialExponent,coord]]]];
];
exponent
]

binomialDivides[binomialSupport[{1,2,-3}],binomialSupport[{3,4,0}],proj[3,3]]

binomialDivides[binomialSupport[{1,3,1}],binomialSupport[{1,2,0}],proj[3,3]]

binomialOrder

Returns 1 if v+ > u+ wrt. monomialOrder, given as weight matrix, and projection, -1 if v+ < u+ and 0 otherwise.

SyntaxInformation[binomialOrder]={"ArgumentsPattern"->{_,_,_,_}};
binomialOrder[v_List,u_List,monomialOrder_List,projection_List]:=Module[{i,j,coord,x},
dbgAssert[Length[v[[ToricGB`Private`binomialExponent]]]===Length[u[[ToricGB`Private`binomialExponent]]],binomialOrder,"Lengths of v and u differ."];
For[i=1,i<=Length[monomialOrder],i++,
dbgAssert[Length[v[[ToricGB`Private`binomialExponent]]]===Length[monomialOrder[[i]]],binomialOrder,"monomialOrder has wrong number of columns."];
For[x=0;j=1,j<=Length[projection[[ToricGB`Private`vectorList]]],j++,
coord=projection[[ToricGB`Private`vectorList,j]];
If[v[[ToricGB`Private`binomialExponent,coord]]>0,x+=v[[ToricGB`Private`binomialExponent,coord]]*monomialOrder[[i,coord]]];
If[u[[ToricGB`Private`binomialExponent,coord]]>0,x-=u[[ToricGB`Private`binomialExponent,coord]]*monomialOrder[[i,coord]]];
];
If[x≠0,Return[If[x>0,1,-1]]];
];
0
];

binomialOrder[binomialSupport[{3,5,0}],binomialSupport[{3,4,0}],lexOrder[3],proj[3,3]]

binomialOrder[binomialSupport[{-3,5,-1}],binomialSupport[{-4,5,1}],lexOrder[3],proj[3,3]]

binomialOrder[binomialSupport[{3,5,-5}],binomialSupport[{3,4,-3}],lexOrder[3],positionsToBitMask[{1,3},3]]

binomialStandardFormQ

Returns True if v+ >= v- wrt. monomialOrder and projection.

SyntaxInformation[binomialStandardFormQ]={"ArgumentsPattern"->{_,_,_}};
binomialStandardFormQ[v_List,monomialOrder_List,projection_List]:=(binomialOrder[v,binomialNegate[v],monomialOrder,projection]>=0)

binomialStandardFormQ[binomialSupport[{3,-4,0}],lexOrder[3],proj[3,3]]

binomialStandardFormQ[binomialSupport[{3,-4,0}],lexOrder[3],positionsToBitMask[{2,3},3]]

binomialStandardForm

Returns u in {v, -v} such that isBinomialStandardForm[u], i.e. u+ > u- wrt. monomialOrder and projection.

SyntaxInformation[binomialStandardForm]={"ArgumentsPattern"->{_,_,_}};
binomialStandardForm[v_List,monomialOrder_List,projection_List]:=If[binomialStandardFormQ[v,monomialOrder,projection],v,binomialNegate[v]]

binomialStandardForm[binomialSupport[{3,-4,0}],lexOrder[3],positionsToBitMask[{2,3},3]]

binomialReducesBy

Internal`binomialHead=1;
Internal`binomialTail=2;

Check whether binomial v reduces with respect to binomials in F which have to be in standard form wrt. the desired monomialOrder and projection. monomialSelection controls, whether the head of the binomial, the tail or both are checked. startidx indicates that only the part of F starting with this index has to be checked. The binomials in F contain markers, i.e. are of the form {u,1} if active and {u,0} if inactive. The return value is the index of the first reducer binomial in F or 0 if none is reducable. The sign of the return value indicates whether it reduces v+ (positive) or v- (negative).

SyntaxInformation[binomialReducesBy]={"ArgumentsPattern"->{_,_,_,_,_.,_.}};

binomialReducesByVEC[F_List,v_List, monomialOrder_List,monomialSelection_Integer,startidx_Integer:1,endidx_Integer:-1]:=Module[{u=binomialNegate[v],i,end=Length[F[[ToricGB`Private`basisList]]]},
If[endidx≥0,end=Min[end,endidx]];
dbgLog["FindLessQ[",v,"]"];
cntQuery++;
For[i=startidx,i≤end,i++,
If[F[[ToricGB`Private`basisList,i,ToricGB`Private`binomialMarker]]=!=1,Continue[]];
If[BitAnd[monomialSelection,Internal`binomialHead]≠0,
cntReduceQDivQ++;
If[binomialDividesQ[F[[ToricGB`Private`basisList,i]],v,F[[ToricGB`Private`basisProjection]]],Return[i]]];
If[BitAnd[monomialSelection,Internal`binomialTail]≠0,
cntReduceQDivQ++;
If[binomialDividesQ[F[[ToricGB`Private`basisList,i]],u,F[[ToricGB`Private`basisProjection]]],Return[-i]]];
];
0
]

SyntaxInformation[callbackReducesByKDL]={"ArgumentsPattern"->{_,_,_}};
SetAttributes[callbackReducesByKDL,HoldFirst];
callbackReducesByKDL[F_Symbol,query_List,v_List]:=If[F[[ToricGB`Private`basisList,v[[ToricGB`Private`binomialIndex]],ToricGB`Private`binomialMarker]]=!=0,v[[ToricGB`Private`binomialIndex]]]

SetAttributes[binomialReducesByKDL,HoldFirst];
binomialReducesByKDL[F_Symbol,v_List, monomialOrder_List,monomialSelection_Integer]:=Module[{i},
If[BitAnd[monomialSelection,Internal`binomialHead]≠0,
i=queryList[F,binomialPositivePart[v],1,callbackReducesByKDL];
If[i=!={},Return[i[[1,1]]]];
];
If[BitAnd[monomialSelection,Internal`binomialTail]≠0,
i=queryList[F,binomialNegativePart[v],1,callbackReducesByKDL];
If[i=!={},Return[-i[[1,1]]]];
];
0
];

binomialReduce

Reduces binomial v with respect to binomials in F which have to be in standard form wrt. the desired monomialOrder and projection. If grading is not Null, apply criterionTailQ. monomialSelection controls, whether the head of the binomial, the tail or both are reduced. startidx indicates that only the part of F starting with this index has to be checked the first time. The binomials in F contain markers, i.e. are of the form {u,1} if active and {u,0} if inactive.

SyntaxInformation[binomialReduce]={"ArgumentsPattern"->{_,_,_,_,_,_.}};

binomialReduceVEC[F_List,v_List,monomialOrder_List,grading_List,monomialSelection_Integer,startidx_Integer:1]:=Module[{u=v,w=binomialNegate[v],i=startidx,j,c},
dbgPrint[ToricGB`Debug`binomialReduce,"binomialReduceVEC[",F,", ",v,", ",monomialOrder,", ",grading,", ",monomialSelection,", ",startidx];
dbgPrint[ToricGB`Debug`binomialReduce,"FindLess[",u,"]"];
For[j=startidx,j≤Length[F[[ToricGB`Private`basisList]]],j++;i=Mod[i,Length[F[[ToricGB`Private`basisList]]]]+1,
If[BitAnd[monomialSelection,Internal`binomialHead]≠0,
cntReduceQDivQ++;
If[binomialDividesQ[F[[ToricGB`Private`basisList,i]],u,F[[ToricGB`Private`basisProjection]]],
cntReduceDiv++;
If[grading[[ToricGB`Private`vectorSupport]]≠0&&criterionTailQ[u,F[[ToricGB`Private`basisList,i]],grading,F[[ToricGB`Private`basisProjection]]],dbgPrint[ToricGB`Debug`geometricBuchberger,"Drop ",v," during binomialReduce because of criterionTail."];Return[{}]];
c=binomialDivides[F[[ToricGB`Private`basisList,i]],u,F[[ToricGB`Private`basisProjection]]];
u=binomialAdd[u,F[[ToricGB`Private`basisList,i]],-c];
If[binomialZeroQ[u,F[[ToricGB`Private`basisProjection]]],Return[{}]];
w=binomialNegate[u];
j=0;
dbgPrint[ToricGB`Debug`binomialReduce,"FindLess[",u,"]"];
];
];
If[BitAnd[monomialSelection,Internal`binomialTail]≠0,
cntReduceQDivQ++;
If[binomialDividesQ[F[[ToricGB`Private`basisList,i]],w,F[[ToricGB`Private`basisProjection]]],
cntReduceDiv++;
If[grading[[ToricGB`Private`vectorSupport]]≠0&&criterionTailQ[w,F[[ToricGB`Private`basisList,i]],grading,F[[ToricGB`Private`basisProjection]]],dbgPrint[ToricGB`Debug`geometricBuchberger,"Drop ",v," during binomialReduce because of criterionTail."];Return[{}]];
c=binomialDivides[F[[ToricGB`Private`basisList,i]],w,F[[ToricGB`Private`basisProjection]]];
u=binomialAdd[u,F[[ToricGB`Private`basisList,i]],c];
If[binomialZeroQ[u,F[[ToricGB`Private`basisProjection]]],Return[{}]];
w=binomialNegate[u];
j=0;
dbgPrint[ToricGB`Debug`binomialReduce,"FindLess[",u,"]"];
];
];
];
binomialStandardForm[u,monomialOrder,F[[ToricGB`Private`basisProjection]]]
];

SyntaxInformation[callbackReduceKDL]={"ArgumentsPattern"->{_,_,_}};
SetAttributes[callbackReduceKDL,HoldFirst];
callbackReduceKDL[F_Symbol,query_List,v_List]:=(If[F[[ToricGB`Private`basisList,v[[ToricGB`Private`binomialIndex]],ToricGB`Private`binomialMarker]]=!=0,v])

SyntaxInformation[binomialReduceKDL]={"ArgumentsPattern"->{_,_,_,_,_}};
SetAttributes[binomialReduceKDL,HoldFirst];
binomialReduceKDL[F_Symbol,v_List,monomialOrder_List,grading_List,monomialSelection_Integer]:=Module[{u=v,w,changed=True},
While[changed,
changed=False;
dbgPrint[ToricGB`Debug`binomialReduce,"FindLess[",u,"]"];
If[BitAnd[monomialSelection,Internal`binomialHead]≠0,
w=queryList[F,binomialPositivePart[u],1,callbackReduceKDL,False,True];
If[w=!={},
w=w[[1,1]];
dbgPrint[ToricGB`Debug`binomialReduce,"-> ",w];
If[grading[[ToricGB`Private`vectorSupport]]≠0&&criterionTailQ[u,w,grading,F[[ToricGB`Private`basisProjection]]],dbgPrint[ToricGB`Debug`geometricBuchberger,"Drop ",v," during binomialReduce because of criterionTail."];Return[{}]];
dbgLog[u,w,binomialDivides[w,u,F[[ToricGB`Private`basisProjection]]],binomialAdd[u,w,-binomialDivides[w,u,F[[ToricGB`Private`basisProjection]]]]];
u=binomialAdd[u,w,-binomialDivides[w,u,F[[ToricGB`Private`basisProjection]]]];
If[binomialZeroQ[u,F[[ToricGB`Private`basisProjection]]],Return[{}]];
changed=True]];
If[BitAnd[monomialSelection,Internal`binomialTail]≠0,
w=queryList[F,binomialNegativePart[u],1,callbackReduceKDL,False,True];
If[w=!={},
w=w[[1,1]];
dbgPrint[ToricGB`Debug`binomialReduce,"-> ",w];
If[grading[[ToricGB`Private`vectorSupport]]!=0&&criterionTailQ[binomialNegate[u],w,grading,F[[ToricGB`Private`basisProjection]]],dbgPrint[ToricGB`Debug`geometricBuchberger,"Drop ",v," during binomialReduce because of criterionTail."];Return[{}]];
dbgLog[u,w,binomialDivides[w,-u,F[[ToricGB`Private`basisProjection]]],binomialAdd[u,w,binomialDivides[w,-u,F[[ToricGB`Private`basisProjection]]]]];
u=binomialAdd[u,w,binomialDivides[w,-u,F[[ToricGB`Private`basisProjection]]]];
If[binomialZeroQ[u,F[[ToricGB`Private`basisProjection]]],Return[{}]];
changed=True]];
];
binomialStandardForm[u,monomialOrder,F[[ToricGB`Private`basisProjection]]]
]

binomialReduceVEC[{binomialSupport[{{1,-1,-1},{0,1,-1}}],proj[3,3]},binomialSupport[{1,1,0}],lexOrder[3],{{0,0,0},0},3]

F=newKDList[proj[3,3]];
insertKDList[F,binomialSupport[{1,-1,-1}]];
insertKDList[F,binomialSupport[{0,1,-1}]];
binomialReduceKDL[F,binomialSupport[{1,1,0}],lexOrder[3],{{0,0,0},0},3]
binomialReducesByKDL[F,binomialSupport[{1,1,0}],lexOrder[3],3]

F={binomialSupport[{{1,-1,-1},-{0,1,-1}}],positionsToBitMask[{1,3},3]};
binomialReduceVEC[F,binomialSupport[{1,1,0}],lexOrder[3],{{0,0,0},0},3]

F=newKDList[proj[1,3]];
insertKDList[F,binomialSupport[{1,-1,-1}]];
insertKDList[F,binomialSupport[-{0,1,-1}]];binomialReduceKDL[F,binomialSupport[{1,1,0}],lexOrder[3],{{0,0,0},0},3]

VectorToBinomial

VectorToBinomial converts the vector v to the binomial (x^v+) - (x^v-) where v+ contains the positive entries of v and v- = (-v)+. Also applicable to list of vectors.

SyntaxInformation[VectorToBinomial]={"ArgumentsPattern"->{_,_}};VectorToBinomial[v_List,variables_List]:=If[Head[v[[1]]]===List,Map[(VectorToBinomial[#,variables])&,v],
assert[Length[v]===Length[variables],VectorToBinomial,"Lengths of v and variables differ."];Numerator[Together[Apply[Times,variables^v]-1]]];

VectorToBinomial[{1,2,0},vars[3]]

VectorToBinomial[{{1,2,-4},{3,-6,1},{1,0,-1}},vars[3]]

BinomialToVector

Provides the inverse function to VectorToBinomial. Might inverse the sign, though.

SyntaxInformation[BinomialToVector]={"ArgumentsPattern"->{_,_}};BinomialToVector[F_List,variables_List]:=Map[(BinomialToVector[#,variables])&,F];
BinomialToVector[p_,variables_List]:=Module[{v=Table[0,{Length[variables]}],q,r,c,e},
q=Collect[p,variables];
assert[Head[q]===Plus,BinomialToVector,"Only binomials may be given as input."];
r=First[q]; q=Rest[q];
assert[Head[q]=!=Plus,BinomialToVector,"Only binomials may be given as input."];
e=Exponent[r,variables];
c=Extract[CoefficientList[r,variables],e+1];
assert[c===1||c===-1,BinomialToVector,"All coefficients must be +-1."];
v+=e;
q*=c;
e=Exponent[q,variables];
c=Extract[CoefficientList[q,variables],e+1];
assert[c===-1,BinomialToVector,"The second coefficient of the binomial must have the opposite sign of the first coefficient."];
v-e
];

BinomialToVector[VectorToBinomial[{1,2,0},vars[3]],vars[3]]

BinomialToVector[VectorToBinomial[{{1,2,-4},{3,-6,1},{1,0,-1}},vars[3]],vars[3]]

Criteria for Buchberger's Algorithm

ToricGB`Private`criterionSupport=1;
ToricGB`Private`criterionMinimalSet=2;
ToricGB`Private`criterionTail=4;

criterionSupportQ

A pair {v,u} can be dropped if the intersection of the supports of u+ and v+ (w.r.t. projection) is empty. This is Buchbergers well-known criterion reformulated for vectors (see Hemmecke-Malkin).

SyntaxInformation[criterionSupportQ]={"ArgumentsPattern"->{_,_,_}};criterionSupportNQ[v_List, u_List, projection_List]:=Apply[And,Map[((v[[ToricGB`Private`binomialExponent,#]] <= 0||u[[ToricGB`Private`binomialExponent,#]] <= 0))&,projection[[ToricGB`Private`vectorList]]]]
criterionSupportBQ[v_List, u_List, projection_List]:=(BitAnd[v[[ToricGB`Private`binomialPosSupp]],u[[ToricGB`Private`binomialPosSupp]],projection[[ToricGB`Private`vectorSupport]]]==0)

criterionSupportNQ[binomialSupportN[{0,1,0,1}],binomialSupportN[{4,0,-4,1}],proj[4,4]]

criterionSupportNQ[binomialSupportN[{1,1,0,1}],binomialSupportN[{1,0,0,0}],positionsToBitMask[{2,3,4},4]]

criterionSupportBQ[binomialSupportB[{0,1,0,1}],binomialSupportB[{4,0,-4,1}],proj[4,4]]

criterionSupportBQ[binomialSupportB[{1,1,0,1}],binomialSupportB[{1,0,0,0}],positionsToBitMask[{2,3,4},4]]

criterionMinimalSetQ

It suffices to consider pairs that generate the module of syzygies. criterionMinimalSet drops pairs to get a pretty good approximation. It was originally given by Gebauer and Möller and reformulated to this form by Hemmecke and Malkin. The binomials in F contain markers, i.e. are of the form {u,1} if active and {u,0} if inactive.

SyntaxInformation[criterionMinimalSetQ]={"ArgumentsPattern"->{_,_,_}};

criterionMinimalSetVECQ[F_List,i_Integer,j_Integer]:=Module[{w,k},
dbgAssert[i<j,criterionMinimalSetQ,"i must be smaller than j."];
cntQuery++;
w=binomialLCM[F[[ToricGB`Private`basisList,i]],F[[ToricGB`Private`basisList,j]]];
For[k=1,k<i,k++,cntMSDivQ++;If[F[[ToricGB`Private`basisList,k,ToricGB`Private`binomialMarker]]===1&&binomialDividesQ[F[[ToricGB`Private`basisList,k]],w,F[[ToricGB`Private`basisProjection]]],Return[True]]];
For[k=i+1,k<j,k++,cntMSDivQ++;If[F[[ToricGB`Private`basisList,k,ToricGB`Private`binomialMarker]]===1&&binomialDividesQ[F[[ToricGB`Private`basisList,k]],w,F[[ToricGB`Private`basisProjection]]]&&!binomialEqualQ[binomialLCM[F[[ToricGB`Private`basisList,k]],F[[ToricGB`Private`basisList,j]]],w,F[[ToricGB`Private`basisProjection]]],Return[True]]];
For[k=j+1,k≤Length[F[[ToricGB`Private`basisList]]],k++,cntMSDivQ++;If[F[[ToricGB`Private`basisList,k,ToricGB`Private`binomialMarker]]===1&&binomialDividesQ[F[[ToricGB`Private`basisList,k]],w,F[[ToricGB`Private`basisProjection]]]&&!binomialEqualQ[binomialLCM[F[[ToricGB`Private`basisList,k]],F[[ToricGB`Private`basisList,i]]],w,F[[ToricGB`Private`basisProjection]]]&&!binomialEqualQ[binomialLCM[F[[ToricGB`Private`basisList,k]],F[[ToricGB`Private`basisList,j]]],w,F[[ToricGB`Private`basisProjection]]],Return[True]]];
False
];

SetAttributes[callbackMinimalSetKDL,HoldFirst];
callbackMinimalSetKDL[F_Symbol,query_List,v_List,i_Integer,j_Integer]:=
(If[F[[ToricGB`Private`basisList,v[[ToricGB`Private`binomialIndex]],ToricGB`Private`binomialMarker]]===1,
If[v[[ToricGB`Private`binomialIndex]]<i,True,
If[i<v[[ToricGB`Private`binomialIndex]]&&v[[ToricGB`Private`binomialIndex]]<j,If[!binomialEqualQ[binomialLCM[v,F[[ToricGB`Private`basisList]][[j]]],query,F[[ToricGB`Private`basisProjection]]],True],
If[j<v[[ToricGB`Private`binomialIndex]],If[!binomialEqualQ[binomialLCM[v,F[[ToricGB`Private`basisList]][[i]]],query,F[[ToricGB`Private`basisProjection]]]&&!binomialEqualQ[binomialLCM[v,F[[ToricGB`Private`basisList]][[j]]],query,F[[ToricGB`Private`basisProjection]]],True]
]
]
]
])

SetAttributes[criterionMinimalSetKDLQ,HoldFirst];
criterionMinimalSetKDLQ[F_Symbol,i_Integer,j_Integer]:=Module[{w},
dbgAssert[i<j,criterionMinimalSetQ,"i must be smaller than j."];
w=binomialLCM[F[[ToricGB`Private`basisList]][[i]],F[[ToricGB`Private`basisList]][[j]]];
queryList[F,w,1,Function[{KD,query,it},callbackMinimalSetKDL[KD,query,it,i,j],{HoldAll}]]=!={}
];

criterionTailQ

This criterion is specific to homogeneous ideals (w.r.t. some grading) that are saturated w.r.t. the product of the variables. If the S-polynomial spol reduces to some polynomial recuct of lower grading, the S-polynomial can be dropped. When called by geometricBuchberger, it checks, whether the reduction of the actual S-polynomial to v-u already is a decrease of the grading. This function is also called by monomialReduce. Here it checks, whether the reduction of u by v to u-v (or vice versa) decreases the grading. This criterion is special to toric ideals and was given by Hemmecke and Malkin.

SyntaxInformation[criterionTailQ]={"ArgumentsPattern"->{_,_,_,_}};
criterionTailNQ[u_List,v_List,grading_List,projection_List]:=(dbgAssert[Length[v[[ToricGB`Private`binomialExponent]]]===Length[u[[ToricGB`Private`binomialExponent]]],criterionTailQ,"Lengths of v and u differ."];
Apply[Or,Map[((v[[ToricGB`Private`binomialExponent,#]]<0&&u[[ToricGB`Private`binomialExponent,#]]<0&&grading[[ToricGB`Private`vectorList,#]]>0))&,projection[[ToricGB`Private`vectorList]]]])
criterionTailBQ[u_List,v_List,grading_List,projection_List]:=(dbgAssert[Length[v[[ToricGB`Private`binomialExponent]]]===Length[u[[ToricGB`Private`binomialExponent]]],criterionTailQ,"Lengths of v and u differ."];
BitAnd[v[[ToricGB`Private`binomialNegSupp]],u[[ToricGB`Private`binomialNegSupp]],grading[[ToricGB`Private`vectorSupport]],projection[[ToricGB`Private`vectorSupport]]]≠0)

grad1={{1,2,0},binomialSupportB[{1,2,0}][[ToricGB`Private`binomialPosSupp]]};
criterionTailNQ[binomialSupportN[{-1,0,-1}],binomialSupportN[{-2,0,0}],grad1,proj[3,3]]

grad2={{0,1,1},binomialSupportB[{0,1,1}][[ToricGB`Private`binomialPosSupp]]};criterionTailNQ[binomialSupportN[{-1,-1,0}],binomialSupportN[{-1,0,-1}],grad2,proj[3,3]]

criterionTailBQ[binomialSupportB[{-1,0,-1}],binomialSupportB[{-2,0,0}],grad1,proj[3,3]]

criterionTailBQ[binomialSupportB[{-1,-1,0}],binomialSupportB[{-1,0,-1}],grad2,proj[3,3]]

Algorithms

latticePositiveElements

(*latticeDegree[v_List,k_Integer]:=Sum[v[[i]],{i,k}];*)

(*liftLatticeElement[v_List,L_List]:=Module[{B,k,i,u},
B=Map[Reverse,HermiteDecomposition[Map[Reverse,L]][[2]]];
Print[B];
k=Length[v];
u=LinearSolve[Transpose[B[[1;;k,1;;k]]],v];
Print[u];
ClearAll[x];
u=Append[v,x];
u[[k]]=Minimize[{x,x≥0&&A[[All,1;;k]].u==0},x,Integers][[1]];
u
];*)

(*liftLatticeElement[{2,0},{{1,0,1},{3,0,1}}]*)

(*latticePositiveElements[A_List,L_List]:=Module[{B,i,j,k,F,G,n=Length[L[[1]]],u,v,spol,pairs,maxdeg,deg,curdeg,coord},
B=HermiteDecomposition[L][[2]];
Print[B];
F={{B[[1,1]]}};
For[k=2,k<=n,k++,
G=Map[(liftLatticeElement[#,A])&,F];
If[k<=Length[B],G=Join[G,{B[[k,1;;k]],-B[[k,1;;k]]}]];
Print[k,G];
coord=Table[x,{x,k-1}];
ClearAll[pairs];
maxdeg=0;
For[j=1,j≤Length[G],j++,
For[i=1,i<j,i++,
If[G[[i,k]]*G[[j,k]]≥0,Continue[]];
spol=G[[i,k]]+G[[j,k]];
deg=latticeDegree[spol,k-1];
maxdeg=Max[deg,maxdeg];
If[Head[pairs[deg]]=!=List,pairs[deg]={spol},pairs[deg]=Append[pairs[deg],spol]];
];
];
For[curdeg=1,curdeg≤maxdeg,curdeg++,
For[i=1,i≤Length[pairs[deg]],i++,
spol=pairs[deg][[i]];
For[j=1,j≤Length[G],j++,
If[Apply[And,Thread[G[[j,coord]]≤spol[[coord]]]]&&Abs[G[[j,k]]]≤Abs[spol[[k]]]&&G[[j,k]]*spol[[k]]>=0,Break[]];
];
If[j>Length[G],
G=Append[G,spol];
For[j=1,j<Length[G],j++,
If[G[[j,k]]*G[[Length[G],k]]≥0,Continue[]];
spol=G[[j,k]]+G[[Length[G],k]];
deg=latticeDegree[spol,k-1];
maxdeg=Max[deg,maxdeg];
If[Head[pairs[deg]]=!=List,pairs[deg]={spol},pairs[deg]=Append[pairs[deg],spol]];
];
];
];
];
];
];*)

(*latticePositiveElements[A_List]:=Module[{B,i,j,k,F,G,n=Length[L[[1]]],u,v,spol,pairs,maxdeg,deg,curdeg,coord},
];*)

(*A={{1,-3,4}};
latticePositiveElements[A,LatticeByNullSpace[A]];*)

intermediateReduce

During the Buchberger algorithm sometimes one gets cascading sequences of polynomials where the later ones reduce the previous ones. In order to keep the intermediate bases small, an interreduction is performed. In this implementation, intermediateReduce is called for each new polynomial. Instead of reducing each polynomial wrt. all other polynomials, a quick scan is performed for polynomials that are reducible by the new polynomial. For those a full reduction is done. They are then treated as new and therefore polynomials which are reducible by them are searched.

SetAttributes[callbackIntermediateReduce,HoldFirst];
callbackIntermediateReduce[F_Symbol,query_List,v_List]:=(If[F[[ToricGB`Private`basisList,v[[ToricGB`Private`binomialIndex]],ToricGB`Private`binomialMarker]]=!=0,v[[ToricGB`Private`binomialIndex]]]);

SyntaxInformation[intermediateReduce]={"ArgumentsPattern"->{_,_,_,_}};
SetAttributes[intermediateReduce,HoldFirst];
intermediateReduce[F_Symbol,min_Integer,new_Integer,monomialOrder_List]:=Module[{changed=newQueue[],results,i,j,u,v,e,cur,nullGrading={Table[0,{Length[F[[ToricGB`Private`basisList,1]]]}],0}},
dbgPrint[ToricGB`Debug`geometricBuchberger,"intermediateReduce: F=",F[[ToricGB`Private`basisList]],", min=",min,", new=",new,", monomialOrder=",monomialOrder,", projection=",F[[ToricGB`Private`basisProjection]]];
enQueue[changed,{new,F[[ToricGB`Private`basisList,new]]}];
While[!queueEmpty[changed],
{i,u}=deQueue[changed];
If[F[[ToricGB`Private`basisList,i]]=!=u,Continue[]];
(*u=F[[ToricGB`Private`basisList,new]]; i=new; *)
results=queryList[F,u,Infinity,callbackIntermediateReduce,True];
If[results=!={},results=results[[1]]];
For[j=1,j≤Length[results],j++,
If[results[[j]]==i,Continue[]];
cntReduceDiv++;
e=binomialDivides[u,F[[ToricGB`Private`basisList,results[[j]]]],F[[ToricGB`Private`basisProjection]]];
dbgAssert[e>0,intermediateReduce,"Query returned element ",F[[ToricGB`Private`basisList,results[[j]]]]," not divivisible by ",u,". F=",F[[ToricGB`Private`basisList]],", results=",results];
v=binomialAdd[F[[ToricGB`Private`basisList,results[[j]]]],u,-e];
F[[ToricGB`Private`basisList,results[[j]],ToricGB`Private`binomialMarker]]=0;
v=binomialReduce[F,v,monomialOrder,nullGrading,BitOr[Internal`binomialHead,Internal`binomialTail]];
dbgPrint[ToricGB`Debug`geometricBuchberger,F[[ToricGB`Private`basisList,results[[j]]]]," reduces to ",v];
If[!binomialZeroQ[v,F[[ToricGB`Private`basisProjection]]],
If[results[[j]]≤min,
deleteList[F,results[j]];
insertList[F,v];
enQueue[changed,{Length[F[[ToricGB`Private`basisList]]],F[[ToricGB`Private`basisList,Length[F[[ToricGB`Private`basisList]]]]]}];
,
changeList[F,results[[j]],v];
enQueue[changed,{results[[j]],F[[ToricGB`Private`basisList,results[[j]]]]}];
];
];
];
];
clearQueue[changed];
];

geometricBuchberger

Performs the Buchberger algorithm for the toric ideal I=VectorToBinomial[F] doing all computations with vectors (of the exponents). The result is a Gröbner basis (in vector form) of an ideal containing I and contained in the radical of I.
In the main loop, for all pairs {v,u} in the approximation to the Gröbner basis G the S-polynomial S(v,u)=v-u is calculated and reduced w.r.t. G. If the remainder is non-zero, it is added to G. Different criteria tell which pairs {v,u} may be omitted.
There are several criteria for dropping pairs which can be controlled by the bitmask criteria.
Finally the Gröbner basis is reduced. Head-reducible binomials are superfluous.

SyntaxInformation[geometricBuchberger]={"ArgumentsPattern"->{_,_,_,_,_}};geometricBuchberger[F_List, monomialOrder_List, gradingP_List,projection_List,criteria_Integer] :=Module[{G=newList[projection],H,spol,i, j,k,l,v,u,grading=gradingP,n=Length[F[[1]]]},
If[BitAnd[criteria,ToricGB`Private`criterionTail]==0,grading={Table[0,{n}],0}];
dbgPrint[ToricGB`Debug`geometricBuchberger,"geometricBuchberger[",F,", ",monomialOrder,", ",grading,", ",projection,", ",criteria,"]"];
dbgLog["geometricBuchberger[",F,", ",monomialOrder,", ",grading,", ",projection,", ",criteria,"]"];
For[i=1,i<=Length[F],i++,
cntInsert++;
insertList[G,binomialStandardForm[F[[i]],monomialOrder,projection]];
];
timePrereduce+=AbsoluteTiming[For[i=1,i<=Length[G[[ToricGB`Private`basisList]]],i++,
If[G[[ToricGB`Private`basisList,i,ToricGB`Private`binomialMarker]]===0,Continue[]];
cntReduce0++;
intermediateReduce[G,0,i,monomialOrder];
]][[1]];
dbgPrint[ToricGB`Debug`geometricBuchberger,"Prereduced: G=",G[[ToricGB`Private`basisList]]];
timeGroebner+=AbsoluteTiming[
For[j=1,j≤Length[G[[ToricGB`Private`basisList]]],j++,
If[G[[ToricGB`Private`basisList,j,ToricGB`Private`binomialMarker]]===0,Continue[]];
j=addTime[rebuildList[G,j],timeRebuild];
For[i=1,i<j,i++,
If[G[[ToricGB`Private`basisList,i,ToricGB`Private`binomialMarker]]===0,Continue[]];
v=G[[ToricGB`Private`basisList,i]]; u=G[[ToricGB`Private`basisList,j]];
cntCS++;
If[BitAnd[criteria,ToricGB`Private`criterionSupport]≠0&&criterionSupportQ[v,u,projection],dbgPrint[ToricGB`Debug`geometricBuchberger,"Drop ",{i,j}," because of criterionSupportQ."];Continue[]];
cntCT++;
If[BitAnd[criteria,ToricGB`Private`criterionTail]≠0&&criterionTailQ[v,u,grading,projection],dbgPrint[ToricGB`Debug`geometricBuchberger,"Drop ",{i,j}," because of criterionTailQ."];Continue[]];
cntCM++;
If[BitAnd[criteria,ToricGB`Private`criterionMinimalSet]≠0&&addTime[criterionMinimalSetQ[G,i,j],timeCM],dbgPrint[ToricGB`Debug`geometricBuchberger,"Drop ",{i,j}," because of criterionMinimalSetQ."];Continue[]];
spol=binomialAdd[v,u,-1];
cntReduce1++;
spol=binomialReduce[G,spol,monomialOrder,grading,BitOr[Internal`binomialHead,Internal`binomialTail]];
If[!binomialZeroQ[spol,projection],
cntInsert++;
insertList[G,spol];
If[ToricGB`intermediateReduce,addTime[intermediateReduce[G,j+1,Length[G[[ToricGB`Private`basisList]]],monomialOrder],timeInterRed]];
];
dbgPrint[ToricGB`Debug`geometricBuchberger,"Pair ",{i,j},": ",binomialAdd[v,u,-1]," reduces to ",spol];
dbgPrint[ToricGB`Debug`geometricBuchberger,"G=",G[[ToricGB`Private`basisList]]];
];
]][[1]];
If[ToricGB`Debug`profile,maxExponent=Map[(Max[#])&,Abs[Transpose[binomialStrip[G[[ToricGB`Private`basisList]]]]]];sizeG=Length[G[[ToricGB`Private`basisList]]]];
(* reduce Gröbner basis *)
dbgPrint[ToricGB`Debug`geometricBuchberger,"Reduce..."];
timeRedGB+=AbsoluteTiming[H=Reap [For[i=1,i≤Length[G[[ToricGB`Private`basisList]]],i++,
If[G[[ToricGB`Private`basisList,i,ToricGB`Private`binomialMarker]]===0,Continue[]];
v=G[[ToricGB`Private`basisList,i]];
G[[ToricGB`Private`basisList,i,ToricGB`Private`binomialMarker]]=0;
If[(*True||*)!ToricGB`intermediateReduce,
cntReduce2Q++;
If[binomialReducesBy[G,v,monomialOrder,Internal`binomialHead]≠0,
dbgPrint[ToricGB`Debug`geometricBuchberger,"Drop G[[",i,"]]=",v," because of binomialReduceQ"];
Continue[]];
];
cntReduce2++;
v=binomialReduce[G,v,monomialOrder,{Table[0,{n}],0},Internal`binomialTail];
dbgPrint[ToricGB`Debug`geometricBuchberger,G[[ToricGB`Private`basisList,i]]," reduces to ",v];
If[binomialZeroQ[v,projection],
deleteList[G,i,v];
,
v[[ToricGB`Private`binomialMarker]]=1;
changeList[G,i,v];
Sow[v];
];
]][[2,1]]][[1]];
clearList[G];
dbgAssert[dbgGroebnerBasisQ[H,monomialOrder,projection],geometricBuchberger,"G should be a Gröbner basis now. F=",F, ", monomialOrder=",monomialOrder,", grading=",gradingP,", projection=",projection,", criteria=",criteria];
dbgPrint[ToricGB`Debug`geometricBuchberger,"G=",H];
If[ToricGB`Debug`profile,
Print["sizeG=",sizeG,", cntInsert=",cntInsert,", cntQuery=",cntQuery,", maxExponent=",maxExponent];
cntInsert=0;cntQuery=0;
];
H
];

F = {{1, 2, -3}, {0, -2, 4}, {-2, 3, 1}}; dbgGroebnerBasisUnique[VectorToBinomial[binomialStrip[geometricBuchberger[binomialSupport[F], lexOrder[3],{{0,0,0},0},proj[3,3],3]], vars[3]],vars[3]]

dbgGroebnerBasisUnique[GroebnerBasis[VectorToBinomial[F, vars[3]], vars[3]],vars[3]]

Attention

The is not result equal since geometricBuchberger constructs a Gröbner basis of an ideal I < J < Sqrt[I].

F={{13,-64,0,29},{-5,24,1,-10}};
order={{8,(13/8),0,0},{1,0,0,0},{0,1,0,0},{0,0,1,0},{0,0,0,1}};
dbgGroebnerBasisUnique[VectorToBinomial[binomialStrip[geometricBuchberger[binomialSupport[F], order,{{0,0,0,0},0},proj[2,4],7]],vars[4]],vars[4]]

profile[geometricBuchberger[binomialSupport[{{1,0,38,9,-19,13},{0,1,366,78,-184,126},{0,0,1123,239,-565,386},{239,0,98,239,-21,19},{1003,1025,0,1025,263,141},{471,544,263,544,0,169}}],lexOrder[6],{{0,0,0,0,0,0},0},proj[6,6],3]]

VEC, bitmasks, geometricBuchberger, profiling, LatticeReduce, part intermediateReduce

VEC, bitmasks, geometricBuchberger, profiling, LatticeReduce, intermediateReduce

profile[geometricBuchberger[binomialSupport[{{1,1,-1,1,-2},{0,10,1,21,-46},{0,0,3,0,2},{2,2,1,2,-2},{42,22,7,0,42}}],lexOrder[5],{{0,0,0,0,0},0},proj[5,5],3]]

VEC, bitmasks, geometricBuchberger, profiling, LatticeReduce, part intermediateReduce

VEC, bitmasks, geometricBuchberger, profiling, LatticeReduce, intermediateReduce

projectionOrdering

Computes an ordering such that the Gröbner basis G with respect to this ordering and projection is also a Gröbner basis with respect to this ordering and Sort[Append[projection,i]].

SyntaxInformation[projectionOrdering]={"ArgumentsPattern"->{_,_,_,_}};
projectionOrdering[L_List,projection_List,i_Integer,monomialOrder_List]:=Module[{x,c,A,b,n=Length[L[[1]]],projectionComplement},
c=Table[1,{n}];
projectionComplement=Complement[Range[n],projection[[ToricGB`Private`vectorList]]];
A=FlattenAt[Append[L,Map[(UnitVector[n,#])&,projectionComplement]],Length[L]+1];
b=-Array[({L[[#,i]],0})&,Length[L]];
b=FlattenAt[Append[b,Table[{0,0},{Length[projectionComplement]}]],Length[b]+1];
x=LinearProgramming[c,A,b];
x=x*Apply[LCM,Denominator[x]];
Prepend[monomialOrder,x]
(*Prepend[Prepend[Array[-(UnitVector[n,#])&,n],Table[1,{n}]],x]*)
];

L={{1,-5,-3,3},{0,8,4,-5}};
order=projectionOrdering[L,positionsToBitMask[{1,2,3},4],4,{}]

L.order[[1]]

L={{1,6,-7,-4,2,-15},{0,7,2,-8,9,-19},{0,0,23,-10,19,-7}};
order=projectionOrdering[L,positionsToBitMask[{1,2,3,4,5},6],6,{}]

L.order[[1]]

HemmeckeToricGroebnerBasis

The integer matrix A defines a lattice L = { v integer vector: A.v=0 }. (see also LatticeNullSpace)
The corresponding monomial ideal is defined by I = < (x^v+) - (x^v-) : v in L >.
ToricGroebnerBasis computes a Gröbner basis G (given as list of polynomials) of I .

LatticeReduce is used as preprocessing such that the vectors have small norms. Then a projection is found such that the lift of this projection is unique and we easily can compute a generating set (aka Markov basis). For this projection there always exists a totally positive element u. Adding this to the basis G ensures that it is a generating set. Now we compute the unbounded variables and lift the generating set wrt. them. Finally we lift wrt. the bounded variables and compute a Gröbner basis.

SyntaxInformation[HemmeckeToricGroebnerBasis]={"ArgumentsPattern"->{_,_,_}};
HemmeckeToricGroebnerBasis[B_List,variables_List,monomialOrder_List]:=Module[{L,G,grading,projection,i,u,v,list,time,H,criteria,order,n=Length[B[[1]]],unbounded,A},
assert[Length[variables]==n,HemmeckeToricGroebnerBasis,"The number of variables provided is wrong."];
A=LatticeReduce[B];
L=LatticeByNullSpace[A];
If[L==={},Return[{}]];
L=LatticeReduce[L];
dbgPrint[ToricGB`Debug`HemmeckeToricGroebnerBasis,"HemmeckeToricGroebnerBasis: A=",A,", L=",L, "monomialOrder=",monomialOrder];
G=binomialSupport[L];
(* find projection and suitable generating set *)
projection=latticeZeroKernelProjection[L];
{u,list}=latticeMaximalPositiveElements[L,projection]; (* Could be done faster using Integer-Gauß since projection of L is full rank *)
assert[Apply[And,Thread[u[[ToricGB`Private`binomialExponent,projection[[ToricGB`Private`vectorList]]]]>0]],HemmeckeToricGroebnerBasis,"No totally positive element in projection of lattice found."];
dbgPrint[ToricGB`Debug`HemmeckeToricGroebnerBasis,"Initial projection=",projection,", G=",G,", u=",u];
G=Join[G,list];
(* find globally bounded / unbounded variables *)
{u,list}=latticeMaximalPositiveElements[L,proj[n,n]];
unbounded=Select[Complement[Table[i,{i,1,n}],projection[[ToricGB`Private`vectorList]]],(u[[ToricGB`Private`binomialExponent,#]]≤0)&];
dbgPrint[ToricGB`Debug`HemmeckeToricGroebnerBasis,"Unbounded variables=",unbounded,", u=",u];
(* lift generating set from projection to all unbounded variables *)
While[Length[unbounded]>0,
i=First[unbounded];unbounded=Rest[unbounded];
v=latticePositiveElement[L,projection,i];
dbgPrint[ToricGB`Debug`HemmeckeToricGroebnerBasis,"G=",G,", projection=", projection,", coordinate=",i,": latticePositiveElement=",u];
If[v=!={},
If[!MemberQ[binomialStrip[G],binomialStrip[v]],G=Append[G,v]],
grading=maximalGrading[L,projection];
criteria=ToricGB`Private`criterionSupport;
If[True||Count[grading[[ToricGB`Private`vectorList]],x_/;x>0]<=Length[variables]/2,criteria=BitOr[criteria,ToricGB`Private`criterionMinimalSet]];
If[grading[[ToricGB`Private`vectorSupport]]≠0,criteria=BitOr[criteria,ToricGB`Private`criterionTail]];
{time,G}=AbsoluteTiming[geometricBuchberger[G,projectionOrdering[L,projection,i,monomialOrder],grading,
projection,criteria]];
dbgPrint[ToricGB`Debug`HemmeckeToricGroebnerBasis,"geometricBuchberger: G=",G,", order=",projectionOrdering[L,projection,i,monomialOrder],", grading=",grading,", projection=",projection,", criteria=",criteria," : time = ",time];
dbgAssert[dbgGroebnerBasisQ[G,projectionOrdering[L,projection,i,monomialOrder],positionsToBitMask[Sort[Append[projection[[ToricGB`Private`vectorList]],i]],n]],HemmeckeToricGroebnerBasis,"G should be a Gröbner basis now."];
];
projection=positionsToBitMask[Sort[Append[projection[[ToricGB`Private`vectorList]],i]],n];
];
(* lift generating set to bounded variables *)
If[!binomialZeroQ[u,projection],
G=Join[G,list];
dbgPrint[ToricGB`Debug`HemmeckeToricGroebnerBasis,"Final lift with u=",u];
];
(* compute Gröbner basis *)
projection=proj[n,n];
grading=maximalGrading[L,projection];
criteria=ToricGB`Private`criterionSupport;
If[True||Count[grading[[ToricGB`Private`vectorList]],x_/;x>0]<=Length[variables]/2,criteria=BitOr[criteria,ToricGB`Private`criterionMinimalSet]];
If[grading[[ToricGB`Private`vectorSupport]]≠0,criteria=BitOr[criteria,ToricGB`Private`criterionTail]];
dbgAssert[dbgSaturatedQ[G],HemmeckeToricGroebnerBasis,"Ideal should be saturated by now. (A=",A,")"];
If[ToricGB`useToricGroebnerBasis,
{time,H}=AbsoluteTiming[GroebnerBasis`ToricGroebnerBasis[VectorToBinomial[binomialStrip[G],variables],variables]];
dbgPrint[ToricGB`Debug`HemmeckeToricGroebnerBasis,"GroebnerBasis`ToricGroebnerBasis: time=",time];
H,
dbgPrint[ToricGB`Debug`HemmeckeToricGroebnerBasis,"geometricBuchberger: G=",G,", order=",monomialOrder,", grading=", grading,", projection=",projection,", criteria=",criteria," : time = ",time];
{time,H}=AbsoluteTiming[geometricBuchberger[G,monomialOrder,grading,projection,criteria]];
dbgPrint[ToricGB`Debug`HemmeckeToricGroebnerBasis,"geometricBuchberger: G=",G,", order=",monomialOrder,", grading=", grading,", projection=",projection,", criteria=",criteria," : time = ",time];
VectorToBinomial[binomialStrip[H],variables]]
];

Example 1

B[1]={{2,1,1}};
dbgGroebnerBasisUnique[HemmeckeToricGroebnerBasis[B[1],vars[3],lexOrder[3]],vars[3]]

dbgGroebnerBasisUnique[GroebnerBasis`ToricGroebnerBasis[Append[VectorToBinomial[LatticeByNullSpace[B[1]],vars[3]],1-t*Apply[Times,vars[3]]],vars[3],t],vars[3]]

Example 2

B[2]={{-6,-2,9,2},{-1,3,-10,1}};
dbgGroebnerBasisUnique[HemmeckeToricGroebnerBasis[B[2],vars[4],lexOrder[4]],vars[4]]

dbgGroebnerBasisUnique[
GroebnerBasis`ToricGroebnerBasis[Append[VectorToBinomial[LatticeByNullSpace[B[2]],vars[4]],1-t*Apply[Times,vars[4]]],vars[4],t],vars[4]]

Example 3

B[3]={{-3,8,-5,-8}};
dbgGroebnerBasisUnique[HemmeckeToricGroebnerBasis[B[3],vars[4],lexOrder[4]],vars[4]]

dbgGroebnerBasisUnique[GroebnerBasis`ToricGroebnerBasis[Append[VectorToBinomial[LatticeByNullSpace[B[3]],vars[4]],1-t*Apply[Times,vars[4]]],vars[4],t],vars[4]]

Example 4

B[4]={{-3,-10,-8}};
dbgGroebnerBasisUnique[HemmeckeToricGroebnerBasis[B[4],vars[3],lexOrder[3]],vars[3]]

dbgGroebnerBasisUnique[GroebnerBasis`ToricGroebnerBasis[Append[VectorToBinomial[LatticeByNullSpace[B[4]],vars[3]],1-t*Apply[Times,vars[3]]],vars[3],t],vars[3]]

Example 5

B[5]={{-9,-1,3,-5,-5,10},{5,0,-2,3,5,1}};
dbgGroebnerBasisUnique[HemmeckeToricGroebnerBasis[B[5],vars[6],lexOrder[6]],vars[6]]

dbgGroebnerBasisUnique[GroebnerBasis`ToricGroebnerBasis[Append[VectorToBinomial[LatticeByNullSpace[B[5]],vars[6]],1-t*Apply[Times,vars[6]]],vars[6],t],vars[6]]

Example 6

B[6]={{2,-9,-4}};
dbgGroebnerBasisUnique[HemmeckeToricGroebnerBasis[B[6],vars[3],lexOrder[3]],vars[3]]

dbgGroebnerBasisUnique[GroebnerBasis`ToricGroebnerBasis[Append[VectorToBinomial[LatticeByNullSpace[B[6]],vars[3]],1-t*Apply[Times,vars[3]]],vars[3],t],vars[3]]

Example 7

B[7]={{-1,3,-5,-5,10}};
dbgGroebnerBasisUnique[HemmeckeToricGroebnerBasis[B[7],vars[5],lexOrder[5]],vars[5]]

dbgGroebnerBasisUnique[GroebnerBasis`ToricGroebnerBasis[Append[VectorToBinomial[LatticeByNullSpace[B[7]],vars[5]],1-t*Apply[Times,vars[5]]],vars[5],t],vars[5]]

Example 8

B[8]={{3,9,7,-4},{6,2,7,7}};
dbgGroebnerBasisUnique[HemmeckeToricGroebnerBasis[B[8],vars[4],lexOrder[4]],vars[4]]

dbgGroebnerBasisUnique[
GroebnerBasis`ToricGroebnerBasis[Append[VectorToBinomial[LatticeByNullSpace[B[8]],vars[4]],1-t*Apply[Times,vars[4]]],vars[4],t],vars[4]]

Example 9

B[9]={{-10,5,-9,-8,10,9},{-6,-2,-4,-3,4,2},{-3,3,-3,5,7,2}};
dbgGroebnerBasisUnique[HemmeckeToricGroebnerBasis[B[9],vars[6],lexOrder[6]],vars[6]]

dbgGroebnerBasisUnique[GroebnerBasis`ToricGroebnerBasis[Append[VectorToBinomial[LatticeByNullSpace[B[9]],vars[6]],1-t*Apply[Times,vars[6]]],vars[6],t],vars[6]]

Example 10

B[10]={{1,6,-5,-5},{0,-6,-10,-9}};
dbgGroebnerBasisUnique[HemmeckeToricGroebnerBasis[B[10],vars[4],lexOrder[4]],vars[4]]

dbgGroebnerBasisUnique[GroebnerBasis`ToricGroebnerBasis[Append[VectorToBinomial[LatticeByNullSpace[B[10]],vars[4]],1-t*Apply[Times,vars[4]]],vars[4],t],vars[4]]

Example 11

B[11]={{-10,9,9}};
dbgGroebnerBasisUnique[HemmeckeToricGroebnerBasis[B[11],vars[3],lexOrder[3]],vars[3]]

dbgGroebnerBasisUnique[GroebnerBasis`ToricGroebnerBasis[Append[VectorToBinomial[LatticeByNullSpace[B[11]],vars[3]],1-t*Apply[Times,vars[3]]],vars[3],t],vars[3]]

Example 12

B[12]={{-8,0,-7,10,-9,1},{9,2,3,-10,-1,-4},{9,0,-2,-9,-3,7}};
dbgGroebnerBasisUnique[HemmeckeToricGroebnerBasis[B[12],vars[6],lexOrder[6]],vars[6]]

dbgGroebnerBasisUnique[GroebnerBasis`ToricGroebnerBasis[Append[VectorToBinomial[LatticeByNullSpace[B[12]],vars[6]],1-t*Apply[Times,vars[6]]],vars[6],t],vars[6]]

Example 13

B[13]={{1,-7,-2,10,3},{-9,7,-4,10,6}};
dbgGroebnerBasisUnique[HemmeckeToricGroebnerBasis[B[13],vars[5],lexOrder[5]],vars[5]]

dbgGroebnerBasisUnique[GroebnerBasis`ToricGroebnerBasis[Append[VectorToBinomial[LatticeByNullSpace[B[13]],vars[5]],1-t*Apply[Times,vars[5]]],vars[5],t],vars[5]]

Example 14

B[14]={{-5,4,-1,-3,-3,10},{0,3,1,-7,7,-1},{7,-2,10,8,8,4}};
profile[dbgGroebnerBasisUnique[HemmeckeToricGroebnerBasis[B[14],vars[6],lexOrder[6]],vars[6]]]

KDL semibin, bitmasks, geometricBuchberger, profiling, LatticeReduce, intermediateReduce, FindInstance

VEC, bitmasks, geometricBuchberger, profiling, LatticeReduce, intermediateReduce, FindInstance

VEC, bitmasks, geometricBuchberger, profiling, LatticeReduce, part intermediateReduce, FindInstance

VEC, bitmasks, geometricBuchberger, profiling, LatticeReduce, part intermediateReduce, Minimize

kdl semibin, bitmasks, geometricBuchberger, profiling, LatticeReduce, part intermediateReduce, Minimize

VEC, bitmasks, geometricBuchberger, profiling, LatticeReduce, part intermediateReduce

kdl semibin, LatticeReduce+bounded variables

vec

kdl fullbin + vec

kdl nobin + vec (=kdl semibin + vec)

kdl fullbin

kdl semibin

kdl nobin

{31.8-10.3,132.0-106.3,51.0-22.6,124.4-105.5,43.8-26.3,43.2-26.3}

Example 15

B[15]={{1,-9,2,-6,-7,2},{-8,-1,7,7,-9,8},{10,-3,9,5,5,10}};
profile[dbgGroebnerBasisUnique[HemmeckeToricGroebnerBasis[B[15],vars[6],lexOrder[6]],vars[6]]]

KDL semibin, bitmasks, geometricBuchberger, profiling, LatticeReduce, intermediateReduce, FindInstance

VEC, bitmasks, geometricBuchberger, profiling, LatticeReduce, intermediateReduce, FindInstance

VEC, bitmasks, geometricBuchberger, profiling, LatticeReduce, part intermediateReduce, FindInstance

VEC, bitmasks, geometricBuchberger, profiling, LatticeReduce, part intermediateReduce, Minimize

kdl semibin, bitmasks, geometricBuchberger, profiling, LatticeReduce, part intermediateReduce, Minimize

VEC, bitmasks, geometricBuchberger, profiling, LatticeReduce, part intermediateReduce

kdl semibin, LatticeReduce+bounded variables

vec

kdl fullbin + vec

kdl nobin + vec (=kdl semibin + vec)

kdl fullbin

kdl semibin

kdl nobin

{55.3-7.1,64.3-2.2,61.8-4.0,59.9-2.2,58.0-4.0,53.1-4.0}

Example 16

B[16]={{-5,6,8,3,5,3},{7,9,-9,8,-3,10},{-2,4,-9,-3,-8,-5}};
profile[dbgGroebnerBasisUnique[HemmeckeToricGroebnerBasis[B[16],vars[6],lexOrder[6]],vars[6]]]

KDL semibin, bitmasks, geometricBuchberger, profiling, LatticeReduce, intermediateReduce, FindInstance

VEC, bitmasks, geometricBuchberger, profiling, LatticeReduce, intermediateReduce, FindInstance

VEC, bitmasks, geometricBuchberger, profiling, LatticeReduce, part intermediateReduce, FindInstance

VEC, bitmasks, geometricBuchberger, profiling, LatticeReduce, part intermediateReduce, Minimize

kdl semibin, bitmasks, geometricBuchberger, profiling, LatticeReduce, part intermediateReduce, Minimize

VEC, bitmasks, geometricBuchberger, profiling, LatticeReduce, part intermediateReduce

kdl semibin, LatticeReduce+bounded variables

vec

kdl fullbin + vec

kdl nobin + vec (=kdl semibin + vec)

kdl fullbin

kdl semibin

kdl nobin

{56.4-10.3,130.8-76.6,58.6-8.5,108.4-75.9,36.9-8.4,35.3-8.4}

Example 17

B[17]={{-2,-3,-7,6,-5,0},{-3,2,-4,9,-1,9}};
profile[dbgGroebnerBasisUnique[HemmeckeToricGroebnerBasis[B[17],vars[6],lexOrder[6]],vars[6]]]

KDL semibin, bitmasks, geometricBuchberger, profiling, LatticeReduce, intermediateReduce, FindInstance

VEC, bitmasks, geometricBuchberger, profiling, LatticeReduce, intermediateReduce, FindInstance

VEC, bitmasks, geometricBuchberger, profiling, LatticeReduce, part intermediateReduce, FindInstance

VEC, bitmasks, geometricBuchberger, profiling, LatticeReduce, part intermediateReduce, Minimize

kdl semibin, bitmasks, geometricBuchberger, profiling, LatticeReduce, part intermediateReduce, Minimize

VEC, bitmasks, geometricBuchberger, profiling, LatticeReduce, part intermediateReduce

kdl semibin, LatticeReduce+bounded variables

kdl fullbin + vec

kdl nobin + vec (=kdl semibin + vec)

kdl fullbin

kdl semibin

kdl nobin

{547.7-525.2,97.1-79.2,80.1-72.7,85.9-78.4,86.1-78.4}

Example 18

B[18]={{912489947,1256837063,1523174407,4535072249,6927997597,8532266279}};
order[18]={{912489947,1256837063,1523174407,4535072249,6927997597,8532266279},{-1,0,0,0,0,0},{0,-1,0,0,0,0},{0,0,-1,0,0,0},{0,0,0,-1,0,0},{0,0,0,0,-1,0}};
profile[HemmeckeToricGroebnerBasis[B[18],vars[6],order[18]]]

KDL semibin, bitmasks, geometricBuchberger, profiling, LatticeReduce, intermediateReduce, FindInstance

VEC, bitmasks, geometricBuchberger, profiling, LatticeReduce, intermediateReduce, FindInstance

VEC, bitmasks, geometricBuchberger, profiling, LatticeReduce, part intermediateReduce, FindInstance

VEC, bitmasks, geometricBuchberger, profiling, LatticeReduce, part intermediateReduce, Minimize

kdl semibin, bitmasks, geometricBuchberger, profiling, LatticeReduce, part intermediateReduce, Minimize