/* Instructions:

1. The command Sol2(n,p,P) returns the Fq-rank of Bn(Z/n1+..+Z/nk) where q is the first prime >=P and p is input as a sequence [n1, n2,...,nk]

2. The command Check(a,m,P) returns the image of the symbol a under quotient map to P\otimes Bm(Z/n1+.+Z/nk) in field P where 
a is input as a symbol [ [a_1^1,...a_1^k], ... ,[a_m^1,..,a_m^k] ];
m is input as a sequence [n1,..,nk];
P is input as a string "p" or "Q": when p stands for Fp and "Q" stands for rational field.

for example:  Check([[6],[9],[12]],[16],"Q");  will return (   0 -1/2  1/2), which is the image of (6,9,12) in Q\otimes B3(Z/16);
  Check([[10],[12],[13]],[16],"Q"); will return (   0 -1/2  1/2), which is the image of (10,12,13) in Q\otimes B3(Z/16);
  This implies (6,9,12)=(19,12,13) in Q\otimes B3(Z/16);
  Check([[1,0],[0,1]],[2,2],"2"); will return (1,1), which is the image of ((1,0),(0,1)) in F2\otimes B2(Z/6+Z/18).

*** Note: The main purpose of this program is to check whether a symbol, or a linear combination of symbols are killed in the group of equivariant type.
if the Check command returns zero, it means the input is killed in that group. So to simplified the computation, we are dropping the generation condition
in construction of Bn. So the image and the rank above is not the real one. However since the symbols satisfying the generation conditions form a subgroup.
Dropping the generation condition will not disturb the result whether a symbol actually lies in the relation space. 

*/





function Index(a,b)
if #a eq 1 then
return a[1];
else
return Index(Remove(a,#a), Remove(b,#b))*b[#b]+a[#a];
end if;
end function;





function Position(a,p)
L:=#a;
a:=[Index(a[i],p): i in [1..L]];
p:=&*[p[i]: i in [1..#p]];
asorted:=Sort(a);
c:=Binomial(p-1+L,L);
index:=c-(&+[Binomial(p-2+L-asorted[i]-i+1,L-i+1):i in [1..L]]);
return index;
end function;




function FullList(n)
if #n eq 1 then
return {@ [i]:i in [0..n[1]-1] @};
else
tempresult:=FullList(Remove(n,#n));
result:=&join [&join[{@tempresult[i] cat [j]@}: j in [0..n[#n]-1]]:i in [1..#tempresult]];
return result;
end if;
end function;




function MinusOne(a,p);
i:=#a;
while a[i] eq 0 do
a[i]:=p[i]-1;
i:=i-1;
end while;

a[i]:=[a[i]-1];
return a;
end function;



function Eshort(n,grp)
k:=2;//Subcases k=2 generate all relations
p:=&*[grp[i]:i in [1..#grp]];
nvars:=Binomial(p-1+n,n);
neqs:=1+Binomial(p-1+k,k)*Binomial(p-1+n-k,n-k);

//Build the relation matrix
M:=SparseMatrix(RationalField(), neqs, nvars);

currentrow:=1;
M[currentrow,1]:=1;//Put (0..0)=0

achoice:=IndexedSetToSet(FullList(grp)) join {[grp[i]-1:i in [1..#grp-1]] cat [grp[#grp]]};

for aset in Subsets(achoice,k) do
alist:=Sort(SetToSequence(aset));
a:=[alist[1],MinusOne(alist[2],grp)];

bchoice:=IndexedSetToSet(FullList(grp)) join &join [{[grp[i]-1:i in [1..#grp-1]] cat [grp[#grp]+j]}: j in [0..n-k-2] ];
for bset in Subsets(bchoice, n-k) do
currentrow:=currentrow+1;
blist:=Sort(SetToSequence(bset));
b:=[];
for coordinate in [1..n-k] do
addition:=blist[coordinate];
for numbers in [1..coordinate-1] do
addition:=MinusOne(addition,grp);
end for;
b:=b cat [addition];
end for;

j:=Position(a cat b, grp);
M[currentrow,j]:=1;

a1:=[(a[1][i]-a[2][i]) mod grp[i]: i in [1..#grp]];
a2:=[(a[2][i]-a[1][i]) mod grp[i]: i in [1..#grp]];
if a[1] eq a[2] then
anew:=[a1,a[2]];
j:=Position(anew cat b, grp);
M[currentrow,j]:=M[currentrow,j]-1;
else
anew:=[a1,a[2]];
j:=Position(anew cat b, grp);
M[currentrow,j]:=M[currentrow,j]-1;
anew:=[a[1],a2];
j:=Position(anew cat b, grp);
M[currentrow,j]:=M[currentrow,j]-1;
end if;
end for;
end for;

return M;
end function;




function Sol2(n,p,P)
N:=&*[p[i]: i in [1..#p]];
return Binomial(N-1+n,n)-Rank(ChangeRing(Eshort(n,p), FiniteField(NextPrime(P))));
end function;





function Check(a,m,P);
if P ne "Q" then
field:=FiniteField(StringToInteger(P));
else
field:=RationalField();
end if;
n:=#a;//Read the dimension for B
p:=&*[m[i]:i in [1..#m]];
d:=Binomial(p-1+n,n);//Find the dimension
v:=Matrix(field, 1, d, [<1,Position(a,m),1>])[1];//Create the vector for a 
E:=ChangeRing(Matrix(Eshort(n,m)), field);
r:=RowSpace(E);//Get the relation space
R:=RSpace(field,d);
Bn,f:=quo<R|r>;//Take the quotient and the map to get Bn
return f(v);//return the image of v in Bn
end function;




/* Sample running 
Check([[6],[9],[12]],[16],"Q");
Check([[10],[12],[13]],[16],"Q"); 
Check([[1,0],[0,1]],[2,2],"2"); 
*/