Algorithm for expressing Gröbner basis in terms of ideal generators

Question

Given a polynomial ring and an ideal $$A \supset I = (f_1, ..., f_m)$$ there are plenty of implementations of an algorithm (e.g. Buchberger's) that produces a Gröbner basis $$G = (g_1, ..., g_n)$$

and additionally decide whether another element $a\in A$ is contained in $I$ and how to express it in terms of the Gröbner basis: $$a = \sum_i a_ig_i$$

Sympy does this by Q, r = G.reduce(f)

What I am looking for is an algorithm to express the elements of the Gröbner basis in terms of the ideal generators:

$$g_i = \sum_j b_{ij} f_j$$

Then by composition one could express elements in terms of the ideal generators.

Presumably it would not be too hard to modify the existing Gröbner algorithms to do this, but I am hoping for a computer algebra package in which someone has already done it, and preferably one that is free and open source.

[Edit]: I have attempted to implement the method in ahulpke's comment, but it does not work in this case:

x,y,z, z0,z1,z2 = S.symbols('x,y,z,z0,z1,z2')

A = 3*x**2 + 1
B = 2*y**2
C = (x**3+1) - (y**2)

G = S.groebner([C, A, B], order='lex')     # These generate the unit ideal
G0 = S.groebner([C+z0, A+z1, B+z2], order='lex')
Q0, r0 = G0.reduce(x**0)    # ([0, 0, 0, 0, 0], 1)
```

Answer 1

As an answer, since it is a bit too long for a comment. Also (the system I'm most familiar with -- sorry!) in GAP, which admittedly is not the optimal system for Grobner bases.

No, you don't need to decompose in the extended Grobner basis, but you identify those generators that after specializing the slack variables to 0 give you the Grobner basis. For example:

Set up the problem

gap> r:=PolynomialRing(Rationals,["x","y","z","z0","z1","z2"]);
Rationals[x,y,z,z0,z1,z2]
gap> AssignGeneratorVariables(r);
#I  Assigned the global variables [ x, y, z, z0, z1, z2 ]
gap> A:=3*x^2+1;
3*x^2+1
gap> B:=3*y^2;
3*y^2
gap> C:=x^3+1-y^2;
x^3-y^2+1

Calculate Grobner basis with and without extra variables:

gap> ReducedGrobnerBasis([A,B,C],MonomialLexOrdering());
[ 1 ]
gap> red:=ReducedGrobnerBasis([A-z0,B-z1,C-z2],MonomialLexOrdering());
[ z0^3-3*z0^2-3*z1^2-18*z1*z2-27*z2^2+3*z0+18*z1+54*z2-28, y^2-1/3*z1,
  x*z1+3*x*z2-1/3*z0^2-3*x+2/3*z0-1/3, x*z0-x-z1-3*z2+3, x^2-1/3*z0+1/3 ]

Now see which of the terms give us (up to scaling) the basis elements we wanted (i.e. set the slack variables to zero)

gap> List(red,x->Value(x,[z0,z1,z2],[0,0,0]));
[ -28, y^2, -3*x-1/3, -x+3, x^2+1/3 ]

Now we can read off:

gap> A^3-3*A^2-3*B^2-18*B*C-27*C^2+3*A+18*B+54*C;
28

That is the Grobner basis generator $1$ can be written as $$ \frac{1}{28}(A^3-3A^2-3B^2-18BC-27C^2+3A+18B+54C) $$