diff --git a/Algebraic_foundations/include/CGAL/extended_euclidean_algorithm.h b/Algebraic_foundations/include/CGAL/extended_euclidean_algorithm.h
index e17bd6a8e26..cfa53ce3899 100644
--- a/Algebraic_foundations/include/CGAL/extended_euclidean_algorithm.h
+++ b/Algebraic_foundations/include/CGAL/extended_euclidean_algorithm.h
@@ -23,46 +23,53 @@
 #define CGAL_EXTENDED_EUCLIDEAN_ALGORITHM_H 1
 
 #include <CGAL/basic.h>
+#include <vector>
 
 CGAL_BEGIN_NAMESPACE
 
-template< class NT > 
-NT extended_euclidean_algorithm(const NT& a_, const NT& b_, NT& u, NT& v){
+// EEA computing the normalized gcd
+// Modern Computer Algebra (Hardcover)
+// by Joachim von zur Gathen (Author), Jürgen Gerhard (Author) 
+// Publisher: Cambridge University Press; 2 edition (September 1, 2003)
+//  Language: English
+// ISBN-10: 0521826462
+// ISBN-13: 978-0521826464
+// pp.: 55
+ 
+template< class AS > 
+AS extended_euclidean_algorithm(const AS& f, const AS& g, AS& s_, AS& t_){
+    typename Algebraic_structure_traits<AS>::Integral_division idiv;
+    typename Algebraic_structure_traits<AS>::Div div;
+    typename Algebraic_structure_traits<AS>::Unit_part unit_part; 
 
-    typedef Algebraic_structure_traits<NT> AST;
-    typename AST::Div_mod div_mod;
-    typename AST::Unit_part unit_part;
-    typename AST::Integral_division idiv;
+    std::vector<AS> p,r,s,t,q; 
+    p.push_back(unit_part(f));   
+    r.push_back(idiv(f,p[0]));   
+    s.push_back(idiv(AS(1),p[0])); 
+    t.push_back(AS(0));  
+    q.push_back(AS(0));
     
-    NT unit_part_a(unit_part(a_));
-    NT unit_part_b(unit_part(b_));
+    p.push_back(unit_part(g));
+    r.push_back(idiv(g,p[1]));
+    s.push_back(AS(0));
+    t.push_back(idiv(AS(1),p[1]));  
     
-    NT a(idiv(a_,unit_part_a));
-    NT b(idiv(b_,unit_part_b));
-    
-    NT x(0),y(1),last_x(1),last_y(0);
-    NT temp, quotient; 
-  
-    //TODO: unroll to avoid swapping 
-    while (b != 0){
-        temp = b;
-        div_mod(a,b,quotient,b);
-        a = temp;
-     
-        temp = x;
-        x = last_x-quotient*x;
-        last_x = temp;
-      
-        temp = y;
-        y = last_y-quotient*y;
-        last_y = temp;
+    int i = 1; 
+    while(!is_zero(r[i])){
+        q.push_back(div(r[i-1],r[i]));
+        r.push_back(r[i-1]-q[i]*r[i]);
+        p.push_back(unit_part(r[i+1]));
+        r[i+1] = idiv(r[i+1],p[i+1]);
+        s.push_back(idiv(s[i-1]-q[i]*s[i],p[i+1]));
+        t.push_back(idiv(t[i-1]-q[i]*t[i],p[i+1])); 
+        i++;
     }
-    u = last_x * unit_part_a;
-    v = last_y * unit_part_b;
-
-    CGAL_precondition(unit_part(a) == NT(1));
-    CGAL_precondition(a == a_*u + b_*v);
-    return a; 
+    
+    s_=s[i-1];
+    t_=t[i-1];
+    AS h = r[i-1];
+    CGAL_precondition( h == f*s_ + g*t_);
+    return h; 
 }
 
 CGAL_END_NAMESPACE