diff --git a/Polynomial/include/CGAL/Polynomial_traits_d.h b/Polynomial/include/CGAL/Polynomial_traits_d.h
index 438fc18d9c0..a059bebf942 100644
--- a/Polynomial/include/CGAL/Polynomial_traits_d.h
+++ b/Polynomial/include/CGAL/Polynomial_traits_d.h
@@ -422,61 +422,36 @@ public:
       return typename CT::Cast()(p);
     }
   };
-};
 
-// Evaluate_homogeneous_func for recursive homogeneous evaluation of a
-//  polynomial, used by Polynomial_traits_d_base for polynomials.
-template< class Polynomial, int d = CGAL::Polynomial_traits_d< Polynomial>::d >
-struct Evaluate_homogeneous_func;
+  struct Substitute_homogeneous{
+  public:
+    // this is the end of the recursion
+    // begin contains the homogeneous variabel 
+    // hdegree is the remaining degree 
+    template <class Input_iterator>
+    typename 
+    CGAL::Coercion_traits<
+        typename std::iterator_traits<Input_iterator>::value_type,
+                                   Innermost_coefficient>::Type
+    operator()(
+        const Innermost_coefficient& p, 
+        Input_iterator begin, 
+        Input_iterator end,
+        int hdegree){ 
+      
+      typedef typename std::iterator_traits<Input_iterator>::value_type 
+        value_type;
+      typedef CGAL::Coercion_traits<Innermost_coefficient,value_type> CT;
+      typename CT::Type result =   
+        typename CT::Cast()(CGAL::ipower(*begin++,hdegree)) 
+        * typename CT::Cast()(p);
 
-template< class Polynomial >
-struct Evaluate_homogeneous_func< Polynomial, 1 > {
-  typedef typename CGAL::Polynomial_traits_d< Polynomial > PT;
-  typedef typename PT::Coefficient Coefficient;
-  typedef typename PT::Innermost_coefficient ICoeff;
-  typedef typename CGAL::Polynomial_traits_d< Coefficient > PTC;
-        
-  template< class Input_iterator >
-  ICoeff operator()( const Polynomial& p,
-      Input_iterator begin,
-      Input_iterator end,
-      int total_degree,
-      const ICoeff& v ) const {
-    // TODO: remove --end in case of Input_iterator
-    --end;
-    CGAL_precondition( begin == end );
-    return p.evaluate_homogeneous( (*end), v, total_degree );
-  }
-        
-}; 
-
-template< class Polynomial, int d >
-struct Evaluate_homogeneous_func {
-  typedef typename CGAL::Polynomial_traits_d< Polynomial > PT;
-  typedef typename PT::Coefficient Coefficient;
-  typedef typename PT::Innermost_coefficient ICoeff;
-  typedef typename CGAL::Polynomial_traits_d< Coefficient > PTC;
-        
-  template< class Input_iterator >
-  ICoeff operator()( const Polynomial& p,
-      Input_iterator begin,
-      Input_iterator end,
-      int total_degree,
-      const ICoeff& v ) const {
-    CGAL_precondition( begin != end );
-    //typename PT::Evaluate evaluate;
-    typename PT::Degree degree;            
-    Evaluate_homogeneous_func< Coefficient > eval_hom;
-    // TODO: remove --end in case of Input_iterator
-    --end;
-            
-    std::vector< ICoeff > cv;  
-    for( int i = 0; i <= degree(p); ++i ) {
-      cv.push_back( eval_hom( p[i], begin, end, total_degree - i, v ) );
+      CGAL_precondition(end == begin);
+      CGAL_precondition(hdegree >= 0);
+      return result;
     }
-    
-    return (CGAL::Polynomial< ICoeff >( cv.begin(), cv.end() )).evaluate((*end));
-  }
+  };
+  
 };
 
 // Now the version for the polynomials with all functors provided by all 
@@ -862,78 +837,67 @@ public:
     }
   };
 
-  //       Evaluate;
+  // Evaluate;
   struct Evaluate
