Integer Representation

doc/multiprecision.qbk

@@ -5998,7 +5998,19 @@ Open question - what should be the default - int32_t or int64_t?  (done 2012/09/
 [[Why not abstract out addition/multiplication algorithms?]
    [This was deemed not to be practical: these algorithms are intimately
    tied to the actual data representation used.]]
-]
+ [[How do I choose between Boost.Multiprecision cpp_bin_50 and cpp_dec_50?]
+   [Unless you have a specific reason to choose `cpp_dec_`, then the default choice should be `cpp_bin_`, 
+   for example using the convenience `typedefs` like `boost::multiprecision::cpp_bin_50` or `boost::multiprecision::cpp_bin_100`.
+
+  In general, both work well and give the same results and at roughly the same speed with `cpp_dec_50` sometimes faster.
+
+  `cpp_dec_` was developed first paving the way for `cpp_bin_`. 
+  `cpp_dec_` has several guard digits and is not rounded at all, using 'brute force' to get the promised number of decimal digits correct,
+  but making it difficult to reason about precision and computational uncertainty, for example see [*https://svn.boost.org/trac10/ticket/12133].
+  It also has a fast but imprecise division operator giving surprising results sometimes, see [*https://svn.boost.org/trac10/ticket/11178].
+
+  `cpp_bin_` is correctly/exactly rounded making it possible to reason about both the precision and rounding of the results.]]
+] [/variablelist]
 
 [endsect]
 

include/boost/multiprecision/detail/default_ops.hpp

@@ -824,6 +824,15 @@ inline int eval_get_sign(const T& val)
    return val.compare(static_cast<ui_type>(0));
 }
 
+template<class U, class T>
+inline void assign_components_imp(T& result, const long long& v1, const U& v2, const mpl::int_<number_kind_modulus>&)
+{
+    result.m_params = v2;
+    result = v1;
+    eval_multiply(result, v2.R2().backend());
+    //eval_multiply(result, v1);
+}
+
 template <class T, class V, class U>
 inline void assign_components_imp(T& result, const V& v1, const U& v2, const mpl::int_<number_kind_rational>&)
 {
@@ -837,7 +846,6 @@ template <class T, class V, class U, int N>
 inline void assign_components_imp(T& result, const V& v1, const U& v2, const mpl::int_<N>&)
 {
    typedef typename component_type<number<T> >::type component_number_type;
-
    component_number_type x(v1), y(v2);
    assign_components(result, x.backend(), y.backend());
 }

include/boost/multiprecision/detail/number_base.hpp

@@ -1546,7 +1546,8 @@ enum number_category_type
    number_kind_floating_point = 1,
    number_kind_rational = 2,
    number_kind_fixed_point = 3,
-   number_kind_complex = 4
+   number_kind_complex = 4,
+   number_kind_modulus = 5
 };
 
 template <class Num, bool, bool>

include/boost/multiprecision/inverse_euclid.hpp

@@ -0,0 +1,103 @@
+#ifndef CRYPTO3_INVERSE_EUCLID_HPP
+#define CRYPTO3_INVERSE_EUCLID_HPP
+
+#include <boost/multiprecision/cpp_int.hpp>
+
+namespace nil {
+    namespace crypto3 {
+        using namespace boost::multiprecision;
+
+        template<typename Backend, expression_template_option ExpressionTemplates>
+        inline Backend eval_inverse_euclid(const Backend &n, const Backend &mod) {
+
+            typedef typename default_ops::double_precision_type<Backend>::type double_type;
+            typedef typename boost::multiprecision::detail::canonical<unsigned int, double_type>::type ui_type;
+
+            using default_ops::eval_is_zero;
+            using default_ops::eval_lt;
+            using default_ops::eval_eq;
+            using default_ops::eval_add;
+            using default_ops::eval_subtract;
+            using default_ops::eval_get_sign;
+            using default_ops::eval_right_shift;
+            using default_ops::eval_integer_modulus;
+
+            if (eval_is_zero(mod)) {
+                BOOST_THROW_EXCEPTION(std::invalid_argument("inverse_mod: arguments must be non-zero"));
+            }
+            if (eval_get_sign(mod) || eval_get_sign(n)) {
+                BOOST_THROW_EXCEPTION(std::invalid_argument("inverse_mod: arguments must be non-negative"));
+            }
+
+            if (eval_is_zero(n) || (!eval_integer_modulus(n, 2) && !eval_integer_modulus(mod, 2))) {
+                return 0;
+            } // fast fail checks
+
+            number<Backend, ExpressionTemplates> u = mod, v = n;
+            number<Backend, ExpressionTemplates> a = 1, b = 0, c = 0, d = 1;
+
+            while (!eval_is_zero(u)) {
+                const ui_type u_zero_bits = lsb(u);
+                eval_right_shift(u, u_zero_bits);
+
+                for (std::size_t i = 0; i != u_zero_bits; ++i) {
+                    if (eval_integer_modulus(a, 2) || eval_integer_modulus(b, 2)) {
+                        eval_add(a, n);
+                        eval_subtract(b, mod);
+                    }
+                    eval_right_shift(a, 1);
+                    eval_right_shift(b, 1);
+                }
+
+                const ui_type v_zero_bits = lsb(v);
+                eval_right_shift(v, v_zero_bits);
+                for (size_t i = 0; i != v_zero_bits; ++i) {
+                    if (eval_integer_modulus(c, 2) || eval_integer_modulus(d, 2)) {
+                        eval_add(c, n);
+                        eval_subtract(d, mod);
+                    }
+                    eval_right_shift(c, 1);
+                    eval_right_shift(d, 1);
+                }
+
+                if (!eval_lt(u, v)) {
+                    eval_subtract(u, v);
+                    eval_subtract(a, c);
+                    eval_subtract(b, d);
+                } else {
+                    eval_subtract(v, u);
+                    eval_subtract(c, a);
+                    eval_subtract(d, b);
+                }
+            }
+
+            if (!eval_eq(v, 1)) {
+                return 0;
+            } // no modular inverse
+
+            while (eval_lt(d, 0)) {
+                eval_add(d, mod);
+            }
+            while (d >= mod) {
+                eval_subtract(d, mod);
+            }
+
+            return d;
+        }
+
+        /**
+         * Modular inversion using extended binary Euclidian algorithm
+         * @param x a positive integer
+         * @param modulus a positive integer
+         * @return y st (x*y) % modulus == 1 or 0 if no such value
+         * @note Non-const time
+         */
+        template<typename Backend, expression_template_option ExpressionTemplates>
+        inline number<Backend, ExpressionTemplates> inverse_euclid(const number<Backend, ExpressionTemplates> &n,
+                                                                   const number<Backend, ExpressionTemplates> &mod) {
+            return number<Backend, ExpressionTemplates>(eval_inverse_euclid(n.backend(), mod.backend()));
+        }
+    }
+}
+
+#endif //CRYPTO3_INVERSE_EUCLID_HPP

