
ZNullSpace[mat_] := (* checked *)
  (* solve A*x == 0 for x in ZZ^r, return a list of generators of the solution lattice *)
  Extract[First[#], Position[Last[#], {0..}]]&[HermiteDecomposition[Transpose[mat]]]

LatticeIntersection[lattice1_, lattice2_] := (* checked *)
  Transpose[lattice1].#& /@ (Take[#, Length[lattice1]]& /@ ZNullSpace[Transpose[Join[lattice1, lattice2]]])

LatticeMemberQ[lattice_, p_] := Module[ {x},
  Length[FindInstance[Thread[p - Array[x, Length[lattice]].lattice == 0], Array[x, Length[lattice]], Integers]] > 0
];

ConeMemberQ[cone_, p_] := Module[ {x},
  CylindricalDecomposition[Join[Thread[Array[x, Length[cone]] >= 0], Thread[-p +
  Array[x, Length[cone]].cone == 0]], Array[x, Length[cone]]]=!=False
];

LoneleyQ[cone_, lattice_, p_] :=
  (* true iff p in cone and there does not exist x in lattice \ {0} such that p + x in cone *)
  Module[ {x, y},
    ConeMemberQ[cone, p] &&
    Length[FindInstance[Join[
      Thread[Array[y, Length[cone]] >= 0],
      Thread[p + Array[x, Length[lattice]].lattice == Array[y, Length[cone]].cone]
    ], Join[Array[x, Length[lattice]], Array[y, Length[cone]]], Integers, 2]] =!= 2
  ];

NonnegativeCone[mat_] :=
  (* solve A*x >= 0 for x in NN^r, return a list of generators of the solution cone *)
  (* only works when the resulting cone is line free, i.e., when ker A = {0} *)
  Module[ {x,red,m,i,j,u},
    red = Reduce[Thread[(mat.Array[x,Length[First[mat]]])>=0], Array[x,Length[First[mat]]], Integers];
    If[ And@@Table[!FreeQ[red, x[i]==0],{i,1,Length[First[mat]]}], Return[ {} ]];
    m = Max[Cases[red, C[j_] -> j, Infinity]];
    If[ Head[red] =!= And || First[red] =!= Element[Alternatives@@Array[C, m], Integers],
      Throw["encountered unexpected output format of Reduce (1)."];
    ];
    red = Rest[red];
    Do[ If[ red[[i]] =!= (C[i] >= 0), Throw["encountered unexpected output format of Reduce (2)."] ], {i, 1, m} ];
    red = Drop[red, m];
    red = Table[Cases[red, x[i] == u_ -> u], {i, 1, Length[First[mat]]}];
    If[ Union[Length /@ red] =!= {1}, Throw["encountered unexpected output format of Reduce (3)."] ];
    red = First /@ red;
    Table[red /. C[i] -> 1 /. C[_] -> 0, {i, 1, m}]
  ];

ConeIntersection[cone1_, cone2_] :=
  DeleteCases[Union[(Transpose[cone1].Take[#, Length[cone1]])& /@
    NonnegativeCone[Join[Transpose[Join[cone1, -cone2]], -Transpose[Join[cone1, -cone2]], IdentityMatrix[Length[cone1]+Length[cone2]]]]], {0..}]

CountLonelyPoints[cone_, lattice_, bound_:Infinity] :=
  Module[ {b, todo, bounds, n, next, e},
    (* returns min(n, bound) where n is the number of lonely points of the 
       cone wrt the lattice *)
    If[ ConeIntersection[cone, Join[lattice, -lattice]] =!= {}, Return[0] ];
    b = False;
    Do[
      If[ConeIntersection[cone, Join[lattice, -lattice, {c}]] === {c}, b = True; Break[]];
    , {c, cone}];
    If[ b, Return[Min[bound, Infinity]] ];
    (* finitely many; fglm style *)
    todo = IdentityMatrix[Length[cone]];
    bounds = {};
    n = 1; (* origin is lonely *)
    While[ Length[todo] > 0,
      next = First[todo]; todo = Rest[todo];
      If[ LoneleyQ[cone, lattice, next.cone],
        Print[next.cone];
        n += 1; (* new lonely point discovered *)
        If[ n >= bound, Break[] ];
        Do[ (* put its neighbors to the todo list *)
          b = True;
          Do[
            If[ ConeMemberQ[cone, (next + e - q).cone], b = False; Break[] ]
          , {q, bounds} ];
          If[ b, AppendTo[todo, next + e]]
        , {e, IdentityMatrix[Length[cone]]}];
      , (* not a lonely point, register it as a stopper *)
        AppendTo[bounds, next];
      ];
      todo = Union[todo];
    ];
    Return[n];
  ];

CountLonelyPoints[lattice_, bound_:Infinity] :=
  Module[ {d, n, basis},
    If[ Length[lattice] === 0, Return[Min[bound, Infinity]] ];
    d = Length[First[lattice]];
    n = 0; basis = Append[IdentityMatrix[d], 0 Range[d]];
    Do[
      n += CountLonelyPoints[DeleteCases[(#-b)& /@ basis, {0..}], lattice, bound - n];
      Print[DeleteCases[(#-b)& /@ basis, {0..}]];
      If[ n >= bound, Break[] ];
    , {b, basis}];
    Return[Min[n, bound]];
  ];