-    :public Unary_function<Polynomial_d,Innermost_coefficient>{
+    :public Unary_function<Polynomial_d,Coefficient>{
+    // Evaluate with respect to one variable 
     Coefficient
-    operator()(const Polynomial_d& p, Innermost_coefficient x, int i = (d-1)) 
-      const {
-      if(i == (d-1) )
-        return p.evaluate(x);
-      else{
-        return Move()(p,i).evaluate(x);
-      }
+    operator()(const Polynomial_d& p, Coefficient x) const {
+      return p.evaluate(x);
     }
-        
-    template< class Input_iterator >
-    Innermost_coefficient 
-    operator()( 
-        const Polynomial_d& p, 
-        Input_iterator begin, 
-        Input_iterator end ) const 
+    Coefficient 
+    operator()(const Polynomial_d& p, Coefficient x, int i ) const {
+      return Move()(p,i).evaluate(x);
+    } 
+#define ICOEFF typename First_if_different<Innermost_coefficient, Coefficient>::Type 
+    Coefficient operator()
+      ( const Polynomial_d& p, const ICOEFF& x) const 
     {
-      CGAL_precondition( begin != end );
-            
-      typename PT::Evaluate evaluatePoly;
-      typename PTC::Evaluate evaluateCoeff;
-      // TODO: remove --end in case of Input_iterator
-      --end;
-      return evaluateCoeff( evaluatePoly( p, (*end) ), begin, end );            
-    }  
+      return p.evaluate(x);
+    }
+    Coefficient operator()
+      (const Polynomial_d& p, const ICOEFF& x, int i) const 
+    {
+      return Move()(p,i,PT::d-1).evaluate(x);
+    }
+#undef ICOEFF      
   };
-    
-  //       Evaluate_homogeneous;
+  
+  // Evaluate_homogeneous;
   struct Evaluate_homogeneous{
     typedef Coefficient           result_type;  
     typedef Polynomial_d          first_argument_type;
-    typedef Innermost_coefficient second_argument_type;
-    typedef Innermost_coefficient third_argument_type;
-    typedef Arity_tag< 3 >         Arity;
-        
-    Coefficient
-    operator()(
-        const Polynomial_d& p, 
-        Innermost_coefficient a, 
-        Innermost_coefficient b,
-        int i = (PT::d-1) ) const {
-      if (i == (d-1) )
-        return p.evaluate_homogeneous(a,b);
-      else
-        return Move()(p,i,PT::d-1).evaluate_homogeneous(a,b);
+    typedef Coefficient           second_argument_type;
+    typedef Coefficient           third_argument_type;
+    typedef Arity_tag< 3 >        Arity;
+       
+    Coefficient operator()(
+        const Polynomial_d& p, Coefficient a, Coefficient b) const 
+    {
+      return p.evaluate_homogeneous(a,b);
     }
-        
-    template< class Input_iterator >
-    Innermost_coefficient operator()( const Polynomial_d & p,
-        Input_iterator begin,
-        Input_iterator end ) const {
-      typename PT::Total_degree total_degree;
-      typename PT::Evaluate_homogeneous eval_hom;
-      return eval_hom( p, begin, end, total_degree(p) );
+    Coefficient operator()(
+        const Polynomial_d& p, Coefficient a, Coefficient b, int i) const 
+    {
+      return Move()(p,i,PT::d-1).evaluate_homogeneous(a,b);
     }
-        
-    template< class Input_iterator >
-    Innermost_coefficient operator()( const Polynomial_d& p,
-        Input_iterator begin,
-        Input_iterator end,
-        int total_degree ) const {
-      --end;
-      Evaluate_homogeneous_func< Polynomial_d > eval_hom;
-      return eval_hom( p, begin, end, total_degree, (*end) );
+
+#define ICOEFF typename First_if_different<Innermost_coefficient, Coefficient>::Type 
+    Coefficient operator()
+      ( const Polynomial_d& p, const ICOEFF& a, const ICOEFF& b) const 
+    {
+      return p.evaluate_homogeneous(a,b);
     }
-        
+    Coefficient operator()
+      (const Polynomial_d& p, const ICOEFF& a, const ICOEFF& b, int i) const 
+    {
+      return Move()(p,i,PT::d-1).evaluate_homogeneous(a,b);
+    }
+#undef ICOEFF    
+
   };
     
-  //       Is_zero_at;
+  // Is_zero_at;
   struct Is_zero_at {
   private:
     typedef Algebraic_structure_traits<Innermost_coefficient> AST;
@@ -946,8 +910,8 @@ public:
         const Polynomial_d& p, 
         Input_iterator begin, 
         Input_iterator end ) const {
-      typename PT::Evaluate evaluate;            
-      return( CGAL::is_zero( evaluate( p, begin, end ) ) );
+      typename PT::Substitute substitute;            
+      return( CGAL::is_zero( substitute( p, begin, end ) ) );
     } 
   };
     
@@ -960,13 +924,11 @@ public:
     typedef BOOL result_type;
         
     template< class Input_iterator >
-    BOOL operator()( 
-        const Polynomial_d& p, 
-        Input_iterator begin, 
-        Input_iterator end ) const 
+    BOOL operator()
+      ( const Polynomial_d& p, Input_iterator begin, Input_iterator end ) const 
     {
-      typename PT::Evaluate_homogeneous evaluate_homogeneous;
-      return( CGAL::is_zero( evaluate_homogeneous( p, begin, end ) ) );
+      typename PT::Substitute_homogeneous substitute_homogeneous;
+      return( CGAL::is_zero( substitute_homogeneous( p, begin, end ) ) );
     }
   };
 
@@ -984,13 +946,11 @@ public:
         Input_iterator begin, 
         Input_iterator end ) const 
     {
-      typename PT::Evaluate evaluate;
-      return CGAL::sign( evaluate( p, begin, end ) );
+      typename PT::Substitute substitute;
+      return CGAL::sign( substitute( p, begin, end ) );
     }
   };
   
