ring rr = (2,X), (A), lp;
minpoly = X^3 + X + 1;

ideal vec_A = 0, 1, X, X+1, X^2, X^2+1, X^2+X, X^2+X+1;
ideal vec_Z = 0, 1, X^2+X+1, X^2+X+1, X^2+1, X+1, X^2+1, X^2+1;
/* (vec_A[i], vec_Z[i]) is a pair of (x,y) coordinates */

poly f = 0; /* final function */
poly prod; /* recording term of product for each (x,y) */
int i, k;

for(k = 1; k <= 8; k = k + 1)
{
  prod = 1;
  /* Initialization */
  for(i = 1; i <= 8; i = i + 1)
  {
    if(i != k)
    {
      prod = prod*(A - vec_A[i])/(vec_A[k] - vec_A[i]);
      /* term \frac{\prod_{i \neq k} (x - x_i)}{\prod_{i \neq k} (x_k - x_i)} */
    }
  }
  f = f + prod*vec_Z[k];
  /* term \sum_{k=1}^N prod \cdot y_k */
}
"Function after interpolation: ",f;