i:=3;

##############################################################
##############################################################
##  test examples:

if i=1 then
  g:=Group([[-1,0],[0,-1]],[[0,1],[1,0]]);
  gamma:=[1,0];
fi;

if i=2 then
  g:=Group([[-1,1],[-1,0]],[[0,1],[1,0]]);
  gamma:=[-1,1];
fi;

if i=3 then
  g:=Group(
   [[1,0,-1],[0,1,-1],[0,0,-1]],
   [[0,0,1],[1,0,0],[0,1,0]],
   [[0,-1,1],[0,-1,0],[1,-1,0]],
   [[0,1,-1],[1,0,-1],[0,0,-1]]
   );
  gamma:=[-1,-2,1];
fi;

###############################################################
###############################################################
##  The program:

Read("multinv.g");

el:=GroupOps.Elements(g);;
dim:=Length(el[1]);;
I:=IdentityMat(dim);;

##############################################################
# the reflections in  G :

refl:=Filtered(el, i->RankMat(i-I)=1);

###########################################################
# the root system Phi (actually only half of it as we don't
# include the negatives here):

Phi:=List(refl, i->SolveHomogenousEquations(TransposedMat(i)+I)[1]);

#############################################################
# the form :

form:=Sum(el, i->i*TransposedMat(i));

###############################################################
# selecting a base  Delta  of the root system  Phi :
# For now, there is a handpicked regular vector  gamma for each
# test example (see above).

Phiplus:=[];
for x in Phi do
  prod:=gamma*form*x;
  if prod>0 then Add(Phiplus,x); fi;
  if prod<0 then Add(Phiplus,-x); fi;
  if prod=0 then Print("BAD VECTOR 'gamma' !!!!!!!","\n");fi;
od;

Sums:=List(Cartesian(Phiplus,Phiplus), a -> a[1]+a[2]);
Delta:=Filtered(Phiplus, x -> not x in Sums);

#############################################################
# the fundamental weights :

wts:=List(Delta, y->List([1..dim], j->2*(I[j]*form*y)/(y*form*y)));
weights:=TransposedMat(wts^(-1));

###############################################################
# Finding monoid generators :

stretchers:=List(weights, v->Lcm(List(v, j -> Denominator(j))));

coeff:=Cartesian(List([1..dim], i->[0..stretchers[i]]));
coeff:=Sublist(coeff, [2..Length(coeff)]);

gens:=List(coeff, c-> Sum([1..dim], i -> c[i]*weights[i]));
gens:=Filtered(gens, w->ForAll(w,IsInt));

sums:=List(Cartesian(gens,gens), i -> i[1]+i[2]);
generators:=Filtered(gens, w -> not w in sums);
coefficients:=generators*TransposedMat(wts);

#########################################################
# orbits of fundamental weights:

orbits:=List(weights, i->Set(List(el,j->i*j)));

##########################################################
# Some output :

Print("\n","\n","The fundamental weights:", weights, "\n");
Print("\n","\n","The monoid generators:", generators,"\n");
Print("\n","\n","The coefficients (in terms of the fundamental weights):",
  coefficients,"\n");
Print("\n","\n","The orbits of the fundamental weights:", orbits,"\n");