[ci skip] Modular exponentiation, modular multiplicative inverse, extended Euclidean algorithm, discrete logarithm.

include/boost/integer/discrete_log.hpp

@@ -0,0 +1,140 @@
+/*
+ *  (C) Copyright Nick Thompson 2018.
+ *  Use, modification and distribution are subject to the
+ *  Boost Software License, Version 1.0. (See accompanying file
+ *  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+ *
+ *  Two methods of computing the discrete logarithm over the multiplicative group of integers mod p.
+ */
+
+#ifndef BOOST_INTEGER_DISCRETE_LOG_HPP
+#define BOOST_INTEGER_DISCRETE_LOG_HPP
+#include <limits>
+#include <unordered_map>
+#include <boost/optional.hpp>
+#include <boost/integer/floor_sqrt.hpp>
+#include <boost/integer/modular_multiplicative_inverse.hpp>
+#include <boost/integer/modular_exponentiation.hpp>
+#include <boost/integer/common_factor.hpp>
+
+namespace boost { namespace integer {
+
+// base^^x = a mod p <-> x = log_base(a) mod p
+template<class Z>
+boost::optional<Z> trial_multiplication_discrete_log(Z base, Z arg, Z p)
+{
+    using std::numeric_limits;
+    static_assert(numeric_limits<Z>::is_integer,
+                 "The discrete log works on integral types.\n");
+
+    if (base <= 1)
+    {
+        throw std::logic_error("The base must be > 1.\n");
+    }
+    if (p < 3)
+    {
+        throw std::logic_error("The modulus must be > 2.\n");
+    }
+    if (arg < 1)
+    {
+        throw std::logic_error("The argument must be > 0.\n");
+    }
+    if (base >= p || arg >= p)
+    {
+        throw std::logic_error("Error computing the discrete log: Are your arguments in the wrong order?\n");
+    }
+
+    if (arg == 1)
+    {
+        return 0;
+    }
+    Z s = 1;
+    for (Z i = 1; i < p; ++i)
+    {
+        s = (s * base) % p;
+        if (s == arg)
+        {
+            return i;
+        }
+    }
+    return {};
+}
+
+template<class Z>
+class baby_step_giant_step_discrete_log
+{
+public:
+    baby_step_giant_step_discrete_log(Z base, Z p) : m_p{p}
+    {
+        using std::numeric_limits;
+        static_assert(numeric_limits<Z>::is_integer,
+                      "The baby_step_giant_step discrete log works on integral types.\n");
+
+        if (base <= 1)
+        {
+            throw std::logic_error("The base must be > 1.\n");
+        }
+        if (p < 3)
+        {
+            throw std::logic_error("The modulus must be > 2.\n");
+        }
+        if (base >= p)
+        {
+            throw std::logic_error("Error computing the discrete log: Are your arguments in the wrong order?\n");
+        }
+        m_root_p = floor_sqrt(p);
+        if (m_root_p*m_root_p != p)
+        {
+            m_root_p += 1;
+        }
+
+        auto x = modular_multiplicative_inverse(base, p);
+        if (!x)
+        {
+            throw std::logic_error("The gcd of the b and the modulus is > 1, hence the discrete log is not guaranteed to exist. If you don't require an existence proof, use trial multiplication.\n");
+        }
+        m_inv_base_pow_m = modular_exponentiation(x.value(), m_root_p, p);
+
+        m_lookup_table.reserve(m_root_p);
+        // Now the expensive part:
+        Z k = 1;
+        for (Z j = 0; j < m_root_p; ++j)
+        {
+            m_lookup_table.emplace(k, j);
+            k = k*base % p;
+        }
+
+    }
+
+    Z operator()(Z arg) const
+    {
+        Z ami = m_inv_base_pow_m;
+        Z k = arg % m_p;
+        if(k == 0)
+        {
+            throw std::domain_error("Cannot take the logarithm of a number divisible by the modulus.\n");
+        }
+        for (Z i = 0; i < m_root_p; ++i)
+        {
+            auto it = m_lookup_table.find(k);
+            if (it != m_lookup_table.end())
+            {
+                return (i*m_root_p + it->second) % m_p;
+            }
+            ami = (ami*m_inv_base_pow_m) % m_p;
+            k = k * ami % m_p;
+        }
+        // never should get here . . .
+        return -1;
+    }
+
+private:
+    Z m_p;
+    Z m_root_p;
+    Z m_inv_base_pow_m;
+    std::unordered_map<Z, Z> m_lookup_table;
+};
+
+
+}}
+#endif

include/boost/integer/extended_euclidean.hpp

