--6/30/05
--author:Jessica Sidman
--Modification and annotation of code by Hosten and Smith from pg. 87 of "Computations in Algebraic Geometry with Macaulay 2", David Eisenbud, Daniel Grayson, Micahel Stillman, Bernd Sturmfels, eds., Springer-Verlag, 2002.
--toricIdeal has input a matrix A whose columns consist of 1 followed by lattice points.


toBinomial = (b,S) -> ( --S is a ring and b is an integer vector of length = number of variables in S
     	  pos := 1_S;  --We want to form the binomial x^(b+) - x^(b-). 
	               --pos will be x^(b+) and neg will be x^(b-).  We start by setting each monomial = 1.
	  neg := 1_S;
	  scan(#b, i-> if b_i > 0 then pos = pos*S_i^(b_i)     
	       else if b_i < 0 then neg = neg*(S_i)^(-b_i));
	       pos-neg); --Output the binomial.
	  
toricIdeal = (A) -> (
     m:= (rank source A)-1;  --This tells us how many variables the ring of the toric ideal should have.
     S := (ZZ/32003)(monoid [x_0..x_m, MonomialSize => 16]);--Declare the ambient ring of the toric ideal.
     B := transpose matrix syz A; --Find the generators for the kernel of A.
     J := ideal apply(entries B, b-> toBinomial(b,S)); --Create the ideal whose generators are the binomials 
     	       	    	      	   	     	       --corresponding to the generators of ker(A).
     scan(gens S, r-> J = saturate(J,r));     	       --Saturate the binomial ideal.

     J); --Output the ideal.