[ci skip] Use less verbose naming. Add asserts as verfication of algorithms is a negligible fraction of total runtime. Use boost::multiprecision::powm and boost::multiprecision::sqrt rather than one-offs.

include/boost/integer/discrete_log.hpp

@@ -12,10 +12,10 @@
 #include <limits>
 #include <unordered_map>
 #include <boost/optional.hpp>
-#include <boost/integer/floor_sqrt.hpp>
+#include <boost/format.hpp>
-#include <boost/integer/modular_multiplicative_inverse.hpp>
+#include <boost/multiprecision/integer.hpp>
-#include <boost/integer/modular_exponentiation.hpp>
+#include <boost/integer/common_factor_rt.hpp>
-#include <boost/integer/common_factor.hpp>
+#include <boost/integer/mod_inverse.hpp>
 
 namespace boost { namespace integer {
 
@@ -29,19 +29,28 @@ boost::optional<Z> trial_multiplication_discrete_log(Z base, Z arg, Z p)
 
     if (base <= 1)
     {
-        throw std::logic_error("The base must be > 1.\n");
+        throw std::domain_error("The base must be > 1.\n");
     }
     if (p < 3)
     {
-        throw std::logic_error("The modulus must be > 2.\n");
+        throw std::domain_error("The modulus must be > 2.\n");
     }
     if (arg < 1)
     {
-        throw std::logic_error("The argument must be > 0.\n");
+        throw std::domain_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 (base >= p)
+        {
+            auto e = boost::format("Error computing the discrete log: The base %1% is greater than the modulus %2%. Are the arguments in the wrong order?") % base % p;
+            throw std::domain_error(e.str());
+        }
+        if (arg >= p)
+        {
+            auto e = boost::format("Error computing the discrete log: The argument %1% is greater than the modulus %2%. Are the arguments in the wrong order?") % arg % p;
+            throw std::domain_error(e.str());
+        }
     }
 
     if (arg == 1)
@@ -54,6 +63,8 @@ boost::optional<Z> trial_multiplication_discrete_log(Z base, Z arg, Z p)
         s = (s * base) % p;
         if (s == arg)
         {
+            // Maybe a bit trivial assertion. But still a negligible fraction of the total compute time.
+            BOOST_ASSERT(arg == boost::multiprecision::powm(base, i, p));
             return i;
         }
     }
@@ -61,14 +72,14 @@ boost::optional<Z> trial_multiplication_discrete_log(Z base, Z arg, Z p)
 }
 
 template<class Z>
-class baby_step_giant_step_discrete_log
+class bsgs_discrete_log
 {
 public:
-    baby_step_giant_step_discrete_log(Z base, Z p) : m_p{p}
+    bsgs_discrete_log(Z base, Z p) : m_p{p}, m_base{base}
     {
         using std::numeric_limits;
         static_assert(numeric_limits<Z>::is_integer,
-                      "The baby_step_giant_step discrete log works on integral types.\n");
+                      "The baby-step, giant-step discrete log works on integral types.\n");
 
         if (base <= 1)
         {
@@ -82,18 +93,20 @@ public:
         {
             throw std::logic_error("Error computing the discrete log: Are your arguments in the wrong order?\n");
         }
-        m_root_p = floor_sqrt(p);
+        m_root_p = boost::multiprecision::sqrt(p);
         if (m_root_p*m_root_p != p)
         {
             m_root_p += 1;
         }
 
-        auto x = modular_multiplicative_inverse(base, p);
+        auto x = mod_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");
+            auto d = boost::integer::gcd(base, p);
+            auto e = boost::format("The gcd of the base %1% and the modulus %2% is %3% != 1, hence the discrete log is not guaranteed to exist, which breaks the baby-step giant step algorithm. If you don't require existence proof for all inputs, use trial multiplication.\n") % base % p % d;
+            throw std::logic_error(e.str());
         }
-        m_inv_base_pow_m = modular_exponentiation(x.value(), m_root_p, p);
+        m_inv_base_pow_m = boost::multiprecision::powm(x.value(), m_root_p, p);
 
         m_lookup_table.reserve(m_root_p);
         // Now the expensive part:
@@ -119,17 +132,24 @@ public:
             auto it = m_lookup_table.find(k);
             if (it != m_lookup_table.end())
             {
-                return (i*m_root_p + it->second) % m_p;
+                Z log_b_arg = (i*m_root_p + it->second) % m_p;
+                // This computation of the modular exponentiation is laughably quick relative to computing the discrete log.
+                // Why not put an assert here for our peace of mind?
+                BOOST_ASSERT(arg == boost::multiprecision::powm(m_base, log_b_arg, m_p));
+                return log_b_arg;
             }
             ami = (ami*m_inv_base_pow_m) % m_p;
             k = k * ami % m_p;
         }
         // never should get here . . .
