Study on Improving the Performance of BigInteger.ModPow

[Quicksilver0218]

In my opinion, the performance of BigInteger.ModPow is too bad. Therefore, I have studied the current implementation and tried to find ways to improve it. I would like to share my findings in this discussion. I found an article that mentions some modular exponential algorithms. The current implementation is also included in the article as the right-to-left (LSB-to-MSB) method of modular exponentiation by squaring. Here is the current implementation in BigIntegerCalculator.ModPow.cs

private static Span<uint> PowCore(Span<uint> value, int valueLength,
                                  ReadOnlySpan<uint> power, ReadOnlySpan<uint> modulus,
                                  Span<uint> result, int resultLength,
                                  Span<uint> temp)
{
    // The big modulus pow algorithm for all but
    // the last power limb using square-and-multiply.

    // NOTE: we're using an ordinary remainder here,
    // since the reducer overhead doesn't pay off.

    for (int i = 0; i < power.Length - 1; i++)
    {
        uint p = power[i];
        for (int j = 0; j < 32; j++)
        {
            if ((p & 1) == 1)
            {
                resultLength = MultiplySelf(ref result, resultLength, value.Slice(0, valueLength), ref temp);
                resultLength = Reduce(result.Slice(0, resultLength), modulus);
            }
            valueLength = SquareSelf(ref value, valueLength, ref temp);
            valueLength = Reduce(value.Slice(0, valueLength), modulus);
            p >>= 1;
        }
    }

    return PowCore(value, valueLength, power[power.Length - 1], modulus, result, resultLength, temp);
}

I have tried to rewrite the method to implement other algorithms mentioned in the article and carry out testing on them.
Here is the binary MSB-to-LSB method:

private static Span<uint> PowCore(Span<uint> value, int valueLength,
                                  ReadOnlySpan<uint> power, ReadOnlySpan<uint> modulus,
                                  Span<uint> result, int resultLength,
                                  Span<uint> temp)
{
    bool msbFound = false;
    for (int i = power.Length - 1; i >= 0; i--)
    {
        uint p = power[i];
        for (int _ = 0; _ < 32; _++)
        {
            if (msbFound)
            {
                resultLength = SquareSelf(ref result, resultLength, ref temp);
                resultLength = Reduce(result[..resultLength], modulus);
            }
            if ((p & 0x80000000) != 0)
                if (msbFound)
                {
                    resultLength = MultiplySelf(ref result, resultLength, value[..valueLength], ref temp);
                    resultLength = Reduce(result[..resultLength], modulus);
                }
                else
                {
                    value[..valueLength].CopyTo(result);
                    resultLength = valueLength;
                    msbFound = true;
                }
            p <<= 1;
        }
    }

    return result[..resultLength];
}

Here is the k-ary MSB-to-LSB method, k can be set to different values (k should be a factor of 32 so that all the shifting operations can be done within a uint):

private static Span<uint> PowCore(Span<uint> value, int valueLength,
                                  ReadOnlySpan<uint> power, ReadOnlySpan<uint> modulus,
                                  Span<uint> result, int resultLength,
                                  Span<uint> temp, int k = 1)
{
    Debug.Assert(32 % k == 0);

    uint[][] table = new uint[(1 << k) - 1][];
    table[0] = value[..valueLength].ToArray();
    for (int i = 1; i < table.Length; i++)
    {
        valueLength = MultiplySelf(ref value, valueLength, table[0], ref temp);
        valueLength = Reduce(value[..valueLength], modulus);
        table[i] = value[..valueLength].ToArray();
    }

    bool msbFound = false;
    for (int i = power.Length - 1; i >= 0; i--)
        for (int j = 32 - k; j >= 0; j -= k)
        {
            if (msbFound)
                for (int _ = 0; _ < k; _++)
                {
                    resultLength = SquareSelf(ref result, resultLength, ref temp);
                    resultLength = Reduce(result[..resultLength], modulus);
                }
            uint digit = (uint)(power[i] >>> j & table.Length);
            if (digit != 0)
            {
                uint index = digit - 1;
                if (msbFound)
                {
                    resultLength = MultiplySelf(ref result, resultLength, table[index], ref temp);
                    resultLength = Reduce(result[..resultLength], modulus);
                }
                else
                {
                    table[index].CopyTo(result);
                    resultLength = table[index].Length;
                    msbFound = true;
                }
            }
        }

    return result[..resultLength];
}

After some experiments, I found that the binary MSB-to-LSB method runs faster with a smaller value, a greater exponent or a greater modulus than the current implementation. However, I cannot find the exact threshold. Regarding the k-ary MSB-to-LSB method, its performance is similar to that of the binary MSB-to-LSB method. Maybe I have made an incorrect implementation. Although the MSB-to-LSB algorithm is faster in some cases, the whole ModPow operation is still too slow compared with other languages. e.g. ~10 times slower than Java. I think further improvements can be made by replacing the MultiplySelf / SquareSelf + Reduce parts with other algorithms such as Montgomery modular multiplication.

[tannergooding]

As has been covered in many other posts, the primary issue with BigInteger today is that it is using 32-bit storage on 64-bit systems. This results in significantly more complex/expensive multiplication and division algorithms, it also results in effectively twice as many operations being done for things like simple addition/subtraction, negation, shifting, bitwise operations, etc.

There are additional things that can be improved as well, but the primary thing is a much needed rewrite of the general type to use nuint[] rather than uint[] as the underlying storage, so that we can best take advantage of the hardware capabilities.

Such a rewrite is incredibly expensive and involved, its something that needs top down approval and dedicated time to be worked on by the team. -- CC. @jeffhandley