include/boost/multiprecision/jacobi.hpp

@@ -0,0 +1,82 @@
+#ifndef CRYPTO3_JACOBI_HPP
+#define CRYPTO3_JACOBI_HPP
+
+#include <boost/multiprecision/cpp_int.hpp>
+
+namespace nil {
+    namespace crypto3 {
+        using namespace boost::multiprecision;
+
+        template<typename Backend>
+        inline limb_type eval_jacobi(const Backend &a, const Backend &n) {
+            using default_ops::eval_is_zero;
+            using default_ops::eval_lt;
+            using default_ops::eval_lsb;
+            using default_ops::eval_gt;
+            using default_ops::eval_integer_modulus;
+            using default_ops::eval_modulus;
+            using default_ops::eval_right_shift;
+            using default_ops::eval_get_sign;
+
+            if (eval_get_sign(a) < 0) {
+                BOOST_THROW_EXCEPTION(std::invalid_argument("jacobi: first argument must be non-negative"));
+            }
+            if (!eval_integer_modulus(n, 2) || eval_lt(n, 2)) {
+                BOOST_THROW_EXCEPTION(std::invalid_argument("jacobi: second argument must be odd and > 1"));
+            }
+
+            Backend x = a, y = n;
+            limb_type J = 1;
+
+            while (eval_gt(y, 1)) {
+                eval_modulus(x, y);
+
+                Backend yd2 = y;
+                eval_divide(yd2, 2);
+
+                if (eval_gt(x, yd2)) {
+                    eval_subtract(x, y, x);
+                    if (eval_integer_modulus(y, 4) == 3) {
+                        J = -J;
+                    }
+                }
+                if (eval_is_zero(x)) {
+                    return 0;
+                }
+
+                size_t shifts = eval_lsb(x);
+                eval_right_shift(x, shifts);
+                if (eval_integer_modulus(shifts, 2)) {
+                    limb_type y_mod_8 = eval_integer_modulus(y, 8);
+                    if (y_mod_8 == 3 || y_mod_8 == 5) {
+                        J = -J;
+                    }
+                }
+
+                if (eval_integer_modulus(x, 4) == 3 && eval_integer_modulus(y, 4) == 3) {
+                    J = -J;
+                }
+
+                std::swap(x, y);
+            }
+            return J;
+        }
+
+        /**
+         * Compute the Jacobi symbol. If n is prime, this is equivalent
+         * to the Legendre symbol.
+         * @see http://mathworld.wolfram.com/JacobiSymbol.html
+         *
+         * @param a is a non-negative integer
+         * @param n is an odd integer > 1
+         * @return (n / m)
+         */
+        template<typename Backend, expression_template_option ExpressionTemplates>
+        inline typename std::enable_if<number_category<Backend>::value == number_kind_integer, limb_type>::type jacobi(
+                const number<Backend, ExpressionTemplates> &a, const number<Backend, ExpressionTemplates> &n) {
+            return number<Backend, ExpressionTemplates>(eval_jacobi(a.backend(), n.backend()));
+        }
+    }
+}
+
+#endif //CRYPTO3_JACOBI_HPP

include/boost/multiprecision/mask_bits.hpp

@@ -0,0 +1,37 @@
+#ifndef CRYPTO3_MASK_BITS_HPP
+#define CRYPTO3_MASK_BITS_HPP
+
+#include <boost/multiprecision/cpp_int.hpp>
+#include <boost/multiprecision/detail/number_base.hpp>
+#include <boost/multiprecision/number.hpp>
+
+namespace nil {
+    namespace crypto3 {
+
+        template<typename Backend, typename Integer>
+        void eval_mask_bits(Backend &val, Integer n) {
+            typedef typename boost::multiprecision::limb_type limb_type;
+
+            typedef typename boost::multiprecision::detail::canonical<unsigned, Backend>::type ui_type;
+            static const ui_type zero = 0u;
+
+            if (n == 0) {
+                val = zero;
+                return;
+            }
+
+            const size_t top_word = n / Backend::limb_bits;
+            const limb_type mask = (limb_type(1) << (n % Backend::limb_bits)) - 1;
+
+            if (top_word < val.size()) {
+                const size_t len = val.size() - (top_word + 1);
+                if (len > 0) {
+                    //clear_mem(&val.limbs()[top_word + 1], len); #TODO return this
+                }
+                val.limbs()[top_word] &= mask;
+            }
+        }
+    }
+}
+
+#endif //CRYPTO3_MASK_BITS_HPP