--Algorithm:
--**********
-- Note: Several of these steps below will be done automatically. One now
-- simply has to do step 1 and call the function Inumbers, which will
-- do the rest of the procedure automatically. 


--1. Read in the polytope from the terminal as a list of vertices
--2. Save this in a variable P which is a matrix
--3. Take the transpose of P and call it N
--4. Generate the normal fan and save it in NF

--5. Next generate all the cones in the normal fan using maxCones(NF)
--   and save this list in TopDimCones

--6. Generate a list of all cones of dimension one less than the
--   the top dim. This can be done by going over the list of all 
--   top dim cones and generating from each such cone a list of
--   cones of one dim less than the top dim. Merge the list produced
--   from top dim cone or keep them separate associated to each top dim cone

--7. Generate a list of rays of the normalfan NF by calling rays(NF)

--8. Take a ray and a cone from the above list. We will compute their
--   intersection number here (using a well-known algorithm, for example, as in
--   Oda's book on Toric varieties.

--9. If the ray and the cone forms a top dim cone in the fan then this number is 1 
--   else
--       if dim of the cone is the same as the top dim, then this number is 0
--	 else
--	     there is a ray in the cone which is the same as the given ray.
--	     In this case we proceed as follows.
	     
--	     Find two top dim cones in the normal fan that both contain the
--	     given cone so that if the extra rays in the two top dim cones are
--	     X1 and X2, then X1+X2 belongs to the given cone. Then
--	     write X1+X2 +a1v1+...+ akvk=0 where vi are the rays in the 
--	     original cone with v1 being the given ray. 
	     
--	     Then the intersection number of the given ray and cone = a1.

--10. We have also modified the program to take in a fan as input and produce
--    the corresponding intersection numbers. For this, one needs to 
--    first create a normal toric variety (for eg: P3= projectiveSpace 3)
--    produces the 3-dim projective space. Now take its fan using F=fan(P3),
--    for example. Then call Inumbers2(F) to get a list of triples consisting
--    of rays of the fan, the cones of codim 1 and their interscetion number.
--    This list is saved in the global variable "L". The older function
--    that takes a list of vertices as input is now renamed "Inumbers1(F)".

--    We can use the above feature along with stellarSubdivsion to produce
--    interesting fans and easily compute all their intersection numbers.

--11. May be as a test if our functions are working correctly, we could use
--    our functions along with another one (which needs to be written) that
--    will check ampleness using the Toric Nakai criterion.
--    Then we can re-check ampleness using the built-in function isAmple
--    and see if they agree.



--  Disclaimer: Currently the program has been checked with a
--  number of examples in two and three dimensions, but not in any dimension higher.
--  I wrote these functions only very recently and so there is
--  likely to be plenty of room for improvements.

--List of global variables:
--**************************
-- NF: used by NormalFan this is the normal fan
-- TopdimCones:used by TDCones  this is the set of cones in the normal fan
-- TopConeslist:used by TopCones1  list of matrices whose columns are the rays
--     	    	      of the cones in TopdimCones
-- SmallConeslist: used by SmallCones  This is a list of matrices whose columns are
--     	    	 the rays of the cones of one dimension less than the top dimension
-- Rays:used by InteresectNo1   These are the rays of the normalfan
-- L: this is a list of triples, consisting of a ray of the normalfan, a
--    cone of dimension one less than the top dimension and their intersection number

---List of functions follows
--**************************

--  This function produces the normal fan from a list of vertices of a polytope
--  when these are entered as the rows of a matrix, with each vertex in one row
-- The global variable NF holds the normalfan
	 
   
   NormFan=(VL) ->(
	N:=transpose(VL);
	P:=convexHull(N);
        NF:=normalFan(P);
        R:=rays(NF);
	return(NF, R));

  

--This function takes a matrix as input and produces a list of submatrices obtained by
--removing one column at a time
--The global variable Ltemp contains the resulting list of submatrices
     RemoveCol=(M) ->(
	 nc:=rank(M);
	 Ltemp:={};
	 for i from 0 to (nc-1) do
	   (Mtemp:=submatrix'(M,, {i});  
	   Ltemp=append(Ltemp, Mtemp));
           return(Ltemp));

--This function produces a list of all the topdimensional cones from the fan and stores
-- these in the global variable TopdimCones
   TDCones0=(F) ->(TopdimCones0= maxCones(F));
  

--This functions takes a list of cones (L) and the rays of the fan (R) and returns a list of 
-- submatrices that correspond to the list of cones
   TDCones1 = (L, R) ->(
      n:=#L; print "n="; print n;
      L1={};
      R= submatrix'(R, ,{0})|R_{0};
      for i from 0 to (n-1) do
         (Ltemp:=L#i; print "Ltemp=";  print Ltemp;
         ctemp:=#Ltemp; print "ctemp ="; print ctemp;
         M:=R_{Ltemp#0}; print "M ="; print M;
         for j from 1 to (ctemp-1) do
            (Mtemp:=M|R_{Ltemp#j};
             M=Mtemp; print "Mtemp ="; print Mtemp);
         print "ctemp ="; print ctemp; print "M ="; print M;
         print "L1="; print L1;
         L1=append(L1, M); print "L1="; print L1);
      return(L1));

--This functions takes a list of cones (L) and the rays of the fan (R) and returns a list of 
-- submatrices that correspond to the list of cones
   TDCones2 = (L, R) ->(
      n:=#L; print "n="; print n;
      L1={};
  --  R= submatrix'(R, ,{0})|R_{0};
      for i from 0 to (n-1) do
         (Ltemp:=L#i; print "Ltemp=";  print Ltemp;
         ctemp:=#Ltemp; print "ctemp ="; print ctemp;
         M:=R_{Ltemp#0}; print "M ="; print M;
         for j from 1 to (ctemp-1) do
            (Mtemp:=M|R_{Ltemp#j};
             M=Mtemp; print "Mtemp ="; print Mtemp);
         print "ctemp ="; print ctemp; print "M ="; print M;
         print "L1="; print L1;
         L1=append(L1, M); print "L1="; print L1);
      return(L1));
     

     
    TDCones =(VL) ->(
      NFR:=NormFan(VL);
      NF= NFR#0;
      Rays= NFR#1;
      TDCones0(NF);
      TopdimCones=TDCones1(TopdimCones0, Rays));
   
--This function produces a list of all the cones of dimension one less than the top dimenion
--The global variable SmallConeslist will contain the list of theses cones stored as matrices
-- whose columns will be the rays of these cones  
   SmallCones=(T) ->(
	Conecount=#TopdimCones; print Conecount;
	SmallConeslist={};
        P:=posHull(TopdimCones#0);
	d:=dim(P); print d;
	for i from 0 to (Conecount-1) do 
	       (C=TopdimCones#i; print "C="; print C;
		Ltemp:=RemoveCol(C); print "Ltemp="; print Ltemp;
	        SmallConeslist= join(SmallConeslist,  Ltemp);  
		SmallConeslist=unique(SmallConeslist);
                print "SmallConeslist="; print SmallConeslist));

--This function produces a list of matrices whose columns are the rays of the top dim cones
  TopCones=(T) ->(
	local  R;
	Conecount=#TopdimCones; print Conecount;
	TopConeslist={};
        P:=posHull(TopdimCones#0);
	d:=dim(P);
        print d;
	for i from 0 to (Conecount-1) do 
	       (C=TopdimCones#i; R:= rays(C); print "R="; print R;
		      TopConeslist= append(TopConeslist,  R);   TopConeslist=unique(TopConeslist); 
		      print "TopConeslist="; print TopConeslist));

	      

--This function takes as input a matrix whose rows are the vertices of a polytope and calls
-- the above functions to create a normalfan, its rays which are stored in Rays,
-- and also produces the list of top dimensional cones and the list of cones of one dimension less
-- stored as matrices of rays
  IntersectNo1=(VL) ->(
	TDCones(VL);
	SmallCones(TopdimCones); print "SmallConeslist=";
	print SmallConeslist);


 IntersectNo2=(F)->(
      NF=F; Rays=rays(F);print "rays of fan="; print Rays;
      TDCones0(NF);
      TopdimCones= TDCones1(TopdimCones0, Rays);
      SmallCones(TopdimCones); print "SmalConeslist";
      print SmallConeslist);
      
      
      
     


--This function takes a matrix and produces a list whose entries are the columns 
  MakelistofColumns=(M) ->(
     c:=rank(M);
     Columnlist={};
     for i from 0 to (c-1) do
       Columnlist=append(Columnlist, M_{i});
     return(Columnlist);
     print "Columnlist=";
      print Columnlist);
       
       
       
       
       
--This function tests if the columns of the second matrix are among the columns of the first
-- and returns true or false accordingly    
  Compare=(M1b, M2s)->(
       Cl1:={};
       Cl2:={};
       MakelistofColumns(M1b);
       Cl1=Columnlist;
       MakelistofColumns(M2s);
       Cl2=Columnlist;
       isSubset(Cl2, Cl1));
          
        
       
	
--This function takes list of Topdimensional cones in TopConeslist and a matrix from the SmallConeslist
--and produces a list of matrices corresponding to topdimensional cones that contain the given cone
--of one dimension less	
 SelectallTopCones=(TL, M) -> (
       	Conecount=#TopdimCones; print Conecount;
        i:=0;
	MatchingL:={};
	while (i<=Conecount-1) do
	  (if (Compare(TL#i, M)) 
	      then (C1:=TL#i; i1=i; print "Cone1="; print C1; print "i1="; print i1;
		    MatchingL=append(MatchingL, C1); i=i+1)
	      else i=i+1);
        return(MatchingL));

--Given three matrices as input with the first two representing two topdimensioanl cones and
--the third a cone of one dimension less, it extracts two rays the first of which is called X1
-- and represents the new column in the first matrix and the second ray is called X2 and it represents
--the new column from the second matrix
	
  Extractrays=(C1, C2, Cs)->(
       i:=0;r=rank(C1);
       while ((i<=r-1) and not(Compare(submatrix'(C1,,{i}), Cs))) do
	    i=i+1; print "i="; print i;
       if (i==r) then print "Error" 
           else X1:=C1_{i}; print "X1="; print X1;
	        i=0;r=rank(C2);
                while ((i<=r-1) and not(Compare(submatrix'(C2,,{i}), Cs))) do
	        i=i+1; print "i="; print i;
                if (i==r) then print "Error" 
                 else X2:=C2_{i}; print "X2="; print X2;
                 return(X1, X2));
 
        
 

--This function takes as input a list of all topdimensional cones, i.e TopConeslist, and a 
--an entry of the SmallConeslist and produces two topdimensional cones that contain the
--smaller cone and so that the two extra rays X1 and X2 have the property that X1+X2 is in the
--linear span of the rays of the given entry of the SmallConeslist. The function returns X1 and X2.
  	 
 SelecttwoTopCones=(TL, M) ->(
       MatchingList:=SelectallTopCones(TL, M);
       i:=0;
       j:=1; print "MatchingL="; print MatchingList;
       (Y1,Y2):=Extractrays(MatchingList#i, MatchingList#j, M);
       r2:=rank(M | Y1 + Y2);
       r1:=rank(M);
       print ("r1="); print (r1);
       print ("r2="); print (r2);
       while ((r2>r1) and (#MatchingList >i)) do
         (while ((r2>r1) and (#MatchingList >j)) do
	    (j=j+1; print "j="; print j; print "i="; print i;
	     if ((i<#MatchingList) and (j< #MatchingList)) then
		  ((Y1, Y2)=Extractrays(MatchingList#i, MatchingList#j, M);
                   r2=rank(M | Y1 + Y2);
                   r1=rank(M)));
	  i=i+1);
       if ((i<#MatchingList) and (j< #MatchingList)) then
          (Cone1=MatchingList#i;
	   Cone2=MatchingList#j);
       return(Y1, Y2));
           
             
      

 --This function takes as input the matrix corresponding to an entry of SmallConeslist and a ray from Ray
 --and computes the corresponding interscetion number using the standard algorithm.
 --Though it seems to handle arbitrary dimensions, it has been checked with
-- examples only in 2 dimensions.	 
--So what needs to be fixed are these two functions, it seems
	
--Use Inumber2 since that seems to have correctly calculated the intersection numbers for
--Hirzebruch surface F_2.

    Inumber1=(Cn, Ry)->(
       BCn=posHull(Cn|Ry);
       bcn=dim(BCn);
       j:=0;
       P:=posHull(TopdimCones#0);
       d=dim(P);
       while((j<d-1) and (not(Ry ==Cn_{j}))) do
         j=j+1;
       if (j==(d-1)) then 
           ( Conecount=#TopdimCones;
             i:=0;
             Q:=posHull(TopdimCones_i);q:=dim(Q);
             while((i< Conecount) and not(BCn == Q) and (bcn==q)) do
               (i=i+1;  Q=posHull(TopdimCones_i); q=dim(Q));
	     if i<Conecount then IN=1 else IN=0; print "IN="; print IN)
	   else 
	     ((X1, X2):=SelecttwoTopCones(TopdimCones, Cn);
	      D1=det(transpose(Ry)*(X1+X2));
	      D= det(transpose(Ry)*(Ry));
	      IN=-D1/D);
        return (IN));
	

 Inumber2=(Cn, Ry)->(
       BCn=posHull(Cn|Ry);
       bcn=dim(BCn);
       j:=0; print "j="; print j;
       P:=posHull(TopdimCones#0);
       d=dim(P); print "dim(P)="; print d;
       while ((j<d-1) and not(Cn_{j} == Ry)) do
         ( j=j+1; print ("jin while loop="); print j);
       print "final value of j ="; print j;
       if (j==(d-1)) then 
           ( Conecount=#TopdimCones;
             i:=0;
             Q:=posHull(TopdimCones_i);
             q:=dim(Q);
             while((i< Conecount) and not(BCn == Q) and (bcn==q)) do
               (i=i+1; print ("i in while loop="); print (i); 
                if i<Conecount then Q=posHull(TopdimCones_i); q=dim(Q));
	     if i<Conecount then IN=1 else IN=0; print "IN="; print IN)
	   else 
             ( print "last part";
	      (X1, X2):=SelecttwoTopCones(TopdimCones, Cn);
              print "X1="; print X1; print "X2="; print X2;
              A:=promote(transpose(Cn)*Cn, QQ); 
              b:=promote(transpose(Cn)*(X1+X2), QQ);
              x:=inverse(A)*b;
	      IN=-det(x^{j}));
        return (IN));
	

 Inumber21=(Cn, Ry)->(
       BCn=posHull(Cn|Ry);
       bcn=dim(BCn);
       j:=0; print "j="; print j;
       P:=posHull(TopdimCones#0);
       d=dim(P); print "dim(P)="; print d;
       while ((j<d-1) and not(Cn_{j} == Ry)) do
         ( j=j+1; print ("jin while loop="); print j);
       print "final value of j ="; print j;
       if (j==(d-1)) then 
           ( Conecount=#TopdimCones;
             i:=0;
             Q:=posHull(TopdimCones_i);
             q:=dim(Q);
             while((i< Conecount) and not(BCn == Q)) do
               (i=i+1; print ("i in while loop="); print (i); 
                if i<Conecount then Q=posHull(TopdimCones_i); q=dim(Q));
	     if i<Conecount then IN=1 else IN=0; print "IN="; print IN)
	   else 
             ( print "last part";
	      (X1, X2):=SelecttwoTopCones(TopdimCones, Cn);
              print "X1="; print X1; print "X2="; print X2;
              A:=promote(transpose(Cn)*Cn, QQ); 
              b:=promote(transpose(Cn)*(X1+X2), QQ);
              x:=inverse(A)*b;
	      IN=-det(x^{j}));
        return (IN));




Inumber3=(Cn, Ry)->(
       j:=0;
       P:=posHull(TopdimCones#0);
       d=dim(P); print "dim(P)="; print d;
       while (not(Ry == Cn_{j}) and (j< d-1)) do
           (j=j+1; print ("jin while loop="); print j);
       print "final value of j ="; print j;
       if (j==(d-1)) then 
           ( BCn=posHull(Cn|Ry); bcn=dim(BCn);
             Conecount=#TopdimCones; print "Conecount="; print Conecount;
             i:=0;
             Q:=posHull(TopdimCones_i);
             q:=dim(Q); print "q="; print q;
             while((i< Conecount) and not(BCn == Q) and (bcn==q)) do
               (i=i+1; print ("i in while loop="); print (i); 
                if i<Conecount then Q=posHull(TopdimCones_i); q=dim(Q));
	     if i<Conecount then IN=1 else IN=0; print "IN="; print IN)
	   else 
              print "last part";
	     ((X1, X2):=SelecttwoTopCones(TopdimCones, Cn);
              print "X1="; print X1; print "X2="; print X2;
              A:=promote(transpose(Cn)*Cn, QQ); 
              b:=promote(transpose(Cn)*(X1+X2), QQ);
              x:=inverse(A)*b;
	      IN=-det(x^{j}));
             
        return (IN));
	
test =(Cn, Ry)->(
j:=0;
while ((j<1) and not(Cn_{j} == Ry)) do
  (j=j+1; print "j in loop"; print j);
print "final value of j ="; print j);

--This function takes as input the list of vertices of a polytope as the rows of a matrix and
--computes all the relevant intersection numbers. It calls on the above functions first to
--generate a list of cones of one dimension less than the top dimension, it generates
--the rays of the normal fan, and produces a list whose entriess are triples consisting
--of a matrix whose columns are the rays of a cone of one dim less than the top dimension,
--a ray of the normal fan and the corresponding intersection number.
-- Of the three functions below, Inumbers1 and Inumbers2 are the correct one.     
-- Inumbers2 takes a fan as input and Inumbers1 takes the vertices of a polytope instead.
     
      
	 
  Inumbers1=(VL)-> (
	 IntersectNo1(VL);
	 smallcount:=#SmallConeslist; print "smallcount="; print smallcount;
	 raycount:= numgens source(Rays) ; print "raycount="; print raycount;
         L={};
         for i from 0 to (raycount-1) do
	    (for j from 0 to  (smallcount-1) do
	       (INtemp:=Inumber2(SmallConeslist#j, Rays_{i});
		X:=(Rays_{i}, SmallConeslist#j,  INtemp);
		L=append(L, X); print "j="; print j;
		print "inside L"; print L;
		j=j+1);
	        print "i="; print i);
                return (L));

	 
  
  Inumbers2=(F)-> (
	 IntersectNo2(F);
	 smallcount:=#SmallConeslist; print "smallcount="; print smallcount;
	 raycount:= numgens source(Rays); print "raycount="; print raycount;
	 L={};
	 for i from 0 to (raycount-1) do
            (print "i in loop ="; print i;
	     for j from 0 to  (smallcount-1) do
                (print "j in loop ="; print j;
	        INtemp:=Inumber21(SmallConeslist_j, Rays_{i});
		X:=(Rays_{i}, SmallConeslist_j,  INtemp);
		L=append(L, X); print "j="; print j;
		print "inside L"; print L;
		j=j+1);
	        print "i="; print i);
         return (L));
	 	 

      
	Inumbers3=(F)-> (
      -- IntersectNo2(F);
	 smallcount:=#SmallConeslist; print "smallcount="; print smallcount;
	 raycount:= numgens source(Rays); print "raycount="; print raycount;
	 L={};
	 for i from 0 to (raycount-1) do
            (print "i in loop ="; print i;
	     for j from 0 to  (smallcount-1) do
                (print "j in loop ="; print j;
	        INtemp:=Inumber21(SmallConeslist_j, Rays_{i});
		X:=(Rays_{i}, SmallConeslist_j,  INtemp);
		L=append(L, X); print "j="; print j;
		print "inside L"; print L;
		j=j+1);
	        print "i="; print i;
                L=append(L, "new divisor"););
         return (L));

      Test1=(n,m)->(
          for j from 0 to n do
               (for i from 0 to m do 
                 (print "i="; print i);
               print  "j="; print j));




      Inumbers31=(F)-> (
      -- IntersectNo2(F);
	 smallcount:=#SmallConeslist; print "smallcount="; print smallcount;
	 raycount:= numgens source(Rays); print "raycount="; print raycount;
	 L={};
         for j from 0 to (smallcount-1) do
	    (for i from 0 to (raycount-1) do
              (print "i in loop ="; print i;
	       INtemp:=Inumber21(SmallConeslist_j, Rays_{i});
		X:=(SmallConeslist_j, Rays_{i}, INtemp);
		L=append(L, X); print "i="; print i;
		print "inside L"; print L;);
	        print "j="; print j;
                L=append(L, "new row"););
         return (L));
		 	 

     

       
 Inumberold=(Cn, Ry)->(
       r1:=rank(Cn); r2:= rank(Cn|Ry);
       BCn=posHull(Cn|Ry);
       if (r2>r1) then 
           ( Conecount=#TopdimCones;
             i:=0;
             while((i< Conecount) and not(BCn == TopdimCones#i)) do
               i=i+1;
	     if i<Conecount then IN=1 else IN=0; print "IN="; print IN)
	   else 
	     ((X1, X2):=SelecttwoTopCones(TopConeslist, Cn);
	      D1=det(transpose(Ry)*(X1+X2));
	      D= det(transpose(Ry)*(Ry));
	      IN=-D1/D);
        return (IN));
	
Inumbersold1=(VL)-> (
	 IntersectNo1(VL);
	 smallcount:=#SmallConeslist; print "smallcount="; print smallcount;
	 raycount:=#Rays; print "raycount="; print raycount;
	 L={};
	 for i from 0 to (raycount-1) do
	    (for j from 0 to  (smallcount-1) do
	       (INtemp:=Inumber1(SmallConeslist#j, Rays#i);
		X:=(Rays#i, SmallConeslist#j,  INtemp);
		L=append(L, X); print "j="; print j;
		print "inside L"; print L;
		j=j+1);
	        print "i="; print i));

---The following function produces the fan for the projective space of dim d
--blown-up at (d-1)-points. It starts with the fan for P^d, where the rays
--are the columns of the id-matrix and the column with all -1s. The last ray is denoted u_{d+1}. Now we 
--add rays u_{d+2}, ..., u_{d+d}, which are given by $u_{d+1+i} = the (d-i)-th column of the identity matrix
-- with a negative sign, for all i=1, ..., d-1. (Recall the columns of the id-matrix are indexed by the integers
--0,..., d-1.


pSpaceBlownup = (d)->(
     V := {toList(d:-1)} | entries id_(ZZ^d);
     Temp:=entries id_(ZZ^d);
     for i from 1 to (d-1) do
       V=append(V, -Temp#(d-i));  
     print ("V="); print (V);
     L1:={};
     for i from 1 to d do
       L1=append(L1, i);
     L2:={0};
     for i from 2 to d do
      L2= append(L2, i);
     F:= {L1, L2};
     for j from 1 to (d-1) do
        (temp:={d+j}; 
          for i from 1 to d-j do
            temp =append(temp, i);
           print("i="); print (i); print ("j="); print (j);
           print("oldtemp="); print (temp);
          if (d-j+2<d+1) then 
              (for i from (d-j+2) to d do
                 temp=append(temp, i));
           print("oldtemp1="); print (temp);
       F=append(F, temp);
       print ("initial F="); print (F));
     for j from 1 to (d-1) do
        ( print ("j0="); print (j);
          temp1:={0,d+j}; 
          print ("j="); print (j);
          print ("outer temp1="); print (temp1);
          for i from 1 to d do
            (temp1 =append(temp1, i);
             print ("i="); print (i));
          print("inner temp1="); print (temp1);
          for i from 2 to d do
            (temp:=temp1;
             print ("inner j="); print (j);
             temp= delete(temp#(d-j+2), temp);
             print("inner temp="); print (temp);
             temp=delete(temp#(i), temp);
             print("final temp="); print (temp);
             F=append(F, temp);
             print ("F="); print (F); print("inner i="); print (i)));
     X := normalToricVariety(V,F);
     return X);



pSpaceBlownup1 = (d)->(
     V := {toList(d:-1)} | entries id_(ZZ^d);
     Temp:=entries id_(ZZ^d);
     for i from 1 to (d-1) do
       V=append(V, -Temp#(d-i));  
     print ("V="); print (V);
     L1:={};
     for i from 1 to d do
       L1=append(L1, i);
     L2:={0};
     for i from 2 to d do
      L2= append(L2, i);
     F:= {L1, L2};
     for j from 1 to (d-1) do
        (temp:={d+j}; 
          for i from 1 to d-j do
            temp =append(temp, i);
           print("i="); print (i); print ("j="); print (j);
           print("oldtemp="); print (temp);
          if (d-j+2<d+1) then 
              (for i from (d-j+2) to d do
                 temp=append(temp, i));
           print("oldtemp1="); print (temp);
       F=append(F, temp);
       print ("initial F="); print (F));
     for j from 1 to (d-1) do
        ( print ("j0="); print (j);
          temp1:={0,d+j}; 
          print ("j="); print (j);
          print ("outer temp1="); print (temp1);
          for i from 1 to d do
            (temp1 =append(temp1, i);
             print ("i="); print (i));
          print("inner temp1="); print (temp1);
          for i from 1 to d do
            (temp:=temp1;
             temp=delete(temp#(i+1), temp);
             print("inner temp="); print (temp);
             print ("inner j="); print (j);
             temp= delete(temp#(d-j+1), temp);
             print("final temp="); print (temp);
             F=append(F, temp);
             print ("F="); print (F); print("inner i="); print (i)));
     X := normalToricVariety(V,F);
     return X);

---Notes on how to use the above functions-----
--1. Start M2
--2. loadPackage("NormalToricVarieties')
--3. load("Documents/M2/intersection_numbers9.m2")

--4. Next read in the various toric varieties, for e.g. X= projectiveSpace(d), hirzebruchSurface(r) etc.
--5. Next generate the fan by saying NFX= fan(X), Rays by rays(NFX).
--6. generate the maxCones by maxCones(NFX)
--7. It may be safer to generate the topdimensional cones and rays directly rather than using IntersectNo2(NFX)
--8. The following are key global variables used: TopdimCones, SmallConeslist, Rays

--9. One can try generating the TopdimCones by calling TDCones0(NFX) and then TDCones1(TopdimCones0, Rays)
--10. It seems the built in functions in Macaulay2 have some bugs at least for blow-up varieties 
--as we saw
--    in the case of Blow-ups of P2 at points.

--11. To make sure everythig is correct, observe that the rays are labelled counterclock-wise 
-- starting at  the x-axis, which will be a ray in all our examples.
--    One has to make sure the columns of Rays are listed in the correct order corresponding to ----the columns in TopdimCones.  

--   If there are errors here, manually correct the global variables TopdimCones, and Rays and call
--   SmallCones(TopdimCones) to generate SmallConeslist.

--12. Finally call Inumbers3(NFX) (in the case of corrections as above) or Inumbers2(NFX) (as in ---the case of Hirzebrucsurfaces.


--13. Regarding weights of the divisors.
---  See Cox, Little, Schenck p. 171, Examples.  They show that one writes
---  each divisor in terms of the basis vectors for the characters.
---  i.e. in the above example, there are 5 columns each corresponding to a toric divisor
     
--14. Then express each of these in terms of the standard basis for Z^3, since Z^= M.
      
--15. Now the weights of the divisors will give a 5x2 matrix, since CH^1(BP3) =Z \oplus Z.
--    The i-th row will be the weight of the i-th toric divisor.

--16. Now the exactness of the sequence 0 -->M --> \oplus _{\rho} Z\rho --> CH^1(BP3)-->0
--    is equivalent to condition that the product of the above two matrices is zero.
--    So the only condition is to chose te weights of the toric divisors such that this condition
--    is satisfied. 
 
--    In fact if we put the weights given in the "NormalToricVariety" package it does work out.

--The rest is copied from Session-March26.txt

P3=projectiveSpace(3);

NFP3=fan(P3);

maxCones(NFP3);

BP3=blowup({0,1,3}, P3);

B1P3= blowup({0, 1, 2}, P3);
NFB1P3=fan(B1P3);
MCB1P3=maxCones(NFB1P3);
RB1P3=rays(NFB1P3);

B2P3=blowup({0,1, 3}, B1P3);
NFB2P3 =fan(B2P3);
MCB2P3=maxCones(NFB2P3);
RB2P3=rays(NFB2P3);

NFBP3=fan(BP3);

M=maxCones(NFBP3);

BBP3 =blowup({0,2,3}, BP3);
NFBBP3=fan(BBP3);

R=rays(NFBBP3);

M2=maxCones(NFBBP3);

R=submatrix'(R,, {0})|R_{0};R=submatrix'(R,, {3})|R_{3};R=submatrix'(R,, {3})|R_{3};
R=submatrix'(R,, {4})|R_{4};Rays=R;

A0= matrix{{1,0, 0}, {0, 1, 0}, {0, 0, 1}};A1= matrix{{1,0, 0}, {0, 1, 0}, {0, 0, -1}};
A2= matrix{{1,0, 0}, {0, 0, -1}, {0, 1, 0}};A3= matrix{{1,-1, 0}, {0, -1, 0}, {0, -1, -1}};
A4= matrix{{1,-1, 0}, {0, -1, -1}, {0, -1, 0}};A5= matrix{{0, 0, -1}, {1, 0, -1}, {0, 1, -1}};
A6= matrix{{0, -1, 0}, {1, -1, 0}, {0, -1, -1}};A7= matrix{{0, -1, 0}, {0, -1, -1}, {1, -1, 0}};
B0=A0; B1=A1; B2=matrix{{1,0,-1},{0,0, -1},{0, 1, -1}}; B3=A3; B4=A5;
B5=matrix{{0, -1, -1}, {1, -1, 0}, {0, -1, 0}}; B6=matrix{{0, 0, -1},{1,0, 0},{0, -1, 0}};
B7=matrix{{-1, 0, -1}, {-1, 0, 0}, {-1, -1, 0}}

D0=matrix{{-1, 1, 0},{-1, 0, 1}, {-1, 0, 0}};D1=matrix{{-1, 1, 0},{-1, 0, -1}, {-1, 0, 0}};
D2=matrix{{-1, 0, 0}, {-1, 1, 0}, {-1, 0, 1}}; D3=matrix{{-1, 0, 0}, {-1, 0, -1}, {-1, 1, 0}};
D4=matrix{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};D5=matrix{{1, -1, 0}, {0, -1, -1},{0, -1, 0}};






TopdimCones ={A0, A1, A2, A3, A4, A5, A6, A7};