Improve performance of BigInteger.Pow/ModPow

- To introduce further performance tweaks, the exponentiation algorithms
  are rewritten and ported to `BigIntegerCalculator`.
- To scale better for bigger numbers a `FastReducer` is in use instead
  of "ordinary" modulo operations, which is based on multiplications.
- Furthermore the newly introduced `FastReducer` triggers a bad corner
  case within the division algorithm, which gets fixed too.
- A performance regression at 64 bits within integer division was found,
  which gets fixed too (no more allocations within that code).
- The test code for threshold values of square / multiply / modpow now
  modifies these thresholds for more thorough testing.

Author: axelheer

src/System.Runtime.Numerics/src/System.Runtime.Numerics.csproj

@@ -19,7 +19,10 @@
   </PropertyGroup>
   <ItemGroup>
     <Compile Include="System\Numerics\BigIntegerCalculator.AddSub.cs" />
+    <Compile Include="System\Numerics\BigIntegerCalculator.BitsBuffer.cs" />
     <Compile Include="System\Numerics\BigIntegerCalculator.DivRem.cs" />
+    <Compile Include="System\Numerics\BigIntegerCalculator.FastReducer.cs" />
+    <Compile Include="System\Numerics\BigIntegerCalculator.PowMod.cs" />
     <Compile Include="System\Numerics\BigIntegerCalculator.SquMul.cs" />
     <Compile Include="System\Numerics\BigInteger.cs" />
     <Compile Include="System\Numerics\BigIntegerBuilder.cs" />

src/System.Runtime.Numerics/src/System/Numerics/BigInteger.cs

@@ -1049,45 +1049,6 @@ public static BigInteger Min(BigInteger left, BigInteger right)
             return right;
         }
 
