diff --git a/Polynomial/test/Polynomial/include/CGAL/new_resultant.h b/Polynomial/test/Polynomial/include/CGAL/new_resultant.h
index 61f231a20e9..0d7f9d7e2f8 100644
--- a/Polynomial/test/Polynomial/include/CGAL/new_resultant.h
+++ b/Polynomial/test/Polynomial/include/CGAL/new_resultant.h
@@ -10,21 +10,21 @@
 CGAL_BEGIN_NAMESPACE
 
 template <class Coeff> 
-inline Coeff new_resultant( CGAL::Polynomial<Coeff>, CGAL::Polynomial<Coeff>, int);
+inline Coeff new_resultant( const CGAL::Polynomial<Coeff>&, const CGAL::Polynomial<Coeff>&, int);
 namespace CGALi{
 
 template <class Coeff> 
-inline Coeff new_resultant_interpolate( CGAL::Polynomial<Coeff>, CGAL::Polynomial<Coeff> );
+inline Coeff new_resultant_interpolate( const CGAL::Polynomial<Coeff>&, const CGAL::Polynomial<Coeff>& );
 template <class Coeff> 
-inline Coeff new_resultant_modularize( CGAL::Polynomial<Coeff>, CGAL::Polynomial<Coeff>, CGAL::Tag_true);
+inline Coeff new_resultant_modularize( const CGAL::Polynomial<Coeff>&, const CGAL::Polynomial<Coeff>&, CGAL::Tag_true);
 template <class Coeff> 
-inline Coeff new_resultant_modularize( CGAL::Polynomial<Coeff>, CGAL::Polynomial<Coeff>, CGAL::Tag_false);
+inline Coeff new_resultant_modularize( const CGAL::Polynomial<Coeff>&, const CGAL::Polynomial<Coeff>&, CGAL::Tag_false);
 template <class Coeff> 
-inline Coeff new_resultant_decompose( CGAL::Polynomial<Coeff>, CGAL::Polynomial<Coeff>, CGAL::Tag_true);
+inline Coeff new_resultant_decompose( const CGAL::Polynomial<Coeff>&, const CGAL::Polynomial<Coeff>&, CGAL::Tag_true);
 template <class Coeff> 
-inline Coeff new_resultant_decompose( CGAL::Polynomial<Coeff>, CGAL::Polynomial<Coeff>, CGAL::Tag_false);
+inline Coeff new_resultant_decompose( const CGAL::Polynomial<Coeff>&, const CGAL::Polynomial<Coeff>&, CGAL::Tag_false);
 template <class Coeff> 
-inline Coeff new_resultant_( CGAL::Polynomial<Coeff>, CGAL::Polynomial<Coeff>);
+inline Coeff new_resultant_( const CGAL::Polynomial<Coeff>&, const CGAL::Polynomial<Coeff>&);
 
 } // namespace CGALi
 
