An element $\notin I$, any ideal is projected to zero. I have rwo unrelated questions but I put them in one question because they occur in the same GAP session:
In following code $x, y, z$ are already defined as indeterminates in the rationals.
In following code:
gap> P := PolynomialRing(Rationals, ["x","y","z"]);;
gap> AssignGeneratorVariables(P);
#I  Assigned the global variables [ x, y, z ]
gap> pols:=[ x^3-3*x-1, x^2+x*y+y^2-3, x+y+z ];
[ x^3-3*x-1, x^2+x*y+y^2-3, x+y+z ]
gap> I := Ideal(P, pols);;
gap> pr := NaturalHomomorphismByIdeal(P, I);
[ x, y, z, 1 ] -> [ 0*(1), (y), (z), (1) ]
gap> Image(pr, x);
0*(1)
gap> x in I;
false
gap> pols[1] in I;
true



*

*Question 1:


How is it possible that an element not in the ideal is projected to zero? 
One could even conclude that since $x^3-3x-1 \in I \implies 1 \in I$
In the same session:
gap> Q := Image(pr);
<ring Rationals,(1),(z),(z2),(y),(yz),(yz2)>
gap> Dimension(Q);   
6
gap> ci := CentralIdempotentsOfAlgebra(Q);;
gap> List(ci, c->Representative(PreImages(pr,c)));
[ 1/3*y^2*z+1/3*y-1/3*z+2/3, -1/3*y^2*z-1/3*y+1/3*z+1/3 ]
gap> p := last[1];;
gap> p in I; p^2-p in I;
false
true



*

*Question 2:


How do I show to somebody else that for $p = 1/3 -z/3 + y/3 - yz^2/3$ we can express $p^2-p$ can be expressed as a combination of generators of I?
 A: The first issue is a bug in NaturalHomomorphismByIdeal. (Apparently this code has been used very little so far.) It will be fixed in the next major release of GAP, I also put file with the correction at
https://www.dropbox.com/s/8dm52o7tglbnmla/polhomfix.g?dl=0
(simply read in the file, this fixes the issue).
It now produces
gap> pr := NaturalHomomorphismByIdeal(P, I);
[ x, y, z, 1 ] -> [ (-1)*(z)+(-1)*(y), (y), (z), (1) ]

For second question, you can use PolynomialDivisionAlgorithm to express elements in terms of a Gröbner basis of the ideal:
gap> b:=p^2-p;
1/9*x^4*y^2+2/9*x^3*y-2/9*x^2*y^2+1/9*x^2*y+1/9*x^2-2/9*x*y+1/9*y^2+1/9*x1/9*\y-2/9
gap> bas:=StoredGroebnerBasis(I);
[ [ x+y+z, y^2+y*z+z^2-3, z^3-3*z-1 ], MonomialGrlexOrdering() ]
gap> PolynomialDivisionAlgorithm(b,bas[1],bas[2]);
[ 0,
  [      1/9*x^3*y^2-1/9*x^2*y^3-    1/9*x^2*y^2*z+1/9*x*y^4+2/9*x*y^3*z+1/9*x*y^2*z^2-1/9*y^5-1/3*y^4*z-1/3*y^3*z^2-  1/9*y^2*z^3+2/9*x^2*y-4/9*x*y^2-2/9*x*y*z+4/9*y^3\
+2/3*y^2*z+2/9*y*z^2+1/9*x*y-1/9*y^2-1/9*y*z+1/9*x-1/3*y-1/9*z+1/9,
  1/9*y^4+1/3*y^3*z+2/9*y^2*z^2-1/9*y*z^3-1/9*y^2-1/9*z^2+1/9*y+1/9*z+1/9,
  1/9*y*z^2-1/9*y+1/9*z-1/9 ] ]
gap> c:=c[2];;
gap> gb:=bas[1];
[ x+y+z, y^2+y*z+z^2-3, z^3-3*z-1 ]
gap> b=c[1]*gb[1]+c[2]*gb[2]+c[3]*gb[3];
true

To trace through how the Groebner basis is computed, you could use
gap> SetInfoLevel(InfoGroebner,3);
gap> ReducedGroebnerBasis(pols,bas[2]);
#I  Spol(2,1)=x^2*y+x*y^2+1
#I  reduces to -y^3+3*y+1
#I  |bas|=4, 5 pairs left
#I  Spol(3,1)=x^2*y+x^2*z+3*x+1
[...]

Though admittedly that is a bit messy. You would have to copy the Groebner basis code to track through the changes done. This is currently not available.
(If the set of polynomials is no a Gröbner basis, the test could fail.)
