C++: Lagrange Polynomial interpolation to interpolate polynomial defined over a field

Question

I need a C++ source code for Lagrange interpolation in a field. So given n pairs(x,y) I can construct a polynomial over a field. My code is as below but it does not output a correct result:

bigint* interpolate (int n, bigint *x, bigint *y,bigint N) {
  int d,i;
  bigint temp,temp2,temp3;
  for (d = 1; d < n; d++) {
    for (i = 0; i < n - d; i++) {
      int j = i + d;

      mul( temp,x[j], y[i]);
      mpz_mul( temp2,x[i], y[i+1] );
      sub( temp2,temp2,temp );
      sub( temp3, x[j], x[i] ); temp3=x[j]-x[i]
      invert(temp3,temp3,N); 
      mul(y[i] ,temp2, temp3 );// y[i]=temp2 * temp3
    }
  }
  return y;
}

Here N is modulus, n is number of pairs, x and y are arrays of big integers containing x_i and f(x_i)=y_i elements.