#//
#//  GAP functions to deal with toric varieties 
#//             and cones in Q^n
#//
#//
#//      thanks to A Hulpke for help
#//
#//wdj, 7-12-2002
#
# bug in in_dual_cone corrected 8-12-2002 
# (thanks to B Siebert for pointing it out)
#
#//////////////////////////////////////////////////////////

#RequirePackage("guava");
#Read("/home/wdj/gapfiles/toric.g"); 
#Read("c:/gap/gapfiles/toric.g"); 
#LogTo("c:/gap/logfiles/toric_gap1.log"); 

#a cone is a list of vectors spanning it 


flatten:=function(L)
#/*
#Takes a lists of lists and returns the elements
#*/
local L0,L1,i;
L0:=ShallowCopy(L[1]);     
if Size(L) = 1 then return L0; fi;  
for i in [2..Size(L)] do
  Append(L0, L[i]);
od;
return Set(L0);
end;
 
delta:=function(i,j) 
 if i = j then return 1; fi; 
 return 0; 
end; 
 
make_basis:=function(dim) 
local e,i,j;
e:=[]; 
for i in [1..dim] do 
 e[i]:=List(delta(i,j), j in [1..dim]); 
od;  
end; 
 
max_vectors:=function(L) 
#L is a list of vectors L:=[v1,...,vn] 
#returns max abs value of all coords
local n,max0,min0,i,max;
 n:=Size(L); 
 max:=0; 
 for i in [1..n] do 
  max0:=Maximum(L[i]); 
  min0:=Minimum(L[i]); 
  max:=Maximum([max,AbsoluteValue(max0),AbsoluteValue(min0)]); 
 od; 
 return max; 
end; 
 
inner_product:=function(v1,v2,A)
#v1 and v2 are n-vectors
#A is nxn pos def matrix
local i,v,inn;
v:=A*v2;
inn:=Sum([1..Size(v1)], i->v1[i]*v[i]);
return inn;
end;

in_dual_cone:=function(v,L) 
local numgens,dim,i,V;
 dim:=Size(v); 
 numgens:=Size(L); 
 V := Rationals^dim; 
 for i in [1..numgens] do 
  if inner_product(v,L[i],IdentityMat(dim)) < 0 then return false; fi; 
 od; 
 return true; 
end; 
 
#//  let L be a list of integral vectors generating a cone sigma 
#//  to find the generators of S_sigma,  
#//  let M be the max of the abs values of the 
#//  coords of the L[i]'s, for each vector v in 
#//  [1..M]^n, test if v in sigma^*. If so, 
#//  add to list of possible generators. 
#//  Once this for loop is finished, one can check 
#//  this list for redundant generators. The  
#//  trick below is to simply omit those elements 
#//  which are of the form d1+d2, where d1 and d2 
#//  are in the integral dual cone and "small" 
 
dual_semigp_gens:=function(L) 
#//L is a list of vectors L:=[v1,...,vn], 
#//v_i in Q^dim, generating the cone sigma 
#//returns semigp gens of S_sigma 
local d,i,dim,V,D,max,v,I,b,DpD,d1,d2,Dgens;
dim:=Size(L[1]); 
V := Rationals^dim; 
max:=max_vectors(L); 
D:=[]; 
I:= Cartesian(List([1..dim],i->[-max..max])); 
for v in I do 
 if in_dual_cone(v,L) then 
  Append(D,[v]); 
 fi; 
od; 
DpD:=[]; 
for d1 in D do 
 if inner_product(d1,d1,IdentityMat(dim)) <> 0 then 
 for d2 in D do 
  if inner_product(d2,d2,IdentityMat(dim)) <> 0 then 
     Append(DpD,[d1+d2]); 
  fi; 
 od; 
 fi; 
od; 
Dgens:=[]; 
for d in D do 
 if not(d in DpD) then 
    Append(Dgens,[d]); 
 fi; 
od; 
return Dgens; 
end; 