-  // TODO: Fix error in evaluate homogenous
-  //       Sign_at_homogeneous;
   struct Sign_at_homogeneous {
     typedef Real_embeddable_traits<Innermost_coefficient> RT;
     typedef typename RT::Sign::result_type SIGN;
@@ -1002,8 +962,8 @@ public:
         const Polynomial_d& p, 
         Input_iterator begin, 
         Input_iterator end) const {
-      typename PT::Evaluate_homogeneous evaluate_homogeneous;
-      return CGAL::sign( evaluate_homogeneous( p, begin, end ) );
+      typename PT::Substitute_homogeneous substitute_homogeneous;
+      return CGAL::sign( substitute_homogeneous( p, begin, end ) );
     }
   };
     
@@ -1344,8 +1304,7 @@ public:
         
     Polynomial_d operator()( Polynomial_d p, const Innermost_coefficient& c,
         int i = (PT::d-1) ) {
-      typename PT::Scale_homogeneous scale_homogeneous;
-            
+      typename PT::Scale_homogeneous scale_homogeneous;          
       return scale_homogeneous( p, c, Innermost_coefficient(1), i );
     }
         
@@ -1553,37 +1512,37 @@ public:
     }
   };
 
-  // substitute every variable by its new value in the iterator range
+    // substitute every variable by its new value in the iterator range
   // begin refers to the innermost/first variable
   struct Substitute{
   public:
     template <class Input_iterator>
-    typename 
-    CGAL::Coercion_traits
-    <typename std::iterator_traits<Input_iterator>::value_type, 
-         Innermost_coefficient>::Type
+    typename CGAL::Coercion_traits<
+         typename std::iterator_traits<Input_iterator>::value_type, 
+         Innermost_coefficient
+    >::Type
     operator()(
-        const Polynomial_d& p_, 
+        const Polynomial_d& p, 
         Input_iterator begin, 
         Input_iterator end){
       typedef typename std::iterator_traits<Input_iterator> ITT;
       typedef typename ITT::iterator_category  Category; 
-      return (*this)(p_,begin,end,Category()); 
+      return (*this)(p,begin,end,Category()); 
     }
         
     template <class Input_iterator>
-    typename 
-    CGAL::Coercion_traits
-    <typename std::iterator_traits<Input_iterator>::value_type, 
-         Innermost_coefficient>::Type
+    typename CGAL::Coercion_traits< 
+         typename std::iterator_traits<Input_iterator>::value_type, 
+         Innermost_coefficient
+    >::Type
     operator()(
-        const Polynomial_d& p_, 
+        const Polynomial_d& p, 
         Input_iterator begin, 
         Input_iterator end,
         std::forward_iterator_tag){
       typedef typename std::iterator_traits<Input_iterator> ITT;
       std::list<typename ITT::value_type> list(begin,end); 
-      return (*this)(p_,list.begin(),list.end()); 
+      return (*this)(p,list.begin(),list.end()); 
     }
         
     template <class Input_iterator>
@@ -1592,7 +1551,7 @@ public:
     <typename std::iterator_traits<Input_iterator>::value_type, 
          Innermost_coefficient>::Type
     operator()(
-        const Polynomial_d& p_, 
+        const Polynomial_d& p, 
         Input_iterator begin, 
         Input_iterator end,
         std::bidirectional_iterator_tag){
@@ -1602,9 +1561,7 @@ public:
       typedef CGAL::Coercion_traits<Innermost_coefficient,value_type> CT;
       typename PTC::Substitute subs; 
 
-      Polynomial_d p = p_;
-      end--;
-      typename CT::Type x = typename CT::Cast()(*end);  
+      typename CT::Type x = typename CT::Cast()(*(--end));  
             
       int i = Degree()(p);
       typename CT::Type y = 
@@ -1616,7 +1573,72 @@ public:
       }
       return y;  
     }        