@@ -32,7 +32,7 @@ namespace CGALi{
 
 template <class IC> 
 inline IC 
-new_resultant_interpolate( CGAL::Polynomial<IC> F, CGAL::Polynomial<IC> G){
+new_resultant_interpolate( const CGAL::Polynomial<IC>& F, const CGAL::Polynomial<IC>& G){
     CGAL_precondition(CGAL::Polynomial_traits_d<CGAL::Polynomial<IC> >::d == 1);
     return CGALi::resultant(F,G); 
 }
@@ -40,9 +40,8 @@ new_resultant_interpolate( CGAL::Polynomial<IC> F, CGAL::Polynomial<IC> G){
 template <class Coeff_2> 
 inline
 CGAL::Polynomial<Coeff_2>  new_resultant_interpolate(
-        CGAL::Polynomial<CGAL::Polynomial<Coeff_2> > F, 
-        CGAL::Polynomial<CGAL::Polynomial<Coeff_2> > G, 
-        int index = 0 ){
+        const CGAL::Polynomial<CGAL::Polynomial<Coeff_2> >& F, 
+        const CGAL::Polynomial<CGAL::Polynomial<Coeff_2> >& G){
     
     typedef CGAL::Polynomial<Coeff_2> Coeff_1;
     typedef CGAL::Polynomial<Coeff_1> POLY;
@@ -50,59 +49,65 @@ CGAL::Polynomial<Coeff_2>  new_resultant_interpolate(
     typedef typename PT::Innermost_coefficient IC; 
 
     CGAL_precondition(PT::d >= 2);
-
-    CGAL_precondition(index >= 0);
-    CGAL_precondition(index <  PT::d);
-        
-    POLY F_ = typename PT::Move()(F, index, 0);
-    POLY G_ = typename PT::Move()(G, index, 0);
-    CGAL_postcondition(index != 0 || F_ == F);
     
     typename PT::Degree degree; 
-    int maxdegree = degree(F_,0)*degree(G_,PT::d-1) + degree(F_,PT::d-1)*degree(G_,0); 
+    typename CGAL::Polynomial_traits_d<Coeff_1>::Degree_vector degree_vector; 
+
+    int maxdegree = degree(F,0)*degree(G,PT::d-1) + degree(F,PT::d-1)*degree(G,0); 
 
     typedef std::pair<IC,Coeff_2> Point; 
     std::vector<Point> points; // interpolation points  
-
-    typename CGAL::Polynomial_traits_d<Coeff_1>::Degree_vector degree_vector; 
-
+    
+   
     int i(0);
     CGAL::Exponent_vector ev_f(degree_vector(Coeff_1()));
     CGAL::Exponent_vector ev_g(degree_vector(Coeff_1()));
+    
+   
     while((int) points.size() <= maxdegree + 1){
         i++;
-        Coeff_1 F_at_i = typename PT::Evaluate()(F_,IC(i));
-        Coeff_1 G_at_i = typename PT::Evaluate()(G_,IC(i));
-        if(degree_vector(F_at_i) >  ev_f || degree_vector(G_at_i) >  ev_g){
+        // timer1.start();
+        Coeff_1 Fat_i = F.evaluate(Coeff_1(i));
+        Coeff_1 Gat_i = G.evaluate(Coeff_1(i));
+        // timer1.stop();
+        
+        if(degree_vector(Fat_i) >  ev_f || degree_vector(Gat_i) >  ev_g){
             points.clear();
-            ev_f  = degree_vector(F_at_i);
-            ev_g  = degree_vector(G_at_i);
+            ev_f  = degree_vector(Fat_i);
+            ev_g  = degree_vector(Gat_i);
             CGAL_postcondition(points.size() == 0);
         }
-        if(degree_vector(F_at_i) ==  ev_f && degree_vector(G_at_i) ==  ev_g){
-            Coeff_2 res_at_i = new_resultant(F_at_i, G_at_i, 0);
+        if(degree_vector(Fat_i) ==  ev_f && degree_vector(Gat_i) ==  ev_g){
+            // timer2.start();
+            Coeff_2 res_at_i = new_resultant_interpolate(Fat_i, Gat_i);
+            // timer2.stop();
             points.push_back(Point(IC(i),res_at_i));
-        }
+        }      
     }
-    Coeff_1 result = CGAL::Interpolator<Coeff_1>(points.begin(),points.end()).get_interpolant();
-    
+   
+    // timer3.start();
+    CGAL::Interpolator<Coeff_1> interpolator(points.begin(),points.end());
+    Coeff_1 result = interpolator.get_interpolant();
+    // timer3.stop();
+
 #ifndef CGAL_NDEBUG
     while((int) points.size() <= maxdegree + 3){
         i++;
-        Coeff_1 F_at_i = typename PT::Evaluate()(F_,IC(i));
-        Coeff_1 G_at_i = typename PT::Evaluate()(G_,IC(i));
+        Coeff_1 Fat_i = typename PT::Evaluate()(F,IC(i));
+        Coeff_1 Gat_i = typename PT::Evaluate()(G,IC(i));
         
-        assert(degree_vector(F_at_i) <= ev_f);
-        assert(degree_vector(G_at_i) <= ev_g);
+        assert(degree_vector(Fat_i) <= ev_f);
+        assert(degree_vector(Gat_i) <= ev_g);
         
-        if(degree_vector(F_at_i) ==  ev_f && degree_vector(G_at_i) ==  ev_g){
-            Coeff_2 res_at_i = new_resultant(F_at_i, G_at_i, 0);
+        if(degree_vector(Fat_i) ==  ev_f && degree_vector(Gat_i) ==  ev_g){
+            Coeff_2 res_at_i = new_resultant(Fat_i, Gat_i, 0);
             points.push_back(Point(IC(i), res_at_i));
         }
     }
+    CGAL::Interpolator<Coeff_1> interpolator_(points.begin(),points.end());
+    Coeff_1 result_= interpolator_.get_interpolant();
     
-    Coeff_1 result_ = CGAL::Interpolator<Coeff_1>(points.begin(), points.end()).get_interpolant();
-    // the interpolate polynomial has to be stable !
+     // the interpolate polynomial has to be stable !
     assert(result_ == result); 
 #endif 
     return result; 
@@ -112,31 +117,119 @@ CGAL::Polynomial<Coeff_2>  new_resultant_interpolate(
 template <class Coeff> 
 inline
 Coeff new_resultant_modularize( 
-        CGAL::Polynomial<Coeff> F, CGAL::Polynomial<Coeff> G, CGAL::Tag_false){
+        const CGAL::Polynomial<Coeff>& F, const CGAL::Polynomial<Coeff>& G, CGAL::Tag_false){
     return new_resultant_interpolate(F,G);
-}
+};
 
 template <class Coeff> 
 inline
 Coeff new_resultant_modularize( 
-        CGAL::Polynomial<Coeff> F, CGAL::Polynomial<Coeff> G, CGAL::Tag_true){
-    return new_resultant_interpolate(F,G);
+        const CGAL::Polynomial<Coeff>& F, const CGAL::Polynomial<Coeff>& G, CGAL::Tag_true){
+    
+    typedef Polynomial_traits_d<CGAL::Polynomial<Coeff> > PT;
+    typedef typename PT::Polynomial_d Polynomial;
+    typedef typename PT::Innermost_coefficient IC;
+    
+    
+    typedef Chinese_remainder_traits<Coeff> CRT;
+    typedef typename CRT::Scalar_type Scalar;
+
+
+    typedef typename CGAL::Modular_traits<Polynomial>::Modular_NT MPolynomial; 
+    typedef typename CGAL::Modular_traits<Coeff>::Modular_NT      MCoeff; 
+    typedef typename CGAL::Modular_traits<Scalar>::Modular_NT     MScalar;  
+        
+    typename CRT::Chinese_remainder chinese_remainder; 
+    typename CGAL::Modular_traits<Coeff>::Modular_image_inv inv_map;
+
+
+    typename PT::Degree_vector                                     degree_vector; 
+    typename CGAL::Polynomial_traits_d<MPolynomial>::Degree_vector mdegree_vector;
+
+    bool solved = false; 
+    int prime_index = 0; 
+    int n = 0;
+    Scalar p,q,pq,s,t; 
+    Coeff R, R_old; 
+    
+    CGAL::Timer timer_evaluate, timer_resultant, timer_cr; 
+    
+    do{
+        MPolynomial mF, mG;
+        MCoeff mR;
+        timer_evaluate.start();
+        do{
+            // select a prime number
+            int current_prime = -1;
+            prime_index++;
+            if(prime_index >= 2000){
+                std::cerr<<"primes in the array exhausted"<<std::endl;
+                assert(false);
+                current_prime = CGALi::get_next_lower_prime(current_prime);
+            } else{
+                current_prime = CGALi::primes[prime_index];
+            }
+            CGAL::Modular::set_current_prime(current_prime);
+            
+            mF = CGAL::modular_image(F);
+            mG = CGAL::modular_image(G);
+            
+        }while( degree_vector(F) != mdegree_vector(mF) || 
+                degree_vector(G) != mdegree_vector(mG));
+        timer_evaluate.stop();
+        
+        timer_resultant.start();
+        n++;
+        mR = new_resultant_interpolate(mF,mG);
+        timer_resultant.stop();
+        timer_cr.start();
+        if(n == 1){ 
+            // init chinese remainder
+            q =  CGAL::Modular::get_current_prime(); // implicit ! 
+            R = inv_map(mR);
+        }else{
+            // continue chinese remainder
+            p = CGAL::Modular::get_current_prime(); // implicit!  
+            R_old  = R ;
+//            chinese_remainder(q,Gs ,p,inv_map(mG_),pq,Gs);             
+//            cached_extended_euclidean_algorithm(q,p,s,t);
+            CGALi_algorithm_M::Cached_extended_euclidean_algorithm
+                <typename CRT::Scalar_type> ceea;
+            ceea(q,p,s,t);
+            pq =p*q;
+            chinese_remainder(q,p,pq,s,t,R_old,inv_map(mR),R);
+            q=pq;
+        }
+        solved = (R==R_old);
+        timer_cr.stop();       
+    } while(!solved);
+        
+    std::cout << "Time Evaluate   : " << timer_evaluate.time() << std::endl; 
+    std::cout << "Time Resultant  : " << timer_resultant.time() << std::endl; 
+    std::cout << "Time Chinese R  : " << timer_cr.time() << std::endl; 
+    // CGAL_postcondition(R == new_resultant_interpolate(F,G));
+    return R;
+    // return new_resultant_interpolate(F,G);
 }
 
 
 template <class Coeff> 
 inline
 Coeff new_resultant_decompose( 
-        CGAL::Polynomial<Coeff> F, CGAL::Polynomial<Coeff> G, CGAL::Tag_false){
+        const CGAL::Polynomial<Coeff>& F, const CGAL::Polynomial<Coeff>& G, CGAL::Tag_false){
+#ifdef CGAL_RESULTANT_NUSE_MODULAR_ARITHMETIC
+    return  new_resultant_modularize(F,G,CGAL::Tag_false());
+#else
     typedef CGAL::Polynomial<Coeff> Polynomial; 
     typedef typename Modular_traits<Polynomial>::Is_modularizable Is_modularizable; 
     return new_resultant_modularize(F,G,Is_modularizable());
+#endif
 }
 
 template <class Coeff> 
 inline
 Coeff new_resultant_decompose( 
-        CGAL::Polynomial<Coeff> F, CGAL::Polynomial<Coeff> G, CGAL::Tag_true){  
+        const CGAL::Polynomial<Coeff>& F, const CGAL::Polynomial<Coeff>& G, CGAL::Tag_true){  
     
     typedef Polynomial<Coeff> POLY;
     typedef typename Fraction_traits<POLY>::Numerator_type Numerator;
@@ -163,7 +256,7 @@ Coeff new_resultant_decompose(
 template <class Coeff> 
 inline
 Coeff new_resultant_( 
-        CGAL::Polynomial<Coeff> F, CGAL::Polynomial<Coeff> G){
+        const CGAL::Polynomial<Coeff>& F, const CGAL::Polynomial<Coeff>& G){
     typedef CGAL::Fraction_traits<Polynomial<Coeff > > FT;
     typedef typename FT::Is_fraction Is_fraction; 
     return new_resultant_decompose(F,G,Is_fraction());
@@ -175,13 +268,13 @@ Coeff new_resultant_(
 template <class Coeff> 
 inline
 Coeff  new_resultant( 
-        CGAL::Polynomial<Coeff> F, 
-        CGAL::Polynomial<Coeff> G, 
+        const CGAL::Polynomial<Coeff>& F_, 
+        const CGAL::Polynomial<Coeff>& G_, 
         int index = CGAL::Polynomial_traits_d< CGAL::Polynomial<Coeff> >::d-1){
     typedef CGAL::Polynomial_traits_d<CGAL::Polynomial<Coeff> > PT;
-    CGAL::Polynomial<Coeff> F_ = typename PT::Move()(F, index, 0);
-    CGAL::Polynomial<Coeff> G_ = typename PT::Move()(G, index, 0);
-    return CGALi::new_resultant_(F_,G_);
+    CGAL::Polynomial<Coeff> F = typename PT::Move()(F_, index, 0);
+    CGAL::Polynomial<Coeff> G = typename PT::Move()(G_, index, 0);
+    return CGALi::new_resultant_(F,G);
 }
     
 CGAL_END_NAMESPACE
diff --git a/Polynomial/test/Polynomial/resultant.cpp b/Polynomial/test/Polynomial/resultant.cpp
new file mode 100644
index 00000000000..bc62327857d
--- /dev/null
+++ b/Polynomial/test/Polynomial/resultant.cpp
@@ -0,0 +1,130 @@
+
+// #define CGAL_RESULTANT_NUSE_MODULAR_ARITHMETIC 1
+
+#include <CGAL/basic.h>
+
+#include <CGAL/Timer.h>
+CGAL::Timer timer1, timer2, timer3; 
+
+#include <cassert>
+#include <iostream>
+
+#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/gen_sparse_polynomial.h>
+#include <CGAL/new_resultant.h>
+
+
+static CGAL::Random my_rnd(346); // some seed 
+
+template<class Polynomial_d> 
+void test_resultant(){
+    
+    typedef typename CGAL::Polynomial_traits_d<Polynomial_d> PT;
+    typedef typename PT::Coefficient Coeff; 
+    for (int k = 0; k < 1; k++){
+        Polynomial_d F2 = 
+            CGAL::generate_sparse_random_polynomial<Polynomial_d>(my_rnd,12); 
+        Polynomial_d F1 = 
+            CGAL::generate_sparse_random_polynomial<Polynomial_d>(my_rnd,12);
+        Polynomial_d G1 = typename PT::Move()(F1,PT::d-1,0);
+        Polynomial_d G2 = typename PT::Move()(F2,PT::d-1,0);
+        
+        CGAL::Timer timer_new;
+        timer_new.start();
+        Coeff F = CGAL::new_resultant(F1,F2,PT::d-1);
+        // Coeff G = CGAL::new_resultant(G1,G2,0);
+        timer_new.stop();
+
+        assert(F==F);
+        std::cout <<"timer1: "<< timer1.time() << std::endl;
+        std::cout <<"timer2: "<< timer2.time() << std::endl;
+        std::cout <<"timer3: "<< timer3.time() << std::endl;
+        timer1.reset();
+        timer2.reset();
+        timer3.reset();
+        std::cout <<"new time: "<< timer_new.time() << std::endl;
+    }
+    std::cout << "end resultant: " << PT::d << std::endl;
+}
+   
+
+template<class Coeff> 
+void test_multi_resultant(){
+    
+    typedef CGAL::Polynomial<Coeff>        Polynomial_1;
+    typedef CGAL::Polynomial<Polynomial_1> Polynomial_2;
+    typedef CGAL::Polynomial<Polynomial_2> Polynomial_3;
+
+    Polynomial_3 Q  =
+        CGAL::generate_polynomial_degree_each<Polynomial_3>(my_rnd,4,100); 
+    Polynomial_2 B1 =
+        CGAL::generate_polynomial_degree_each<Polynomial_2>(my_rnd,6,80); 
+    Polynomial_2 B2 = 
+        CGAL::generate_polynomial_degree_each<Polynomial_2>(my_rnd,6,80); 
+    std::cout << "Q:  "<< Q   << std::endl;
+    std::cout << "B1: "<< B1 << std::endl;
+    std::cout << "B2: "<< B2 <<std::endl;
+    
+    CGAL::Timer timer_new;
+    timer_new.start();
+    Polynomial_2 R1 = CGAL::new_resultant(Q,Polynomial_3(B1),0);
+    Polynomial_1 R2 = CGAL::new_resultant(R1,B2,0);
+    std::cout << std::endl;
+    std::cout << R1.degree() << std::endl;
+    std::cout << R2.degree() << std::endl;
+    timer_new.stop();
+    
+    std::cout <<"timer1: "<< timer1.time() << std::endl;
+    std::cout <<"timer2: "<< timer2.time() << std::endl;
+    std::cout <<"timer3: "<< timer3.time() << std::endl;
+    timer1.reset();
+    timer2.reset();
+    timer3.reset();
+
+    std::cout <<"new time: "<< timer_new.time() << std::endl;
+    std::cout << "end multi resultant: " << std::endl;
+}
+   
+
+int main(){
+
+    //CGAL::force_ieee_double_precision();
+    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;
+
+    test_multi_resultant<Integer>();
+    
+    // test_resultant<Polynomial_1>();
+    // test_resultant<Polynomial_2>();
+    //test_resultant<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<MPolynomial_1>();
+    //test_resultant<MPolynomial_2>();
+    //test_resultant<MPolynomial_3>();
+
+    typedef CGAL::Polynomial<AK::Rational>   RPolynomial_1;
+    typedef CGAL::Polynomial<RPolynomial_1>  RPolynomial_2;
+    typedef CGAL::Polynomial<RPolynomial_2>  RPolynomial_3;
+
+    //test_resultant<RPolynomial_1>();
+    //test_resultant<RPolynomial_2>();
+    //test_resultant<RPolynomial_3>();
+}
+
+ 
+