ideal_affine_toric_variety:=function(L)
local u,V,O,m,d,R,M,i,j,N,B,b0,J,terms,a,b,d0;
u:=dual_semigp_gens(L);
d:=Size(u[1]);
V:=Rationals^d;
O:=0*u[1];
RemoveSet(u,O);
m:=Size(u);
M:=NullMat(m,d);
for i in [1..m] do
 for j in [1..d] do
   M[i][j]:=u[i][j];
 od;
od;
B:=NullspaceMat(M);
b0:=Size(B);
R:=PolynomialRing(Rationals,d); 
if b0 = 0 then return R; fi; 
a:=[];
b:=[];
d0:=Size(B[1]);
for i in [1..b0] do a[i]:=[]; od;
for i in [1..d0] do b[i]:=[]; od;
for i in [1..b0] do
 for j in [1..d0] do
 if B[i][j] >= 0 then a[i][j]:=B[i][j]; else a[i][j]:=0; fi;
 if B[i][j] < 0 then b[i][j]:=-B[i][j]; else b[i][j]:=0; fi;
 od;
od;
terms:=[];
for j in [1..b0] do 
 terms[j]:= Product(List([1..d],i->X(Rationals,i)^a[j][i]))- 
Product(List([1..d0],i->X(Rationals,i)^b[j][i]));
od;
#Print("The ideal generators are :\n ",terms,"\n\n");
J:=Ideal(R,terms);
return J;
end;

embedding_affine_toric_variety:=function(L)
#/* 
#   L is a list gen a cone, as in dual_semigp_gens
#   returns the toroidal embedding of X=affine_toric_variety(L)
#*/
local map0,u,O,i,P1,P2,e,dim,V,R1,R2,hom;
u:=dual_semigp_gens(L);
dim:=Size(u[1]);
V:=Rationals^dim;
O:=Zero(V);
u:=Difference(u,[O]);
map0:=[];
for e in u do
  P1:=Filtered([1..dim],i->e[i] < 0);
  P2:=Filtered([1..dim],i->e[i] >= 0);
  if Size(P1) > 0 then 
    Add(map0, Product(P2,i->X(Rationals,i)^(e[i]))/Product(P1,i->X(Rationals,i+dim)^(-e[i])));
  fi;
  if Size(P1) = 0 then 
    Add(map0, Product(P2,i->X(Rationals,i)^(e[i])));
  fi;
od;
#R1:=PolynomialRing(Rationals,dim);
#R2:=PolynomialRing(Rationals,Size(u));
#hom:=PolynomialRingHomomorphismByVariables(R1,R2,IndeterminantsOfPolynomialRing(R1),map0);
#return [map0,hom];
return map0;
end;

#desing_affine_toric_variety:=function(L)
## 
##   L is a list gen a cone, as in dual_semigp_gens
##   returns the rational map of schemes
##   which resolves the singularites of X=affine_toric_variety(L)
##
#u:=dual_semigp_gens(L);
#V:=Parent(u[1]);
#O:=Rep(V);
#Exclude(~u,O);
#map:=<>;
#d:=Dimension(V);
#A:=create_affine_space(d);
#B:=AffineSpace(Rationals(),#u);
#Z:=Integers();
#for e in u do
#  P1:=[i:i in [1..d]|e[i] lt 0];
#  P2:=[i:i in [1..d]|e[i] ge 0];
#  if #P1 gt 0 and #P2 gt 0 then 
#    map:=Append(map, &*[A.i^(Z!e[i]):i in [1..d]|e[i] ge 0]/&*[A.i^(-Z!e[i]):i in [1..d]|e[i] lt 0]);
#  end if;
#  if #P1 eq 0 then 
#    map:=Append(map, &*[A.i^(Z!e[i]):i in [1..d]|e[i] ge 0]);
#  end if;
#  if #P2 eq 0 then 
#    map:=Append(map, 1/&*[A.i^(-Z!e[i]):i in [1..d]|e[i] lt 0]);
#  end if;
#end for;
#return map<A->B|[x: x in map]>; //works in magma 2.8, not in 2.7
#end function;