+        BOOST_ASSERT(false);
+        // Suppress compiler warnings.
         return -1;
     }
 
 private:
     Z m_p;
+    Z m_base;
     Z m_root_p;
     Z m_inv_base_pow_m;
     std::unordered_map<Z, Z> m_lookup_table;

include/boost/integer/extended_euclidean.hpp

@@ -59,7 +59,13 @@ std::tuple<Z, Z, Z> extended_euclidean(Z m, Z n)
     if (swapped)
     {
         std::swap(u1, u2);
+        BOOST_ASSERT(u2*m+u1*n==u0);
     }
+    else
+    {
+        BOOST_ASSERT(u1*m+u2*n==u0);
+    }
+
     return std::make_tuple(u0, u1, u2);
 }
 

include/boost/integer/floor_sqrt.hpp

@@ -1,34 +0,0 @@
-/*
- *  (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/mod_inverse.hpp

@@ -13,8 +13,13 @@
 namespace boost { namespace integer {
 
 // From "The Joy of Factoring", Algorithm 2.7.
+// The name is a bit verbose. Here's some others names I've found for this function:
+// PowerMod[a, -1, m] (Mathematica)
+// mpz_invert (gmplib)
+// modinv (some dude on stackoverflow)
+// Would modular_inverse be sometimes mistaken as the modular *additive* inverse?
 template<class Z>
-boost::optional<Z> modular_multiplicative_inverse(Z a, Z modulus)
+boost::optional<Z> mod_inverse(Z a, Z modulus)
 {
     using std::numeric_limits;
     static_assert(numeric_limits<Z>::is_integer,
@@ -37,12 +42,13 @@ boost::optional<Z> modular_multiplicative_inverse(Z a, Z modulus)
         return {};
     }
     Z x = std::get<1>(u);
-    // x might not be in the range 0 < x < m, let's fix that:
     x = x % modulus;
+    // x might not be in the range 0 < x < m, let's fix that:
     while (x <= 0)
     {
         x += modulus;
     }
+    BOOST_ASSERT(x*a % modulus == 1);
     return x;
 }
 

include/boost/integer/modular_exponentiation.hpp

@@ -1,39 +0,0 @@
-/*
- *  (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

test/Jamfile.v2

@@ -18,8 +18,7 @@ test-suite integer
         [ run static_min_max_test.cpp ]
         [ run discrete_log_test.cpp ]
         [ run extended_euclidean_test.cpp ]
-        [ run modular_exponentiation_test.cpp ]
+        [ run mod_inverse_test.cpp ]
-        [ run modular_multiplicative_inverse_test.cpp ]
         [ compile integer_traits_include_test.cpp ]
         [ compile integer_include_test.cpp ]
         [ compile integer_mask_include_test.cpp ]

test/discrete_log_test.cpp

@@ -8,10 +8,11 @@
 #define BOOST_TEST_MODULE discrete_log_test
 #include <boost/test/included/unit_test.hpp>
 #include <boost/integer/discrete_log.hpp>
+#include <boost/math/special_functions/prime.hpp>
 
 
 using boost::integer::trial_multiplication_discrete_log;
-using boost::integer::baby_step_giant_step_discrete_log;
+using boost::integer::bsgs_discrete_log;
 
 template<class Z>
 void test_trial_multiplication_discrete_log()
@@ -58,13 +59,52 @@ void test_trial_multiplication_discrete_log()
 template<class Z>
 void test_bsgs_discrete_log()
 {
-    baby_step_giant_step_discrete_log<Z> dl(7, 41);
+    bsgs_discrete_log<Z> dl_7(7, 41);
-    BOOST_CHECK_EQUAL(dl(7), 1);
+    BOOST_CHECK_EQUAL(dl_7(7), 1);
-    BOOST_CHECK_EQUAL(dl(8), 2);
+    BOOST_CHECK_EQUAL(dl_7(8), 2);
-    BOOST_CHECK_EQUAL(dl(15), 3);
+    BOOST_CHECK_EQUAL(dl_7(15), 3);
-    BOOST_CHECK_EQUAL(dl(23), 4);
+    BOOST_CHECK_EQUAL(dl_7(23), 4);
-    BOOST_CHECK_EQUAL(dl(38), 5);
+    BOOST_CHECK_EQUAL(dl_7(38), 5);
-    BOOST_CHECK_EQUAL(dl(20), 6);
+    BOOST_CHECK_EQUAL(dl_7(20), 6);
+}
+
+template<class Z>
+void test_trial_multiplication_with_prime_base()
+{
+    for (Z i = 0; i < boost::math::max_prime; ++i)
+    {
+        Z p = boost::math::prime(i);
+        for (Z j = 2; j < p; ++j)
+        {
+            bsgs_discrete_log<Z> dl_j(j, p);
+            for (Z k = 1; k < p; ++k)
+            {
+                boost::optional<Z> dl = trial_multiplication_discrete_log(j, k, p);
+                // It is guaranteed to exist with the modulus is prime:
+                BOOST_ASSERT(dl);
+                BOOST_CHECK_EQUAL(k, boost::multiprecision::powm(j, dl.value(), p));
+            }
+        }
+    }
+}
+
+
+template<class Z>
+void test_bsgs_with_prime_base()
+{
+    for (Z i = 0; i < boost::math::max_prime; ++i)
+    {
+        Z p = boost::math::prime(i);
+        for (Z j = 2; j < p; ++j)
+        {
+            bsgs_discrete_log<Z> dl_j(j, p);
+            for (Z k = 1; k < p; ++k)
+            {
+                Z dl = dl_j(k);
+                BOOST_CHECK_EQUAL(k, boost::multiprecision::powm(j, dl, p));
+            }
+        }
+    }
 }
 
 
@@ -72,4 +112,6 @@ BOOST_AUTO_TEST_CASE(discrete_log_test)
 {
     test_trial_multiplication_discrete_log<size_t>();
     test_bsgs_discrete_log<int>();
+    test_trial_multiplication_with_prime_base<long long>();
+    test_bsgs_with_prime_base<long long>();
 }

test/extended_euclidean_test.cpp

@@ -17,7 +17,7 @@ using boost::integer::gcd;
 template<class Z>
 void test_extended_euclidean()
 {
-    Z max_arg = 500;
+    Z max_arg = 1000;
     for (Z m = 1; m < max_arg; ++m)
     {
         for (Z n = 1; n < max_arg; ++n)
@@ -36,6 +36,6 @@ BOOST_AUTO_TEST_CASE(extended_euclidean_test)
 {
     test_extended_euclidean<int>();
     test_extended_euclidean<long>();
-    test_extended_euclidean<size_t>();
+    test_extended_euclidean<long long>();
     test_extended_euclidean<int128_t>();
 }

test/mod_inverse_test.cpp

@@ -8,14 +8,14 @@
 #include <boost/test/included/unit_test.hpp>
 #include <boost/multiprecision/cpp_int.hpp>
 #include <boost/integer/common_factor.hpp>
-#include <boost/integer/modular_multiplicative_inverse.hpp>
+#include <boost/integer/mod_inverse.hpp>
 
 using boost::multiprecision::int128_t;
-using boost::integer::modular_multiplicative_inverse;
+using boost::integer::mod_inverse;
 using boost::integer::gcd;
 
 template<class Z>
-void test_modular_multiplicative_inverse()
+void test_mod_inverse()
 {
     Z max_arg = 1000;
     for (Z modulus = 2; modulus < max_arg; ++modulus)
@@ -23,7 +23,7 @@ void test_modular_multiplicative_inverse()
         for (Z a = 1; a < max_arg; ++a)
         {
             Z gcdam = gcd(a, modulus);
-            boost::optional<Z> inv_a = modular_multiplicative_inverse(a, modulus);
+            boost::optional<Z> inv_a = mod_inverse(a, modulus);
             // Should fail if gcd(a, mod) != 1:
             if (gcdam > 1)
             {
@@ -41,8 +41,8 @@ void test_modular_multiplicative_inverse()
 
 BOOST_AUTO_TEST_CASE(extended_euclidean_test)
 {
-    test_modular_multiplicative_inverse<int>();
+    test_mod_inverse<int>();
-    test_modular_multiplicative_inverse<long>();
+    test_mod_inverse<long>();
-    test_modular_multiplicative_inverse<long long>();
+    test_mod_inverse<long long>();
-    test_modular_multiplicative_inverse<int128_t>();
+    test_mod_inverse<int128_t>();
 }

test/modular_exponentiation_test.cpp

@@ -1,38 +0,0 @@
-/*
- *  (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)
- *
- */
-
-#define BOOST_TEST_MODULE modular_exponentiation_test
-#include <boost/test/included/unit_test.hpp>
-#include <boost/multiprecision/cpp_int.hpp>
-#include <boost/integer/modular_exponentiation.hpp>
-
-using boost::multiprecision::int128_t;
-using boost::integer::modular_exponentiation;
-
-template<class Z>
-void test_modular_exponentiation()
-{
-    Z base = 7;
-    Z modulus = 51;
-    Z expected = 1;
-    for (Z exponent = 0; exponent < 10000; ++exponent)
-    {
-        Z x = modular_exponentiation<Z>(base, exponent, modulus);
-        BOOST_CHECK_EQUAL(expected, x);
-        expected = (expected*base) % modulus;
-    }
-}
-
-BOOST_AUTO_TEST_CASE(modular_exponentiation_test)
-{
-    test_modular_exponentiation<int>();
-    test_modular_exponentiation<unsigned>();
-    test_modular_exponentiation<short>();
-    test_modular_exponentiation<size_t>();
-    test_modular_exponentiation<int128_t>();
-}