Finding relations between invariant polynomials Suppose I have an action of a linear reductive group ($GL(2,\mathbb{C})^2$ in this case) on a complex vector space (of dimension $16$) and I want to compute explicitly the ring of invariants of this representation. I already have the candidates for the generators (certain polynomials of $16$ variables) but I need to find the relations between them.
My question is the following - is it possible to force Magma, Sage or other program to find these polynomial relations (hopefully the whole ideal of relations...) for me? I don't expect them to be complicated, but they are way too big to manage by hand. I know the topic is discussed in several books (Derksen & Kemper, Sturmfels etc.) but are there precise, already implemented algorithms doing this or should one do it by himself? 
PS. I saw Algorithms in Invariant Theory topic, but that's not exactly what I'm asking about.
 A: I once used the following MAGMA-code (just for finite groups, I don't know if it works for general reductive groups).
intrinsic Presentation(R::RngInvar) -> SeqEnum
   {A presentation of the invariant ring R.}

   fund := FundamentalInvariants(R);
   prim := PrimaryInvariants(R);
   sec := IrreducibleSecondaryInvariants(R);
   invar := prim cat sec;
   P := PolynomialRing(BaseRing(R), #fund);
   A := Algebra(R);

   invarpres := [];
   for f in invar do
       b,g := HomogeneousModuleTest(fund,[R!1],f);
       Append(~invarpres, g[1]);
   end for;

   rel := RelationIdeal(R);

   phi:=hom<A->P|invarpres>;

   return ideal<P|[phi(r) : r in Basis(rel)]>;

end intrinsic;

Put this into a file, say "pres.m", and then you can do:
> Attach("pres.m");
> G:=MatrixGroup<2,Rationals() | [-1,0,0,-1] >;
> G;
MatrixGroup(2, Rational Field)
Generators:
[-1  0]
[ 0 -1]
> R:=InvariantRing(G);
> Presentation(R);
Ideal of Polynomial ring of rank 3 over Rational Field
Order: Lexicographical
Variables: $.1, $.2, $.3
    Homogeneous
    Basis:
    [
    $.1*$.3 - $.2^2
]

This is an example of the $A_2$-singularity
A: As mentioned in Sturmfel's book, it can be done in any computational environment by finding the Gröbner basis of the invariants with respect to lexicographic ordering. 
For instance, if you have invariants $$f_1=x_1^2+x_2^2,$$ $$f_2=x_1^2x_2^2,$$ $$f_3=x_1^3x_2+x_2^3x_1,$$ finding a Gröbner basis of the ideal $$I=\langle f_1-x_1^2+x_2^2, \; f_2-x_1^2x_2^2,\; f_3-x_1^3x_2+x_2^3x_1 \rangle$$ with lexicographic ordering $x_1<x_2<x_3<f_1<f_2<f_3$ yields the first generator $f_1^2f_2-4f_2^2-f_3^2$ which is the syzygy you are looking for.  Now, if you have other invariant polynomials that you want to rewrite in the fundamental invariants $f_1,f_2$ and $f_3$, it can be done by determining their normal form reduction modulo the Gröbner basis. For example, $x_1^7x_2-x_2^7x_1 \rightarrow f_1^2f_3-f_2f_3$.