toric_points:=function(n,F)
#//returns the points in (F^*)^n
local T,L,i,x0;
T:=[];
for x0 in AsList(F ) do
 if x0<>Zero(F) then Add(T,x0); fi;
od;
L:=Cartesian(List([1..n],i->T));
return L;
end;

toric_code_old:=function(L,F)
##requires guava
#/* This function returns the same toric code as in
#J. P. Hansen, "Toric surfaces and error-correcting codes",
#except that the number of generators is not all the
#integral vectors in the polytope but rather the 
#semigroup generators of S_sigma.
#L is a list of integral vectors generating the cone sigma
#F is the finite field (**not** char F = 2)
#*/
local u,gens,d,V,n,i,Z,v,B0,B,C,C1,t,e,tmpvar;
u:=dual_semigp_gens(L);
gens:=[];
d:=Size(L[1]);
n:=Size(toric_points(d,F));
V:=F^n;         # VectorSpace(F,n);
Z:=Integers;
for v in u do
 Add(gens,v);
od;
B0:=[];
for e in gens do
 tmpvar:=List(toric_points(d,F),t->Product([1..d],i->t[i]^e[i]));
 Add(B0,tmpvar);
od;
B:=AsList(B0);
C:=GeneratorMatCode(One(F)*B,F); ##uses guava command
return C;
end;

toric_code:=function(L,F)
##requires guava
##
#/* This function returns the same toric code as in
#J. P. Hansen, "Toric surfaces and error-correcting codes",
#except that the number of generators is not all the
#integral vectors in the polytope but rather the 
#semigroup generators of S_sigma.
#L is a list of integral vectors (in Hansen's case,
#L is the list of integral vectors in a polytope)
#F is the finite field (**not** char F = 2)
#*/
local u,gens,d,V,n,i,Z,v,B0,B,C,C1,t,e,tmpvar;
gens:=L;
#Print(gens,"\n\n");
d:=Size(gens[1]);
n:=Size(toric_points(d,F));
V:=F^n;         # VectorSpace(F,n);
Z:=Integers;
B0:=[];
e:=gens[1];
for e in gens do
 tmpvar:=List(toric_points(d,F),t->Product([1..d],i->t[i]^e[i]));
 Add(B0,tmpvar);
od;
B:=AsList(B0);
C:=GeneratorMatCode(One(F)*B,F); ##uses guava command
return C;
end;

toric_codewords:=function(L,F)
#//lists the codewords of toric_code
local C,codewords,c;
codewords:=[];
C:=toric_code(L,F);
for c in C do
 Append(codewords,[c]);
od;
return codewords;
end;

divisor_polytope:=function(D,Rays)
#/*
#Rays is the list of vectors in the rays for the fan Delta,
#D is the list of coeffs for the Weil divisor D
#returns the expressions which must be >0
#*/
local n,Inequalities,R,i,j;
 n:=Size(Rays[1]);
 Inequalities:=[];
 R:=PolynomialRing(Rationals,n);
 for i in [1..Size(Rays)] do
  Add(Inequalities,D[i]+Sum([1..n],j->X(Rationals,j)*Rays[i][j]));
 od;
 return Inequalities;
end;

#
polytope_lattice_points:=function(A,Perps)
#
#returns the list of points in P \cap L^\perp
#Perps=[v1,...,vk] is the list of "inward normal" 
#vectors perpendicular to the walls of the polytope P,
#A=[a1,...,ak] is the amount wall is shifted from the origin (ai>=0)
#Eg, x=0, x=a, y=0, y=b has Perps=[[1,0],[-1,0],[0,1],[0,-1]] 
#and A=[0,a,0,b]
#
local R,n,M,Pts,i,j,Inequalities,L0,L1,L2,v,ineq;
R:=ShallowCopy(Perps);
n:=Size(R[1]);
M:=Maximum(List([1..Size(A)],i->AbsoluteValue(A[i])));
L0:=List([1..n],i->[-M..M]);
Pts:=Cartesian(L0);
 L2:=[];
 for v in Pts do
 Inequalities:=[];
 for i in [1..Size(R)] do
  ineq:=A[i]+Sum([1..n],j->v[j]*R[i][j]);
  Add(Inequalities,ineq);
 od;
 if Minimum(Inequalities) >= 0 then Add(L2,v); fi;
 od;