+  };  // substitute every variable by its new value in the iterator range
+
+
+  
+  // begin refers to the innermost/first variable
+  struct Substitute_homogeneous{
+
+    template<typename Input_iterator>
+    struct Result_type{
+      typedef std::iterator_traits<Input_iterator> ITT;
+      typedef typename ITT::value_type value_type;
+      typedef Coercion_traits<value_type, Innermost_coefficient> CT; 
+      typedef typename CT::Type Type; 
+    };
+    
+  public:
+    
+    template <class Input_iterator>
+    typename Result_type<Input_iterator>::Type
+    operator()( const Polynomial_d& p, Input_iterator begin, Input_iterator end){
+      int hdegree = Total_degree()(p);
+
+      typedef std::iterator_traits<Input_iterator> ITT;
+      std::list<typename ITT::value_type> list(begin,end); 
+      
+      // make the homogeneous variable the first in the list
+      list.push_front(list.back());
+      list.pop_back();
+      
+      // reverse and begin with the outermost variable 
+      return (*this)(p, list.rbegin(), list.rend(), hdegree); 
+    }
+      
+    // this operator is undcoumented and for internal use:
+    // the iterator range starts with the outermost variable 
+    // and ends with the homogeneous variable 
+    template <class Input_iterator>
+    typename Result_type<Input_iterator>::Type
+    operator()(
+        const Polynomial_d& p, 
+        Input_iterator begin, 
+        Input_iterator end,
+        int hdegree){
+      
+      
+      typedef std::iterator_traits<Input_iterator> ITT;
+      typedef typename ITT::value_type value_type;
+      typedef Coercion_traits<value_type, Innermost_coefficient> CT; 
+      typedef typename CT::Type Type; 
+      
+      typename PTC::Substitute_homogeneous subsh; 
+      
+      typename CT::Type x = typename CT::Cast()(*begin++);  
+            
+
+      int i = Degree()(p);
+      typename CT::Type y = subsh(Get_coefficient()(p,i),begin,end, hdegree-i);
+       
+      while (--i >= 0){
+        y *= x;
+        y += subsh(Get_coefficient()(p,i),begin,end,hdegree-i);
+      }
+      return y;  
+    }        
   };
+
 };
 
 } // namespace CGALi
diff --git a/Polynomial/test/Polynomial/Polynomial_traits_d.cpp b/Polynomial/test/Polynomial/Polynomial_traits_d.cpp
index a1883678359..65e85d97f73 100644
--- a/Polynomial/test/Polynomial/Polynomial_traits_d.cpp
+++ b/Polynomial/test/Polynomial/Polynomial_traits_d.cpp
@@ -926,33 +926,18 @@ void test_evaluate(){
     assert(evaluate(Polynomial_d(1),ICoeff(0)) == Coeff(1));
     assert(evaluate(Polynomial_d(2),ICoeff(5)) == Coeff(2));
 
-    assert(
-            evaluate(Polynomial_d(Coeff(3),Coeff(2)),ICoeff(0)) == Coeff(3));
-    assert(
-            evaluate(Polynomial_d(Coeff(3),Coeff(2)),ICoeff(1)) == Coeff(5));
-    assert(
-            evaluate(Polynomial_d(Coeff(3),Coeff(2)),ICoeff(2)) == Coeff(7));
+    assert( evaluate(Polynomial_d(Coeff(3),Coeff(2)),ICoeff(0)) == Coeff(3));
+    assert( evaluate(Polynomial_d(Coeff(3),Coeff(2)),ICoeff(1)) == Coeff(5));
+    assert( evaluate(Polynomial_d(Coeff(3),Coeff(2)),ICoeff(2)) == Coeff(7));
     
     for(int i = 0; i < 5; i++){
         int n = my_rnd.get_int(0,PT::d-1);
         Polynomial_d p,q;
         p = generate_sparse_random_polynomial<Polynomial_d>();
-        assert(evaluate(p,ICoeff(3),n) == 
-                evaluate(move(p,n,PT::d-1),ICoeff(3)));
+        assert(evaluate(p,ICoeff(3),n) 
+            == evaluate(move(p,n,PT::d-1),ICoeff(3)));
     }
     
-    Polynomial_d p(Coeff(1), Coeff(2), Coeff(3));
-    std::vector< ICoeff > ev;
-    for(int i = 0; i < PT::d-1; ++i) {
-        ev.push_back(ICoeff(0));
-    }
-    ev.push_back(ICoeff(1));
-    assert(evaluate(p, ev.begin(), ev.end()) == ICoeff(1+2+3));
-    ev.pop_back();
-    ev.push_back(ICoeff(2));  
-    assert(evaluate(p, ev.begin(), ev.end()) == ICoeff(1+4+12));
-    
-    
     std::cerr << " ok "<< std::endl;
 }
 
