From 54d0e4c63e74951e80a23db478214d4681704d6c Mon Sep 17 00:00:00 2001
From: Nick Thompson <nathompson7@protonmail.com>
Date: Thu, 25 Oct 2018 14:28:39 -0600
Subject: [PATCH] [ci skip] Trade out algorithm from 'The Joy of Factoring' to
 Wikipedia's version which reduces the number of required temporaries. In
 fact, the speedup is not large, but the code is more compact, and for larger
 types, the difference becomes more noticeable.

---
 include/boost/integer/extended_euclidean.hpp | 59 +++++++-------------
 1 file changed, 20 insertions(+), 39 deletions(-)

diff --git a/include/boost/integer/extended_euclidean.hpp b/include/boost/integer/extended_euclidean.hpp
index 64d10d9..93c8c54 100644
--- a/include/boost/integer/extended_euclidean.hpp
+++ b/include/boost/integer/extended_euclidean.hpp
@@ -13,9 +13,8 @@
 
 namespace boost { namespace integer {
 
-// From "The Joy of Factoring", Algorithm 2.7.
+// From "The Joy of Factoring", Algorithm 2.7, with a small optimization to remove tmps from Wikipedia.
 // Solves mx + ny = gcd(m,n). Returns tuple with (gcd(m,n), x, y).
-// Is this the natural ordering?, or must people simply have to read the docs?
 
 template<class Z>
 struct euclidean_result_t {
@@ -38,45 +37,27 @@ euclidean_result_t<Z> extended_euclidean(Z m, Z n)
     {
         BOOST_THROW_EXCEPTION(std::domain_error("Arguments must be strictly positive.\n"));
     }
-    bool swapped = false;
-    if (m < n)
-    {
-        swapped = true;
-        std::swap(m, n);
-    }
-    Z u0 = m;
-    Z u1 = 1;
-    Z u2 = 0;
-    Z v0 = n;
-    Z v1 = 0;
-    Z v2 = 1;
-    Z w0;
-    Z w1;
-    Z w2;
-    while(v0 > 0)
-    {
-        Z q = u0/v0;
-        w0 = u0 - q*v0;
-        w1 = u1 - q*v1;
-        w2 = u2 - q*v2;
-        u0 = v0;
-        u1 = v1;
-        u2 = v2;
-        v0 = w0;
-        v1 = w1;
-        v2 = w2;
-    }
-    if (swapped)
-    {
-        std::swap(u1, u2);
-        BOOST_ASSERT(u2*m+u1*n==u0);
-    }
-    else
-    {
-        BOOST_ASSERT(u1*m+u2*n==u0);
+
+    Z s = 0;
+    Z old_s = 1;
+    Z r = n;
+    Z old_r = m;
+
+    while (r != 0) {
+      Z q = old_r/r;
+      Z tmp = r;
+      r = old_r - q*tmp;
+      old_r = tmp;
+
+      tmp = s;
+      s = old_s - q*tmp;
+      old_s = tmp;
     }
 
-    return {u0, u1, u2};
+    Z y = (old_r - old_s*m)/n;
+
+    BOOST_ASSERT(old_s*m+y*n==old_r);
+    return {old_r, old_s, y};
 }
 
 }}