return L2;
end;

#
divisor_polytope_lattice_points:=function(D,Cones,Rays_ordered)
#/*
#returns the list of points in P_D \cap L^\perp
#which paramterize the elements in L(D)
#Rays is the list of vectors in the rays for the fan Delta,
#Note the rays do not determine the cones in a fan.
#
#Rays_ordered is the *ordered* list of rays for the fan Cones
#*/
local R,n,M,Pts,i,j,Inequalities,L0,L1,L2,v,ineq,Rays;
Rays:=Rays_ordered;
R:=ShallowCopy(Rays);
#Print(R,"\n");
n:=Size(R[1]);
M:=Maximum(List([1..Size(D)],i->AbsoluteValue(D[i])));
L0:=List([1..n],i->[-M..M]);
Pts:=Cartesian(L0);
 L2:=[];
  for v in Pts do
 Inequalities:=List([1..Size(R)],i->D[i]+Sum([1..n],j->v[j]*R[i][j])); 
#Inequalities:=[];
 #for i in [1..Size(R)] do
 # ineq:=D[i]+Sum([1..n],j->v[j]*R[i][j]);
 # Add(Inequalities,ineq);
 #od;
#Print(v,"\n",Inequalities,"\n");
 if Minimum(Inequalities) >= 0 then Add(L2,v); fi;
 od;
return L2;
end;

riemann_roch:=function(D,Rays)
#/*computes a basis for the RR space of the
#divisor represented by D
#*/
local P,map0,d,R,e,i,P1,P2;
P:=divisor_polytope_lattice_points(D,Rays);
map0:=[];
d:=Size(Rays[1]);
R:=PolynomialRing(Rationals,d); 
for e in P do
  P1:=Filtered([1..d],i->e[i] < 0);
  P2:=Filtered([1..d],i->e[i] >= 0);
  if Size(P1) > 0 and Size(P2) > 0 then 
    Add(map0, Product(P2,i->X(Rationals,i)^(e[i]))/Product(P1,i->X(Rationals,i)^(-e[i])));
  fi;
  if Size(P1)= 0 then 
    Add(map0, Product(P2,i->X(Rationals,i)^(e[i])));
  fi;
  if Size(P2) = 0 then 
    Add(map0, 1/Product(P1,i->X(Rationals,i)^(-e[i])));
  fi;
od;
return map0;
end;

normal_to_hyperplane:=function(L)
#L is a list of n-1 vectors in Q^n
#returns a vector v s.t. <v,h>=0, for all h in span(L)
local normal;
normal:=NullspaceMat(TransposedMat(L));
return normal[1];
end;

faces:=function(Rays)
#returns all the normals to the faces (hyperplanes of the cone)
local i,normals,n0,N,m,R,dim_sets,dim,v,pos,neg,n1,v0,N0;
normals:=[];
R:=ShallowCopy(Rays);
dim:=Size(R[1]);
v:=Sum(R);
dim_sets:=Combinations(R,dim-1);
for m in dim_sets do
 n0:=normal_to_hyperplane(m);
 if inner_product(v,n0,IdentityMat(dim))>0 then Add(normals,n0);
   else
   Add(normals,-n0);
 fi;
od;
Sort(normals);
N0:=ShallowCopy(normals);
N:=ShallowCopy(normals);
for n1 in N do 
    if Minimum(List(Rays,v0->inner_product(v0,n1,IdentityMat(dim))))<0 then
        N0:=Difference(N0,[n1]); 
    fi;
