diff --git a/Polynomial/include/CGAL/Polynomial/fwd.h b/Polynomial/include/CGAL/Polynomial/fwd.h
index 31ea2545dc2..fea83b67710 100644
--- a/Polynomial/include/CGAL/Polynomial/fwd.h
+++ b/Polynomial/include/CGAL/Polynomial/fwd.h
@@ -97,6 +97,15 @@ template <typename OutputIterator, typename NT > inline
 template<typename OutputIterator, typename NT> OutputIterator
     sturm_habicht_sequence(CGAL::Polynomial<NT> P, OutputIterator out);
 
+  template<typename OutputIterator1,
+    typename OutputIterator2,
+    typename OutputIterator3,
+    typename NT> OutputIterator1
+    sturm_habicht_sequence_with_cofactors(CGAL::Polynomial<NT> P, 
+                                          OutputIterator1 out_stha,
+                                          OutputIterator2 out_f,
+                                          OutputIterator3 out_fx);
+
 
 } // namespace CGALi
 
diff --git a/Polynomial/include/CGAL/Polynomial_traits_d.h b/Polynomial/include/CGAL/Polynomial_traits_d.h
index 9d8e5eef7a9..36a88e9149b 100644
--- a/Polynomial/include/CGAL/Polynomial_traits_d.h
+++ b/Polynomial/include/CGAL/Polynomial_traits_d.h
@@ -1381,8 +1381,28 @@ public:
                                                            out);
         }  
     };
+
+    //       Sturm-Habicht sequence with cofactors
+    struct Sturm_habicht_sequence_with_cofactors {
+        template<typename OutputIterator1,
+                 typename OutputIterator2,
+                 typename OutputIterator3>
+        OutputIterator1 operator()(
+                const Polynomial_d& p, 
+                OutputIterator1 out_stha,
+                OutputIterator2 out_f,
+                OutputIterator3 out_fx,
+                int i = (d-1) ) const {
+            if(i == (d-1) )
+                return CGAL::CGALi::sturm_habicht_sequence_with_cofactors
+                    (p,out_stha,out_f,out_fx);
+            else
+                return CGAL::CGALi::sturm_habicht_sequence_with_cofactors
+                    (Move()(p,i),out_stha,out_f,out_fx);
+        }  
+    };
     
-    //       Sturm-Habicht sequence
+    //       Principal Sturm-Habicht sequence
     struct Principal_sturm_habicht_sequence {
         
         template<typename OutputIterator>
diff --git a/Polynomial/test/Polynomial/sturm_habicht_sequence.cpp b/Polynomial/test/Polynomial/sturm_habicht_sequence.cpp
index 7f05e5cff9c..180f4aa949d 100644
--- a/Polynomial/test/Polynomial/sturm_habicht_sequence.cpp
+++ b/Polynomial/test/Polynomial/sturm_habicht_sequence.cpp
@@ -207,6 +207,23 @@ template<typename ArithmeticKernel> void test_routine() {
       assert(stha[i]==co_f[i]*f + co_fx[i]*fx);
     }
 
+    stha.clear();
+    co_f.clear();
+    co_fx.clear();
+    typename CGAL::Polynomial_traits_d<PNT_1>
+        ::Sturm_habicht_sequence_with_cofactors()(f,
+                                                  std::back_inserter(stha),
+                                                  std::back_inserter(co_f),
+                                                  std::back_inserter(co_fx)
+        );
+    
+    n = static_cast<int>(stha.size());
+    assert(n == static_cast<int>(stha_check.size()));
+    for(int i = 0; i < n; i++) {
+      assert(stha[i]==stha_check[i]);
+      assert(stha[i]==co_f[i]*f + co_fx[i]*fx);
+    }
+
   }
 
   { // Cofactors, defective case