@@ -0,0 +1,67 @@
+/*
+ *  (C) Copyright Nick Thompson 2018.
+ *  Use, modification and distribution are subject to the
+ *  Boost Software License, Version 1.0. (See accompanying file
+ *  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+#ifndef BOOST_INTEGER_EXTENDED_EUCLIDEAN_HPP
+#define BOOST_INTEGER_EXTENDED_EUCLIDEAN_HPP
+#include <tuple>
+#include <limits>
+
+namespace boost { namespace integer {
+
+// From "The Joy of Factoring", Algorithm 2.7.
+// Should the tuple be a named tuple? Is that possible?
+// Solves mx + ny = gcd(m,n). Returns tuple with (gcd(m,n), x, y).
+template<class Z>
+std::tuple<Z, Z, Z> extended_euclidean(Z m, Z n)
+{
+    using std::numeric_limits;
+    static_assert(numeric_limits<Z>::is_integer,
+                  "The extended Euclidean algorithm works on integral types.\n");
+
+    static_assert(numeric_limits<Z>::is_signed,
+                  "The extended Euclidean algorithm only works on signed integer types.\n");
+
+    if (m < 1 || n < 1)
+    {
+        throw 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);
+    }
+    return std::make_tuple(u0, u1, u2);
+}
+
+}}
+#endif

include/boost/integer/floor_sqrt.hpp

@@ -0,0 +1,34 @@
+/*
+ *  (C) Copyright Nick Thompson 2017.
+ *  Use, modification and distribution are subject to the
+ *  Boost Software License, Version 1.0. (See accompanying file
+ *  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+ *
+ *  The integer floor_sqrt doesn't lose precision like a cast does.
+ *  Based on Algorithm 5.9 of "The Joy of Factoring".
+ */
+
+
+#ifndef BOOST_INTEGER_FLOOR_SQRT_HPP
+#define BOOST_INTEGER_FLOOR_SQRT_HPP
+#include <limits>
+
+namespace boost { namespace integer {
+
+template<class Z>
+Z floor_sqrt(Z N)
+{
+    static_assert(std::numeric_limits<Z>::is_integer,
+                  "The floor_sqrt function is for taking square roots of integers.\n");
+
+     Z x = N;
+     Z y = x/2 + (x&1);
+     while (y < x) {
+        x = y;
+        y = (x + N / x)/2;
+     }
+     return x;
+}
+}}
+
+#endif

include/boost/integer/modular_exponentiation.hpp

@@ -0,0 +1,39 @@
+/*
+ *  (C) Copyright Nick Thompson 2018.
+ *  Use, modification and distribution are subject to the
+ *  Boost Software License, Version 1.0. (See accompanying file
+ *  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+#ifndef BOOST_INTEGER_MODULAR_EXPONENTIATION_HPP
+#define BOOST_INTEGER_MODULAR_EXPONENTIATION_HPP
+#include <limits>
+
+namespace boost { namespace integer {
+
+template<class Z>
+Z modular_exponentiation(Z base, Z exponent, Z modulus)
+{
+    using std::numeric_limits;
+    static_assert(numeric_limits<Z>::is_integer,
+                  "Modular exponentiation works on integral types.\n");
+    Z result = 1;
+    if (exponent < 0 || modulus < 0)
+    {
+        throw std::domain_error("Both the exponent and the modulus must be > 0.\n");
+    }
+
+    while (exponent > 0)
+    {
+        if (exponent & 1)
+        {
+            result = (result*base) % modulus;
+        }
+        base = (base*base) % modulus;
+        exponent >>= 1;
+    }
+    return result;
+}
+
+
+}}
+#endif

include/boost/integer/modular_multiplicative_inverse.hpp

@@ -0,0 +1,50 @@
+/*
+ *  (C) Copyright Nick Thompson 2018.
+ *  Use, modification and distribution are subject to the
+ *  Boost Software License, Version 1.0. (See accompanying file
+ *  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+ */
+#ifndef BOOST_INTEGER_MODULAR_MULTIPLICATIVE_INVERSE_HPP
+#define BOOST_INTEGER_MODULAR_MULTIPLICATIVE_INVERSE_HPP
+#include <limits>
+#include <boost/optional.hpp>
+#include <boost/integer/extended_euclidean.hpp>
+
+namespace boost { namespace integer {
+
+// From "The Joy of Factoring", Algorithm 2.7.
+template<class Z>
+boost::optional<Z> modular_multiplicative_inverse(Z a, Z modulus)
+{
+    using std::numeric_limits;
+    static_assert(numeric_limits<Z>::is_integer,
+                  "The modular multiplicative inverse works on integral types.\n");
+    if (modulus < 2)
+    {
+        throw std::domain_error("Modulus must be > 1.\n");
+    }
+    // make sure a < modulus:
+    a = a % modulus;
+    if (a == 0)
+    {
+        // a doesn't have a modular multiplicative inverse:
+        return {};
+    }
+    auto u = extended_euclidean(a, modulus);
+    Z gcd = std::get<0>(u);
+    if (gcd > 1)
+    {
+        return {};
+    }
+    Z x = std::get<1>(u);
+    // x might not be in the range 0 < x < m, let's fix that:
+    x = x % modulus;
+    while (x <= 0)
+    {
+        x += modulus;
+    }
+    return x;
+}
+
+}}
+#endif