od;        
return Set(N0);
end;

in_cone:=function(v,Rays)
#tests if v is in the cone represented by Rays
#(a cone is a list of vectors generating it)
#look at all dim-subsets of Rays and test if 
#there is a positive solns vector
local w,F,dim,R;
R:=ShallowCopy(Rays);
dim:=Size(R[1]);
F:=faces(R);
if Minimum(List(F,w->inner_product(v,w,IdentityMat(dim))))<=0 then 
 return false; 
fi;
return true;
end;

#the proc below takes a fan of cones and returns the 
#number of k-diml cones in the fan.
#each hyperplane of one cone should not form the face
#of any other cone in the fan.
#The idea is this: 
#the fan is represented as a set of max cones.
#for each max cone, look at the k-diml faces
#obtained by taking n choose k subsets of the
#rays describing the cone. **Certain** of these k-subsets yield the
#desired cones

number_of_cones_dim:=function(Cones,k)
local sigma,R,Rays,n,i,j,C,C0,dim,ell,k_cones,sigma0,N,v0,v;
if k=0 then return 1; fi;
C0:=ShallowCopy(Cones);
Rays:=flatten(C0);
R:=ShallowCopy(Rays);
if k=1 then return Size(R); fi;
k_cones:=[];
dim:=Size(R[1]);
if k>dim then Print("dimension error\n"); return 0; fi;
N:=Size(C0);
if k=dim then return N; fi;
for i in [1..N] do
 C:=ShallowCopy(Cones[i]);
 sigma:= Combinations(C,k);
 for ell in [1..Size(sigma)] do
     v0:=Sum(sigma[ell]);
     if not(in_cone(v0,C)) then  Add(k_cones,sigma[ell]); fi;
 od;
od;
return Size(Set(k_cones));
end;

subcones_of_fan:=function(Cones,k)
###returns k-diml subcones of the fan Cones
#only counts those subcones which are *not* interior
local sigma,R,Rays,n,x,i,j,C,C0,dim,ell,k_cones,sigma0,N,v0,v;
C0:=ShallowCopy(Cones);
if k=0 then return 0*C0[1][1]; fi;
Rays:=flatten(C0);
R:=ShallowCopy(Rays);
if k=1 then return List(R,x->[x]); fi;
k_cones:=[];
dim:=Size(R[1]);
if k>dim then Print("dimension error\n"); return 0; fi;
N:=Size(C0);
if k=dim then return C0; fi;
for i in [1..N] do
 C:=ShallowCopy(Cones[i]);
 sigma:= Combinations(C,k);
 for ell in [1..Size(sigma)] do
     v0:=Sum(sigma[ell]);
#Print("\n\n",ell," = ell  ,  v0 =  ",v0,"  C = ",C,"  in? = ",in_cone(v0,C));
  if not(in_cone(v0,C)) then  Add(k_cones,sigma[ell]); fi;
 od; ##excludes the "interior" cones
#Print("\n\n  i=",i,"  k_cones=  ",k_cones);
od;
return Set(k_cones);
end;

is_sublist:=function(L1,L2)
#returns true if L1 is a subset of L2
#L1,L2 are lists
local i,j,ans,x,y;
x:=Set(L1);
y:=Set(L2);
if IsSubset(y,x) then return true; fi;
return false;
end;

star:=function(cone,Cones)
    #computes the star of the cone in a fan
    #ie, the cones tau which have cone as a face
local T,tau,Taus,i,j,sigma,C0,R,Rays,dim;
sigma:=ShallowCopy(cone);
C0:=ShallowCopy(Cones);
Rays:=flatten(C0);
R:=ShallowCopy(Rays);
dim:=Size(R[1]);
Taus:=[];
for j in [1..dim] do
    tau:=subcones_of_fan(C0,j);
    Append(Taus,tau);
