improve MontgomeryForm pow(), add 128 bit lowuop REDC (not enabled yet), add and improve functions to experimental montgomery two_pow

Author: hurchalla

montgomery_arithmetic/include/hurchalla/montgomery_arithmetic/MontgomeryForm.h

@@ -43,6 +43,7 @@ class MontgomeryForm final {
     static_assert(ut_numeric_limits<T>::is_integer, "");
     static_assert(ut_numeric_limits<T>::digits <=
                   ut_numeric_limits<typename MontyType::uint_type>::digits, "");
+    using SV = typename MontyType::squaringvalue_type;
     using RU = typename MontyType::uint_type;
 public:
     using IntegerType = T;
@@ -477,6 +478,7 @@ class MontgomeryForm final {
 
     // Calculates and returns the modular exponentiation of the montgomery value
     // 'base' to the power of (the type T variable) 'exponent'.
+    // Performance note: if your base is 2, two_pow() is much more efficient.
     HURCHALLA_FORCE_INLINE
     MontgomeryValue pow(MontgomeryValue base, T exponent) const
     {

montgomery_arithmetic/include/hurchalla/montgomery_arithmetic/detail/ImplMontgomeryForm.contents

@@ -38,6 +38,7 @@ public:
     using MontgomeryValue = typename MontyType::montvalue_type;
     using CanonicalValue = typename MontyType::canonvalue_type;
     using FusingValue = typename MontyType::fusingvalue_type;
+    using SquaringValue = typename MontyType::squaringvalue_type;
 
     explicit ImplMontgomeryForm(T modulus) : impl(static_cast<U>(modulus)) {}
 
@@ -247,6 +248,17 @@ public:
         return impl.convertIn(a, PTAG());
     }
 
+    HURCHALLA_IMF_MAYBE_FORCE_INLINE CanonicalValue getMontvalueR() const
+    {
+        return impl.getMontvalueR();
+    }
+
+    template <class PTAG> HURCHALLA_IMF_MAYBE_FORCE_INLINE
+    MontgomeryValue twoPowLimited_times_x(size_t exponent, CanonicalValue x) const
+    {
+        return impl.twoPowLimited_times_x(exponent, x, PTAG());
+    }
+
     HURCHALLA_IMF_MAYBE_FORCE_INLINE U getMagicValue() const
     {
         return impl.getMagicValue();
@@ -268,3 +280,28 @@ public:
     {
         return impl.RTimesTwoPowLimited(exponent, magicValue, PTAG());
     }
+
+
+    HURCHALLA_IMF_MAYBE_FORCE_INLINE
+    SquaringValue getSquaringValue(MontgomeryValue x) const
+    {
+        return impl.getSquaringValue(x);
+    }
+
+    HURCHALLA_IMF_MAYBE_FORCE_INLINE
+    SquaringValue squareSV(SquaringValue sv) const
+    {
+        return impl.squareSV(sv);
+    }
+
+    HURCHALLA_IMF_MAYBE_FORCE_INLINE
+    MontgomeryValue squareToMontgomeryValue(SquaringValue sv) const
+    {
+        return impl.squareToMontgomeryValue(sv);
+    }
+
+    HURCHALLA_IMF_MAYBE_FORCE_INLINE
+    MontgomeryValue getMontgomeryValue(SquaringValue sv) const
+    {
+        return impl.getMontgomeryValue(sv);
+    }

montgomery_arithmetic/include/hurchalla/montgomery_arithmetic/detail/MontgomeryFormExtensions.h

@@ -30,14 +30,44 @@ struct MontgomeryFormExtensions final {
     static_assert(ut_numeric_limits<RU>::is_integer, "");
     static_assert(!(ut_numeric_limits<RU>::is_signed), "");
 
+    using CanonicalValue = typename MF::CanonicalValue;
     using MontgomeryValue = typename MF::MontgomeryValue;
+    using SquaringValue = typename MF::SV;
 
     HURCHALLA_FORCE_INLINE
     static MontgomeryValue convertInExtended(const MF& mf, RU a)
     {
         return mf.impl.template convertInExtended<PTAG>(a);
     }
 
+    // note: montvalueR is the Montgomery representation of R.
+    //       In normal integer form it is literally R squared mod N.
+    HURCHALLA_FORCE_INLINE
+    static CanonicalValue getMontvalueR(const MF& mf)
+    {
+        return mf.impl.getMontvalueR();
+    }
+
+    // this first shifts x by exponent, which is equivalent to
+    // multiplying x by 2^exponent, and then it completes the
+    // mont mul as usual by calling REDC.
+    // -- IMPORTANT NOTE -- because (2^exponent) is an integer domain
+    // value rather than a montgomery domain value, the returned
+    // result viewed as an integer value is
+    // REDC((x_int * R) * (2^exponent)) == (x_int * (2^exponent) * R) * R^(-1)
+    // To counteract the inverse R factor, so that you get what most likely
+    // you wanted, being just plain   (x_int * (2^exponent) * R),
+    // you need to ensure that x has an extra factor of R built into it it,
+    // rather than just the normal single factor of x_int * R.  To build an
+    // extra factor of R into x, you first get  montR = getMontvalueR(mf),
+    // and then you do a normal montgomery multiply of x and montR.
+    HURCHALLA_FORCE_INLINE
+    static MontgomeryValue twoPowLimited_times_x(const MF& mf, size_t exponent, CanonicalValue x)
+    {
+        HPBC_CLOCKWORK_PRECONDITION(exponent < ut_numeric_limits<RU>::digits);
+        return mf.impl.template twoPowLimited_times_x<PTAG>(exponent, x);
+    }
+
     // note: magicValue is R cubed mod N  (in normal integer form)
     HURCHALLA_FORCE_INLINE
     static RU getMagicValue(const MF& mf)
@@ -76,6 +106,32 @@ struct MontgomeryFormExtensions final {
         HPBC_CLOCKWORK_PRECONDITION(exponent < ut_numeric_limits<RU>::digits);
         return mf.impl.template RTimesTwoPowLimited<PTAG>(exponent, magicValue);
     }
+
+
+    HURCHALLA_FORCE_INLINE
+    static SquaringValue getSquaringValue(const MF& mf, MontgomeryValue x)
+    {
+        return mf.impl.getSquaringValue(x);
+    }
+
+    HURCHALLA_FORCE_INLINE
+    static SquaringValue squareSV(const MF& mf, SquaringValue sv)
+    {
+        return mf.impl.squareSV(sv);
+    }
+
+    HURCHALLA_FORCE_INLINE
+    static MontgomeryValue
+    squareToMontgomeryValue(const MF& mf, SquaringValue sv)
+    {
+        return mf.impl.squareToMontgomeryValue(sv);
+    }
+
+    HURCHALLA_FORCE_INLINE
+    static MontgomeryValue getMontgomeryValue(const MF& mf, SquaringValue sv)
+    {
+        return mf.impl.getMontgomeryValue(sv);
+    }
 };
 
 

montgomery_arithmetic/include/hurchalla/montgomery_arithmetic/detail/MontyCommonBase.h

@@ -99,10 +99,8 @@ class MontyCommonBase {
         namespace hc = ::hurchalla;
         bool isNegative;
         T result = hc::REDC_incomplete(isNegative, u_hi, u_lo, n_, inv_n_);
-        HPBC_CLOCKWORK_ASSERT2(result == 0 ?
-                               isNegative == false : isNegative == true);
-        // We would expect that result == 0 is usually very rare, so going by
-        // the assert above, the next 'if' should be well predicted.
+        // Because u_hi == 0, we should expect the following 'if' to be very
+        // well predicted.
         if (isNegative)
             result = static_cast<T>(result + n_);
         HPBC_CLOCKWORK_ASSERT2(result==hc::REDC_standard(u_hi,u_lo,n_,inv_n_));
@@ -444,6 +442,32 @@ class MontyCommonBase {
 
 
 
+    // returns (R*R) mod N
+    HURCHALLA_FORCE_INLINE C getMontvalueR() const
+    {
+        HPBC_CLOCKWORK_INVARIANT2(r_squared_mod_n_ < n_);
+        return C(r_squared_mod_n_);
+    }
+    template <class PTAG> HURCHALLA_FORCE_INLINE
+    V twoPowLimited_times_x(size_t exponent, C cx, PTAG) const
+    {
+        static constexpr int digitsT = ut_numeric_limits<T>::digits;
+        int power = static_cast<int>(exponent);
+        HPBC_CLOCKWORK_PRECONDITION2(0 <= power && power < digitsT);
+
+        T tmp = cx.get();
+        HPBC_CLOCKWORK_INVARIANT2(tmp < n_);
+        T u_lo = static_cast<T>(tmp << power);
+        int rshift = digitsT - power;
+        HPBC_CLOCKWORK_ASSERT2(rshift > 0);
+        T u_hi = (tmp >> 1) >> (rshift - 1);
+
+        HPBC_CLOCKWORK_ASSERT2(u_hi < n_);
+        const D* child = static_cast<const D*>(this);
+        V result = child->montyREDC(u_hi, u_lo, PTAG());
+        HPBC_CLOCKWORK_POSTCONDITION2(child->isValid(result));
+        return result;
+    }
     // returns (R*R*R) mod N
     HURCHALLA_FORCE_INLINE T getMagicValue() const
     {

montgomery_arithmetic/include/hurchalla/montgomery_arithmetic/detail/MontyFullRange.h

@@ -61,6 +61,21 @@ struct MontyFRValueTypes {
         template <typename> friend class MontyFullRange;
         HURCHALLA_FORCE_INLINE explicit FV(T a) : V(a) {}
     };
+    // squaring value type - used for square() optimizations (fyi, those
+    // optimizations wouldn't help much or at all for monty types other than
+    // MontyFullRange).
+    struct SV {
+        HURCHALLA_FORCE_INLINE SV() = default;
+     protected:
+        template <typename> friend class MontyFullRange;
+        HURCHALLA_FORCE_INLINE T getbits() const { return bits; }
+        HURCHALLA_FORCE_INLINE T get_subtrahend() const { return subtrahend; }
+        HURCHALLA_FORCE_INLINE SV(T bbits, T subt) :
+                                                bits(bbits), subtrahend(subt) {}
+     private:
+        T bits;
+        T subtrahend;
+    };
 };
 
 
@@ -77,12 +92,14 @@ class MontyFullRange final :
     using typename BC::V;
     using typename BC::C;
     using FV = typename MontyFRValueTypes<T>::FV;
+    using SV = typename MontyFRValueTypes<T>::SV;
  public:
     using MontyTag = TagMontyFullrange;
     using uint_type = T;
     using montvalue_type = V;
     using canonvalue_type = C;
     using fusingvalue_type = FV;
+    using squaringvalue_type = SV;
 
     explicit MontyFullRange(T modulus) : BC(modulus) {}
 
@@ -193,6 +210,59 @@ class MontyFullRange final :
         return add(x, x);
     }
 
+
+    HURCHALLA_FORCE_INLINE SV getSquaringValue(V x) const
+    {
+        return SV(x.get(), 0);
+    }
+
+    HURCHALLA_FORCE_INLINE SV squareSV(SV sv) const
+    {
+        // see squareToHiLo in MontyFullRangeMasked.h for basic ideas of
+        // proof for why u_hi and u_lo are correct
+        namespace hc = ::hurchalla;
+        T a = sv.getbits();
+        T sqlo;
+        T sqhi = hc::unsigned_multiply_to_hilo_product(sqlo, a, a);
+        T a_or_zero = sv.get_subtrahend();
+        T u_hi = static_cast<T>(sqhi - a_or_zero - a_or_zero);
+        T u_lo = sqlo;
+        HPBC_CLOCKWORK_ASSERT2(u_hi < n_);
+
+        bool isNegative;
+        T res = hc::REDC_incomplete(isNegative, u_hi, u_lo, n_, BC::inv_n_);
+        T subtrahend = isNegative ? res : static_cast<T>(0);
+        SV result(res, subtrahend);
+        return result;
+    }
+
+    HURCHALLA_FORCE_INLINE V squareToMontgomeryValue(SV sv) const
+    {
+        // see squareToHiLo in MontyFullRangeMasked.h for basic ideas of
+        // proof for why u_hi and u_lo are correct
+        namespace hc = ::hurchalla;
+        T a = sv.getbits();
+        T sqlo;
+        T sqhi = hc::unsigned_multiply_to_hilo_product(sqlo, a, a);
+        T a_or_zero = sv.get_subtrahend();
+        T u_hi = static_cast<T>(sqhi - a_or_zero - a_or_zero);
+        T u_lo = sqlo;
+        HPBC_CLOCKWORK_ASSERT2(u_hi < n_);
+
+        T res = hc::REDC_standard(
+                               u_hi, u_lo, n_, BC::inv_n_, hc::LowlatencyTag());
+        V result(res);
+        return result;
+    }
+
+    // probably I would not want to use this, instead preferring to get a SV
+    // via squareToMontgomeryValue
+    HURCHALLA_FORCE_INLINE V getMontgomeryValue(SV sv) const
+    {
+        T nonneg_value = sv.get_subtrahend() != 0 ? sv.get() + n_ : sv.get();
+        return V(nonneg_value);
+    }
+
 private:
     // functions called by the 'curiously recurring template pattern' base (BC).
     friend BC;

montgomery_arithmetic/include/hurchalla/montgomery_arithmetic/detail/MontyHalfRange.h

@@ -113,12 +113,14 @@ class MontyHalfRange final :
     using typename BC::V;
     using typename BC::C;
     using FV = typename MontyHRValueTypes<T>::FV;
+    using SV = V;
  public:
     using MontyTag = TagMontyHalfrange;
     using uint_type = T;
     using montvalue_type = V;
     using canonvalue_type = C;
     using fusingvalue_type = FV;
+    using squaringvalue_type = SV;
 
     explicit MontyHalfRange(T modulus) : BC(modulus)
     {
@@ -490,6 +492,28 @@ class MontyHalfRange final :
         return C(result);
     }
 
+
+    HURCHALLA_FORCE_INLINE SV getSquaringValue(V x) const
+    {
+        static_assert(std::is_same<V, SV>::value, "");
+        return x;
+    }
+    HURCHALLA_FORCE_INLINE SV squareSV(SV sv) const
+    {
+        static_assert(std::is_same<V, SV>::value, "");
+        return BC::square(sv, LowlatencyTag());
+    }
+    HURCHALLA_FORCE_INLINE V squareToMontgomeryValue(SV sv) const
+    {
+        static_assert(std::is_same<V, SV>::value, "");
+        return BC::square(sv, LowlatencyTag());
+    }
+    HURCHALLA_FORCE_INLINE V getMontgomeryValue(SV sv) const
+    {
+        static_assert(std::is_same<V, SV>::value, "");
+        return sv;
+    }
+
 private:
     // functions called by the 'curiously recurring template pattern' base (BC).
     friend BC;
@@ -500,8 +524,7 @@ class MontyHalfRange final :
     {
         HPBC_CLOCKWORK_PRECONDITION2(u_hi < n_);  // verifies that (u_hi*R + u_lo) < n*R
         namespace hc = ::hurchalla;
-        bool isNegative;  // ignored
-        T result = hc::REDC_incomplete(isNegative, u_hi, u_lo, n_, BC::inv_n_);
+        T result = hc::REDC_incomplete(u_hi, u_lo, n_, BC::inv_n_);
         resultIsZero = (result == 0);
         V v = V(static_cast<S>(result));
         HPBC_CLOCKWORK_POSTCONDITION2(isValid(v));