@@ -967,34 +952,13 @@ void test_evaluate_homogeneous(){
     assert(evh(Polynomial_d(1),ICoeff(0),ICoeff(2)) == Coeff(1));
     assert(evh(Polynomial_d(2),ICoeff(5),ICoeff(3)) == Coeff(2));
 
-    assert(evh(Polynomial_d(Coeff(3),Coeff(2)),ICoeff(0),ICoeff(1)) ==
-            Coeff(3));
-    assert(evh(Polynomial_d(Coeff(3),Coeff(2)),ICoeff(1),ICoeff(1)) ==
-            Coeff(5));
-    assert(evh(Polynomial_d(Coeff(3),Coeff(2)),ICoeff(2),ICoeff(3)) ==
-            Coeff(9+4));
+    assert(evh( Polynomial_d(Coeff(3),Coeff(2)) , ICoeff(0),ICoeff(1)) 
+        == Coeff(3));
+    assert(evh( Polynomial_d(Coeff(3),Coeff(2)) , ICoeff(1),ICoeff(1)) 
+        == Coeff(5));
+    assert(evh( Polynomial_d(Coeff(3),Coeff(2)) , ICoeff(2),ICoeff(3)) 
+        == Coeff(9+4));
     
-    // hdegree parameter no longer available
-    /*assert(evh(Polynomial_d(Coeff(5),Coeff(7)),ICoeff(2),ICoeff(3),2)
-      == Coeff(5*3*3+7*3*2));*/
-    
-    Polynomial_d p(Coeff(1), Coeff(2), Coeff(3));
-    std::vector< ICoeff > ev;
-    for(int i = 0; i < PT::d-1; ++i) {
-        ev.push_back(ICoeff(0));
-    }
-    ev.push_back(ICoeff(1));
-    ev.push_back(ICoeff(1));
-    assert(evh(p, ev.begin(), ev.end()) == ICoeff(1+2+3));
-    ev.pop_back();
-    ev.push_back(ICoeff(2));  
-    assert(evh(p, ev.begin(), ev.end()) == ICoeff(4+4+3));
-    ev.pop_back();
-    ev.pop_back();
-    ev.push_back(ICoeff(2));  
-    ev.push_back(ICoeff(2));  
-    assert(evh(p, ev.begin(), ev.end()) == ICoeff(4+8+12));
-        
     std::cerr << " ok "<< std::endl;
 }
 
