To describe the solution of the linear system in GAP:
The first step is to get linear equations by evaluating the polynomial at the points. To allow distinguishing the coefficients we introduce the variables as extra polynomial variables:
field:=GF(5);
x:=X(field,"x");
y:=X(field,"y");
z:=X(field,"z");
a:=X(field,"a");
b:=X(field,"b");
c:=X(field,"c");
d:=X(field,"d");
e:=X(field,"e");
f:=X(field,"f");
g:=X(field,"g");
i:=X(field,"i");
k:=X(field,"k");
l:=X(field,"l");
pol:=a*x^3+b*y^3+c*z^3+d*x^2*y+e*x^2*z+f*y^2*x+i*y^2*z+g*z^2*x
+k*z^2*y+l*x*y*z;
Next, define the points, move them over GF(5), and evaluate the polynomial:
pts:=[[1,2,0],[2,0,0],[1,4,1],[-2,-1,-1],[-3,2,1],[1,-1,0],
[2,-3,1],[-3,-2,2],[1,-2,1]];
pts:=List(pts,x->x*One(field));
vals:=List(pts,p->Value(pol,[x,y,z],p));
The result you get are linear equations, but written as polynomials. While it probably is quickest to retype them, the following GAP code constructs the coefficient matrix for these equations. (The list mums
gives the index numbers of the variables):
nums:=List([a,b,c,d,e,f,g,i,k,l],IndeterminateNumberOfUnivariateLaurentPolynomial);
mat:=[];
for v in vals do
ex:=ExtRepPolynomialRatFun(v);
vec:=List([1..Length(nums)],x->Zero(field));
for s in [1,3..Length(ex)-1] do
vec[Position(nums,ex[s][1])]:=ex[s+1];
od;
Add(mat,vec);
od;
Since GAP works with row vectors, we need to have the equation coefficients actually in the columns, thus transpose the matrix. We now calculate a basis for the nullspace, and construct the corresponding polynomials:
mat:=TransposedMat(mat);
sol:=[];
for v in TriangulizedNullspaceMat(mat) do
Add(sol,Value(pol,[a,b,c,d,e,f,g,i,k,l],v));
od;
We get the following three basis elements for the solutions:
$Z(5)^{3}x^{2}y-xy^{2}+xyz+y^{3}-y^{2}z+Z(5)^{3}yz^{2}$,
$Z(5)xyz+Z(5)^{3}y^{2}z+Z(5)yz^{2}+z^{3}$,
$x^{2}z+xyz+Z(5)xz^{2}-y^{2}z+yz^{2}$