od;
T:=[];
for tau in Taus do
    if is_sublist(sigma,tau) then Append(T,[tau]); fi;
od;
return T;
end;

betti_number:=function(Cones,k0)
    #computes the k-th betti number of the toric variety 
    #X(Delta), where Delta is a fan defined by its maximal
    #cones, Cones
### X(Delta) **must be non-singular***
local i,j,n,a,C,k;
if IsOddInt(k0) then return 0; fi;
k:=Int(k0/2);
C:=ShallowCopy(Cones)[1][1];
n:=Size(C);
a:=Sum( [k..n],i->(-1)^(i-k)*Binomial(i,k)*number_of_cones_dim(Cones,n-i));
return a;
end;

cardinality_of_X:=function(Cones,q)
    #finds the number of elements of X(Delta) over the field
    #GF(q), where Delta is a fan defined by its maximal
    #cones, Cones
### X(Delta) **must be non-singular***
local i,j,n,a;
n:=Size(Cones[1][1]);
a:=Sum([0..n],i->(q-1)^(i)*number_of_cones_dim(Cones,n-i));
return a;
end;
 
euler_characteristic:=function(Cones)
    #finds the Euler char of X(Delta), where 
    #Delta is a fan defined by its maximal cones, Cones
    ### X(Delta) **must be non-singular***
local sigma,R,Rays,C0,dim;
C0:=ShallowCopy(Cones);
Rays:=flatten(C0);
R:=ShallowCopy(Rays);
dim:=Size(R[1]);
return number_of_cones_dim(Cones,dim);
end;

Print(" "); 
Print("toric.g for GAP 4.3, wdj, version 7-11-2002\n"); 
Print("available functions: dual_semigp_gens, cart_prod_lists,\n "); 
Print("in_dual_cone, max_vectors, inner_product, toric_points,\n ");
Print("ideal_affine_toric_variety, embedding_affine_toric_variety, \n ");
Print("toric_code, toric_codewords, divisor_polytope,\n");
Print("divisor_polytope_lattice_points, riemann_roch, flatten, \n");  
Print("faces, in_cone, normal_to_hyperplane, number_of_cones_dim , \n");  
Print("subcones_of_fan, star, betti_number, cardinality_of_X, \n");
Print("euler_characteristic, ... \n");  

###########################################################################
#
#                Examples 
#
#dual_semigp_gens([[1,0],[3,4]]);
#J:=ideal_affine_toric_variety([[1,0],[3,4]]);
#GeneratorsOfIdeal(J);
#toric_points(2,GF(5));
#C:=toric_code([[1,0],[3,4]],GF(3));
#Display(GeneratorMat(C));
#divisor_polytope([6,6,0],[[2,-1],[-1,2],[-1,-1]]);
#Cones1:=[[[2,-1],[-1,2]],[[-1,2],[-1,-1]],[[-1,-1],[2,-1]]];Div:=[6,6,0];
#Rays1:=[[2,-1],[-1,2],[-1,-1]];
#riemann_roch(Div,Rays1); 
#P_Div:=divisor_polytope_lattice_points(Div,Rays1);
#Cones2:=[ [ [2,0,0],[0,2,0],[0,0,2] ], [ [2,0,0],[0,2,0],[2,-2,1],[1,2,-2] ] ];
#faces(Cones2[1]); 
#faces(Cones2[2]); 
#number_of_cones_dim(Cones2,1);  
#number_of_cones_dim(Cones2,2);  
#number_of_cones_dim(Cones2,3); 
#star([[1,0]],Cones2);
#Cones3:=[ [ [2,0,0],[0,2,0],[0,0,2] ], [ [2,0,0],[0,2,0],[1,1,-2] ] ];
#star([[2,0,0],[0,2,0]],Cones3);
#subcones_of_fan(Cones2,1);
#subcones_of_fan(Cones3,2); 
#betti_number(Cones1,1);
#cardinality_of_X(Cones1,3);
#euler_characteristic(Cones1);
#
#
###########################################################################