@@ -1036,31 +1000,25 @@ void test_is_zero_at_homogeneous() {
     
     typename PT::Is_zero_at_homogeneous is_zero_at_homogeneous;
     
-    Polynomial_d p(Coeff(-1), Coeff(0), Coeff(1));
+    Polynomial_d p(Coeff(-1), Coeff(0), Coeff(4));
+    std::cout << p << std::endl;
     
     std::vector< ICoeff > cv;
     for(int i = 0; i < PT::d-1; ++i)
         cv.push_back(ICoeff(0));
-        
-    for(int v = 1; v < 5; ++v) {
-        cv.push_back(ICoeff(0));
-        cv.push_back(ICoeff(v));
-        assert(!is_zero_at_homogeneous(p, cv.begin(), cv.end()));
-        
-        cv.pop_back();
-        cv.pop_back();
-        cv.push_back(ICoeff(v));
-        cv.push_back(ICoeff(v));
-        assert(is_zero_at_homogeneous(p, cv.begin(), cv.end()));
-        
-        cv.pop_back();
-        cv.pop_back();
-        cv.push_back(ICoeff(-v));
-        cv.push_back(ICoeff(v));
-        assert(is_zero_at_homogeneous(p, cv.begin(), cv.end()));
-        
-        cv.pop_back();
-        cv.pop_back();
+    
+    for(int v = 2; v < 5; ++v) {
+      cv.push_back(ICoeff(v));
+      cv.push_back(ICoeff(2*v));
+      assert(is_zero_at_homogeneous(p, cv.begin(), cv.end()));
+      cv.pop_back();
+      cv.pop_back();
+
+      cv.push_back(ICoeff(v));
+      cv.push_back(ICoeff(2*-v));
+      assert(is_zero_at_homogeneous(p, cv.begin(), cv.end()));
+      cv.pop_back();
+      cv.pop_back();
     }
     
     std::cerr << " ok" << std::endl;
@@ -1258,15 +1216,7 @@ void test_substitute(){
     assert(Innermost_coefficient(2) 
             == substitute(Polynomial_d(2),vec.begin(),vec.end()));
     assert(Innermost_coefficient(-2) 
-            == substitute(Polynomial_d(-2),vec.begin(),vec.end()));
-    
-    for(int i = 0; i< 5; i++){
-        Polynomial_d p = 
-            generate_sparse_random_polynomial<Polynomial_d>(3); 
-        assert(typename PT::Evaluate()(p,vec.begin(),vec.end()) 
-                == substitute(p,vec.begin(),vec.end()));
-    }
-    
+            == substitute(Polynomial_d(-2),vec.begin(),vec.end())); 
     
     for(int i = 0; i< 5; i++){
         typedef typename PT
@@ -1289,6 +1239,50 @@ void test_substitute(){
     std::cerr << " ok "<< std::endl;
 }
 
+// //       Substitute;
+template <class Polynomial_traits_d>
+void test_substitute_homogeneous(){
+    std::cerr << "start test_substitute_homogeneous "; std::cerr.flush();
+    CGAL_SNAP_CGALi_TRAITS_D(Polynomial_traits_d);
+    typename PT::Substitute_homogeneous substitute_homogeneous;
+    typedef typename PT::Innermost_coefficient Innermost_coefficient;
+    
+    
+    std::vector<Innermost_coefficient> vec;
+    for(int i = 0; i < PT::d; i++){
+        vec.push_back(Innermost_coefficient(i));
+    }
+    vec.push_back(Innermost_coefficient(3));
+    assert(Innermost_coefficient(0) 
+            == substitute_homogeneous(Polynomial_d(0),vec.begin(),vec.end()));
+    assert(Innermost_coefficient(1) 
+            == substitute_homogeneous(Polynomial_d(1),vec.begin(),vec.end()));
+    assert(Innermost_coefficient(2) 
+            == substitute_homogeneous(Polynomial_d(2),vec.begin(),vec.end()));
+    assert(Innermost_coefficient(-2) 
+            == substitute_homogeneous(Polynomial_d(-2),vec.begin(),vec.end())); 
+    
+    for(int i = 0; i< 5; i++){
+        typedef typename PT
+            :: template Rebind<Innermost_coefficient,5>::Other PT_5;
+        typedef typename PT_5::Polynomial_d Polynomial_5;
+        std::vector<Polynomial_5> vec1,vec2;
+        for(int j = 0; j < PT::d+1; j++){
+            vec1.push_back(
+                    generate_sparse_random_polynomial<Polynomial_5>(3));
+        }
+        vec2=vec1;
+        std::swap(vec2[0],vec2[PT::d-1]);
+        Polynomial_d p 
+            = generate_sparse_random_polynomial<Polynomial_d>(3); 
+        assert( substitute_homogeneous(p,vec1.begin(),vec1.end()) == 
+                substitute_homogeneous(typename PT::Swap()(p,0,PT::d-1),
+                        vec2.begin(),vec2.end()));
+    }   
+    
+    std::cerr << " ok "<< std::endl;
+}
+
 
 
 
@@ -1373,6 +1367,8 @@ struct Test_polynomial_traits_d<Polynomial_traits_d, CGAL::Null_tag > {
         test_canonicalize<Polynomial_traits_d>();
         //      Substitute;
         test_substitute<Polynomial_traits_d>();
+        //      Substitute_homogeneous;
+        test_substitute_homogeneous<Polynomial_traits_d>();
         //   private:
         //       Innermost_leading_coefficient;
     }