diff --git a/Polynomial/test/Polynomial/Interpolator.cpp b/Polynomial/test/Polynomial/Interpolator.cpp
new file mode 100644
index 00000000000..1058a5dcd3b
--- /dev/null
+++ b/Polynomial/test/Polynomial/Interpolator.cpp
@@ -0,0 +1,67 @@
+#include <iostream>
+#include <CGAL/basic.h>
+#include <cassert>
+
+#include <CGAL/Arithmetic_kernel.h>
+#include <CGAL/Modular.h>
+#include <CGAL/Polynomial.h>
+#include <CGAL/Polynomial_traits_d.h>
+#include <CGAL/Random.h>
+#include <CGAL/Timer.h>
+
+#include <CGAL/gen_sparse_polynomial.h>
+#include <CGAL/Interpolator.h>
+
+static CGAL::Random my_rnd(346); // some seed 
+
+template<class Polynomial_d> 
+void test_interpolate(){
+    typedef typename CGAL::Polynomial_traits_d<Polynomial_d> PT;
+    typedef typename PT::Innermost_coefficient IC; 
+    typedef typename PT::Coefficient Coeff; 
+   
+    for(int k = 0 ; k < 5; k++){
+        Polynomial_d F = 
+            CGAL::generate_sparse_random_polynomial<Polynomial_d>
+            (my_rnd, 20/PT::d);
+        
+        typedef CGAL::Polynomial_traits_d<Polynomial_d> PT;
+        typedef std::pair<IC,Coeff> Point; 
+        std::vector<Point> points;
+        for(int i = 0; i <= F.degree(); i++){
+            points.push_back(Point(IC(i),typename PT::Evaluate()(F,IC(i))));
+        }
+        CGAL::Interpolator<Polynomial_d> 
+            interpolator(points.begin(),points.end());
+        Polynomial_d F_new = interpolator.get_interpolant();
+        assert(F_new == F);
+    }
+    std::cout << "end interpolate: " << PT::d << std::endl;
+}
+
+int main(){
+
+    CGAL::set_pretty_mode(std::cout);
+
+    typedef CGAL::Arithmetic_kernel AK;
+    typedef AK::Integer Integer; 
+    typedef CGAL::Polynomial<Integer>      Polynomial_1;
+    typedef CGAL::Polynomial<Polynomial_1> Polynomial_2;
+    typedef CGAL::Polynomial<Polynomial_2> Polynomial_3;
+    
+
+    typedef CGAL::Polynomial<CGAL::Modular>  MPolynomial_1;
+    typedef CGAL::Polynomial<MPolynomial_1>  MPolynomial_2;
+    typedef CGAL::Polynomial<MPolynomial_2>  MPolynomial_3;
+
+
+    // test_resultant<Polynomial_3>();
+
+    test_interpolate<Polynomial_1>();
+    test_interpolate<Polynomial_2>();
+    test_interpolate<Polynomial_3>();
+    test_interpolate<MPolynomial_1>();
+    test_interpolate<MPolynomial_2>();
+    test_interpolate<MPolynomial_3>();
+}
+
diff --git a/Polynomial/test/Polynomial/include/CGAL/Interpolator.h b/Polynomial/test/Polynomial/include/CGAL/Interpolator.h
new file mode 100644
index 00000000000..6791b0950a3
--- /dev/null
+++ b/Polynomial/test/Polynomial/include/CGAL/Interpolator.h
@@ -0,0 +1,107 @@
+#ifndef CGAL_INTERPOLATE_H
+#define CGAL_INTERPOLATE_H
+
+CGAL_BEGIN_NAMESPACE
+
+// Class for interpolation of univariate or multivariate polynomials. 
+// The template argument must be a model of concept Polynomial_d
+// 
+//
+template <class Polynomial_d_>
+class Interpolator{
+    typedef CGAL::Polynomial_traits_d<Polynomial_d_> PT;
+    
+public:
+    typedef typename PT::Polynomial_d Polynomial_d; 
+    typedef typename PT::Coefficient Coeff; 
+    typedef typename PT::Innermost_coefficient IC;
+
+private: 
+    typedef typename CGAL::Coercion_traits<Coeff,IC>::Cast IC2Coeff;
+    typedef typename PT::Construct_polynomial Construct;
+
+    std::vector<IC> xvals; 
+    std::vector<Coeff> yvals; 
+    std::vector<Coeff> b; 
+    
+    bool valid; 
+    Polynomial_d interpolant; 
+
+
+    Coeff eval_newton(int n, IC z)
+    {
+        Coeff p(b[n]);
+        for (int i = n-1; i >=0; i--){
+            Coeff tmp(IC2Coeff()((z - xvals[i])));
+            p = p * tmp + b[i];
+        }
+        return p;
+    }
+
+
+    Polynomial_d eval_newton_poly(int n)
+    {    
+        Polynomial_d p(Construct()(b[n]));
+        for (int i = n-1; i >=0; i--) {
+            Polynomial_d tmp = Construct()(IC2Coeff()(-xvals[i]),Coeff(1));
+            p = (p * tmp) + b[i];
+        }
+        return p;
+    }
+    
+public:
+    // constructor from an InputIterator range with value type std::pair<IC,Coeff> 
+    template<class InputIterator>
+    Interpolator(const InputIterator& begin, const InputIterator& end){
+        for(InputIterator it = begin; it != end; it++){
+            add_interpolation_point(*it);
+        }
+    }
+    
+    /*Interpolator(std::vector<IC> xvals_, std::vector<Coeff> yvals_)
+        : valid(false) {
+        CGAL_precondition(xvals_.size() != 0);
+        CGAL_precondition(xvals_.size() == yvals_.size());
+        for(unsigned i = 0; i < xvals_.size();  i++){
+            add_interpolation_point(xvals_[i],yvals_[i]);
+        }
+    }*/
+    
+    // void add_interpolation_point(std::pair<IC,Coeff> point){
+    //    add_interpolation_point(point.first, point.second);
+    // }
+    
+    // void add_interpolation_point(IC xval, Coeff yval){
+    void add_interpolation_point(std::pair<IC,Coeff> point){
+        valid = false;
+//        CGAL_precondition(0 == std::count(xval, xvals.begin(), yvals.end()));
+        xvals.push_back(point.first);
+        yvals.push_back(point.second);
+        
+        Coeff num, den;
+        int k = xvals.size() - 1; 
+        if(k == 0){
+            b.push_back(yvals[0]);
+        }else{
+            num = yvals[k] - eval_newton(k-1,xvals[k]);    
+            den = Coeff(1);            
+            for (int j = 0; j < k; j++) {
+                // (k-j) if xvals's are sequential
+                den *= (xvals[k] - xvals[j]);
+            }
+            b.push_back(num / den);
+        }
+    }
+    
+    Polynomial_d get_interpolant(){ 
+        // TODO: there should be a way to compute new interpolant from old interpolant
+        if(!valid)
+            interpolant = eval_newton_poly(xvals.size()-1);
+        return interpolant; 
+    }
+    
+};
+
+CGAL_END_NAMESPACE
+
+#endif // CGAL_INTERPOLATE_H