n := 31847;
k := FF(31847,1);
X := List([413,432,451,470,489,508,527,546,565,584], i-> FFElt(k,i));
Y := List([25439,14847,24780,5910,12734,12492,12555,28578,20806,21462], i-> FFElt(k,i));
kx := PolyAlg(k, "x");



MakeFun := function(p)

   local L,K,f;
       
      L :=PolyToList(p);
      K := List([1..Length(L)], i-> [L[i], Length(L)-i]);
      f := function(x)
        return(Sum(List(K, i-> x^(i[2])*i[1])));
      end;
      return(f);
end;



LagrangePoly := function(X,Y)
 
  local PX,PList, PKoeff, p;

     PX := List(X, i-> (x-i));
     PList := List([1..Length(PX)], i-> Product(Concatenation(List([1..i-1], j-> PX[j]), List([i+1..Length(PX)], j-> PX[j]))));
     PList := List(PList,k->  PolyMove(k,kx));
     PKoeff := List([1..Length(X)], i-> MakeFun(PList[i])(X[i]));
     PList := List([1..Length(X)], i-> (1/(PKoeff[i]))*PList[i]);
     p := Sum(List([1..Length(PList)], i-> PList[i]*Y[i]));
     return(p);

  end;

p := LagrangePoly(X,Y);
Print("\np(X) = ",p,"\n\n");
Print("Ueberpruefen: \n");
g := MakeFun(p);
Print("---------------------------------------------\n");
Print("i     |   X_i       |  p(X_i)    |    Y_i \n");
Print("---------------------------------------------\n");
for xx in [1..Length(X)] do
   Print(xx,"         ",X[xx], "           ",g(X[xx]),"          ",Y[xx]," \n");
od;