-
-        private static void ModPowUpdateResult(ref BigIntegerBuilder regRes, ref BigIntegerBuilder regVal, ref BigIntegerBuilder regMod, ref BigIntegerBuilder regTmp)
-        {
-            NumericsHelpers.Swap(ref regRes, ref regTmp);
-            regRes.Mul(ref regTmp, ref regVal);   // result = result * value;
-            regRes.Mod(ref regMod);               // result = result % modulus;                       
-        }
-
-        private static void ModPowSquareModValue(ref BigIntegerBuilder regVal, ref BigIntegerBuilder regMod, ref BigIntegerBuilder regTmp)
-        {
-            NumericsHelpers.Swap(ref regVal, ref regTmp);
-            regVal.Mul(ref regTmp, ref regTmp);   // value = value * value;
-            regVal.Mod(ref regMod);               // value = value % modulus;
-        }
-
-        private static void ModPowInner(uint exp, ref BigIntegerBuilder regRes, ref BigIntegerBuilder regVal, ref BigIntegerBuilder regMod, ref BigIntegerBuilder regTmp)
-        {
-            while (exp != 0)          // !(Exponent.IsZero)
-            {
-                if ((exp & 1) == 1)   // !(Exponent.IsEven)
-                    ModPowUpdateResult(ref regRes, ref regVal, ref regMod, ref regTmp);
-                if (exp == 1)         // Exponent.IsOne - we can exit early
-                    break;
-                ModPowSquareModValue(ref regVal, ref regMod, ref regTmp);
-                exp >>= 1;
-            }
-        }
-
-        private static void ModPowInner32(uint exp, ref BigIntegerBuilder regRes, ref BigIntegerBuilder regVal, ref BigIntegerBuilder regMod, ref BigIntegerBuilder regTmp)
-        {
-            for (int i = 0; i < 32; i++)
-            {
-                if ((exp & 1) == 1)   // !(Exponent.IsEven)
-                    ModPowUpdateResult(ref regRes, ref regVal, ref regMod, ref regTmp);
-                ModPowSquareModValue(ref regVal, ref regMod, ref regTmp);
-                exp >>= 1;
-            }
-        }
-
         public static BigInteger ModPow(BigInteger value, BigInteger exponent, BigInteger modulus)
         {
             if (exponent.Sign < 0)
@@ -1098,36 +1059,31 @@ public static BigInteger ModPow(BigInteger value, BigInteger exponent, BigIntege
             exponent.AssertValid();
             modulus.AssertValid();
 
-            int signRes = +1;
-            int signVal = +1;
-            int signMod = +1;
-            bool expIsEven = exponent.IsEven;
-            BigIntegerBuilder regRes = new BigIntegerBuilder(BigInteger.One, ref signRes);
-            BigIntegerBuilder regVal = new BigIntegerBuilder(value, ref signVal);
-            BigIntegerBuilder regMod = new BigIntegerBuilder(modulus, ref signMod);
-            BigIntegerBuilder regTmp = new BigIntegerBuilder(regVal.Size);
+            bool trivialValue = value._bits == null;
+            bool trivialExponent = exponent._bits == null;
+            bool trivialModulus = modulus._bits == null;
 
-            regRes.Mod(ref regMod);   // Handle special case of exponent=0, modulus=1
+            if (trivialModulus)
+            {
+                long bits = trivialValue && trivialExponent ? BigIntegerCalculator.Pow(NumericsHelpers.Abs(value._sign), NumericsHelpers.Abs(exponent._sign), NumericsHelpers.Abs(modulus._sign)) :
+                            trivialValue ? BigIntegerCalculator.Pow(NumericsHelpers.Abs(value._sign), exponent._bits, NumericsHelpers.Abs(modulus._sign)) :
+                            trivialExponent ? BigIntegerCalculator.Pow(value._bits, NumericsHelpers.Abs(exponent._sign), NumericsHelpers.Abs(modulus._sign)) :
+                            BigIntegerCalculator.Pow(value._bits, exponent._bits, NumericsHelpers.Abs(modulus._sign));
 
-            if (exponent._bits == null)
-            {   // exponent fits into an Int32
-                ModPowInner((uint)exponent._sign, ref regRes, ref regVal, ref regMod, ref regTmp);
+                return value._sign < 0 && !exponent.IsEven ? -1 * bits : bits;
             }
             else
-            {   // very large exponent
-                int len = Length(exponent._bits);
-                for (int i = 0; i < len - 1; i++)
-                {
-                    uint exp = exponent._bits[i];
-                    ModPowInner32(exp, ref regRes, ref regVal, ref regMod, ref regTmp);
-                }
-                ModPowInner(exponent._bits[len - 1], ref regRes, ref regVal, ref regMod, ref regTmp);
-            }
+            {
+                uint[] bits = trivialValue && trivialExponent ? BigIntegerCalculator.Pow(NumericsHelpers.Abs(value._sign), NumericsHelpers.Abs(exponent._sign), modulus._bits) :
+                              trivialValue ? BigIntegerCalculator.Pow(NumericsHelpers.Abs(value._sign), exponent._bits, modulus._bits) :
+                              trivialExponent ? BigIntegerCalculator.Pow(value._bits, NumericsHelpers.Abs(exponent._sign), modulus._bits) :
+                              BigIntegerCalculator.Pow(value._bits, exponent._bits, modulus._bits);
 
-            return regRes.GetInteger(value._sign > 0 ? +1 : expIsEven ? +1 : -1);
+                return new BigInteger(bits, value._sign < 0 && !exponent.IsEven);
+            }
         }
 
-        public static BigInteger Pow(BigInteger value, Int32 exponent)
+        public static BigInteger Pow(BigInteger value, int exponent)
         {
             if (exponent < 0)
                 throw new ArgumentOutOfRangeException("exponent", SR.ArgumentOutOfRange_MustBeNonNeg);
@@ -1136,75 +1092,27 @@ public static BigInteger Pow(BigInteger value, Int32 exponent)
             value.AssertValid();
 
             if (exponent == 0)
-                return BigInteger.One;
+                return s_bnOneInt;
             if (exponent == 1)
                 return value;
-            if (value._bits == null)
+
+            bool trivialValue = value._bits == null;
+
+            if (trivialValue)
             {
                 if (value._sign == 1)
                     return value;
                 if (value._sign == -1)
-                    return (exponent & 1) != 0 ? value : 1;
+                    return (exponent & 1) != 0 ? value : s_bnOneInt;
                 if (value._sign == 0)
                     return value;
             }
 
-            int sign = +1;
-            BigIntegerBuilder regSquare = new BigIntegerBuilder(value, ref sign);
-
-            // Get an estimate of the size needed for regSquare and regRes, so we can minimize allocations.
-            int cuSquareMin = regSquare.Size;
-            int cuSquareMax = cuSquareMin;
-            uint uSquareMin = regSquare.High;
-            uint uSquareMax = uSquareMin + 1;
-            if (uSquareMax == 0)
-            {
-                cuSquareMax++;
-                uSquareMax = 1;
-            }
-            int cuResMin = 1;
-            int cuResMax = 1;
-            uint uResMin = 1;
-            uint uResMax = 1;
-
-            for (int expTmp = exponent; ;)
-            {
-                if ((expTmp & 1) != 0)
-                {
-                    MulUpper(ref uResMax, ref cuResMax, uSquareMax, cuSquareMax);
-                    MulLower(ref uResMin, ref cuResMin, uSquareMin, cuSquareMin);
-                }
-
-                if ((expTmp >>= 1) == 0)
-                    break;
-
-                MulUpper(ref uSquareMax, ref cuSquareMax, uSquareMax, cuSquareMax);
-                MulLower(ref uSquareMin, ref cuSquareMin, uSquareMin, cuSquareMin);
-            }
-
-            if (cuResMax > 1)
-                regSquare.EnsureWritable(cuResMax, 0);
-            BigIntegerBuilder regTmp = new BigIntegerBuilder(cuResMax);
-            BigIntegerBuilder regRes = new BigIntegerBuilder(cuResMax);
-            regRes.Set(1);
-
-            if ((exponent & 1) == 0)
-                sign = +1;
-
-            for (int expTmp = exponent; ;)
-            {
-                if ((expTmp & 1) != 0)
-                {
-                    NumericsHelpers.Swap(ref regRes, ref regTmp);
-                    regRes.Mul(ref regSquare, ref regTmp);
-                }
-                if ((expTmp >>= 1) == 0)
-                    break;
-                NumericsHelpers.Swap(ref regSquare, ref regTmp);
-                regSquare.Mul(ref regTmp, ref regTmp);
-            }
+            uint[] bits = trivialValue
+                        ? BigIntegerCalculator.Pow(NumericsHelpers.Abs(value._sign), NumericsHelpers.Abs(exponent))
+                        : BigIntegerCalculator.Pow(value._bits, NumericsHelpers.Abs(exponent));
 
-            return regRes.GetInteger(sign);
+            return new BigInteger(bits, value._sign < 0 && (exponent & 1) != 0);
         }
 
         #endregion public static methods
@@ -1810,8 +1718,8 @@ private static BigInteger Subtract(uint[] leftBits, int leftSign, uint[] rightBi
 
             if (trivialDivisor)
             {
-                uint[] bits = BigIntegerCalculator.Remainder(dividend._bits, NumericsHelpers.Abs(divisor._sign));
-                return new BigInteger(bits, dividend._sign < 0);
+                long bits = BigIntegerCalculator.Remainder(dividend._bits, NumericsHelpers.Abs(divisor._sign));
+                return dividend._sign < 0 ? -1 * bits : bits;
             }
 
             if (dividend._bits.Length < divisor._bits.Length)
@@ -2032,11 +1940,10 @@ private void SetBitsFromDouble(Double value)
         [Pure]
         internal static int Length(uint[] rgu)
         {
-            int cu = rgu.Length;
-            if (rgu[cu - 1] != 0)
-                return cu;
-            Debug.Assert(cu >= 2 && rgu[cu - 2] != 0);
-            return cu - 1;
+            Debug.Assert(rgu[rgu.Length - 1] != 0);
+
+            // no leading zeros
+            return rgu.Length;
         }
 
         internal int _Sign { get { return _sign; } }
@@ -2088,47 +1995,6 @@ private static bool GetPartsForBitManipulation(ref BigInteger x, out uint[] xd,
             return x._sign < 0;
         }
 
-        // Gets an upper bound on the product of uHiRes * (2^32)^(cuRes-1) and uHiMul * (2^32)^(cuMul-1).
-        // The result is put int uHiRes and cuRes.
-        private static void MulUpper(ref uint uHiRes, ref int cuRes, uint uHiMul, int cuMul)
-        {
-            ulong uu = (ulong)uHiRes * uHiMul;
-            uint uHi = NumericsHelpers.GetHi(uu);
-            if (uHi != 0)
-            {
-                if (NumericsHelpers.GetLo(uu) != 0 && ++uHi == 0)
-                {
-                    uHi = 1;
-                    cuRes++;
-                }
-                uHiRes = uHi;
-                cuRes += cuMul;
-            }
-            else
-            {
-                uHiRes = NumericsHelpers.GetLo(uu);
-                cuRes += cuMul - 1;
-            }
-        }
-
-        // Gets a lower bound on the product of uHiRes * (2^32)^(cuRes-1) and uHiMul * (2^32)^(cuMul-1).
-        // The result is put int uHiRes and cuRes.
-        private static void MulLower(ref uint uHiRes, ref int cuRes, uint uHiMul, int cuMul)
-        {
-            ulong uu = (ulong)uHiRes * uHiMul;
-            uint uHi = NumericsHelpers.GetHi(uu);
-            if (uHi != 0)
-            {
-                uHiRes = uHi;
-                cuRes += cuMul;
-            }
-            else
-            {
-                uHiRes = NumericsHelpers.GetLo(uu);
-                cuRes += cuMul - 1;
-            }
-        }
-
         internal static int GetDiffLength(uint[] rgu1, uint[] rgu2, int cu)
         {
             for (int iv = cu; --iv >= 0;)

src/System.Runtime.Numerics/src/System/Numerics/BigIntegerCalculator.AddSub.cs

@@ -197,6 +197,7 @@ private unsafe static void Subtract(uint* left, int leftLength,
                 bits[i] = (uint)digit;
                 carry = digit >> 32;
             }
+
             Debug.Assert(carry == 0);
         }