Fix bug of FastReducer used in BigInteger.ModPow

[key-moon]

There is a bug in the current implementation of FastReducer that causes a runtime error when certain conditions described below are met.

• Modulus is less than or equal to base^{modulus.Length-1}.
• Length of divided value must be sufficiently large when call reduce (this length depending on the modulus).

When these conditions are met, the size of the buffers q1 and q2 used in the algorithm becomes insufficient, causing a assertion error on BigIntegerCalculator.FastReducer.cs#L77 in debug build, or IndexOutOfRangeException on BigIntegerCalculator.PowMod.cs#L401 in Release build.

FastReducer is used only in BigInteger.ModPow, but the specific condition cause above error. Following code is one of the example.

BigInteger.ModPow(BigInteger.One << 1008, 2, BigInteger.One << 992);

This PR fixes this bug by changing the allocated size of the buffer. The validity of the change is explained below.

Edit: After submitting this PR, I found a bug that is related but has a different cause, so I will fix that as well. See the below comments for details.

Explanation

Unless otherwise stated, all variables are non-negative and division is defined as truncated division.

Before we begin, we will define a few constants and functions.

               base := 2^32
      Length(value) := Length of the current array representing the value
ActualLength(value) := Minimum array length that needs to represent value

First, we will look at the algorithm of the constructor.

given m: uint array notation of modulus

let k := Length(m)
let r := base^(2k)

let mu := r / m
let muLength := ActualLength(mu)

Note: The variable k is implicitly defined and does not appear in the actual code.

The above corresponds to the following code.

runtime/src/libraries/System.Runtime.Numerics/src/System/Numerics/BigIntegerCalculator.FastReducer.cs

Lines 26 to 43 in 949fef4

	public FastReducer(uint[] modulus)
	{
	Debug.Assert(modulus != null);
	// Let r = (2^32)^(2k), with (2^32)^k > m
	uint[] r = new uint[modulus.Length * 2 + 1];
	r[r.Length - 1] = 1;
	// Let mu = r / m
	_mu = Divide(r, modulus);
	_modulus = modulus;
	// Allocate memory for quotients once
	_q1 = new uint[modulus.Length * 2 + 2];
	_q2 = new uint[modulus.Length * 2 + 1];
	_muLength = ActualLength(_mu);
	}

Next, we will look at the algorithm of the Reduce function.

given V: value that ActualLength(V) is less than or equal 2k
      ⇔ V is less than base^(2k)

let q1 := V / base^(k-1) * mu
let q2 := q1 / base^(k+1) * m

The above corresponds to the following code.

runtime/src/libraries/System.Runtime.Numerics/src/System/Numerics/BigIntegerCalculator.FastReducer.cs

Lines 45 to 66 in 949fef4

	public int Reduce(uint[] value, int length)
	{
	Debug.Assert(value != null);
	Debug.Assert(length <= value.Length);
	Debug.Assert(value.Length <= _modulus.Length * 2);
	// Trivial: value is shorter
	if (length < _modulus.Length)
	return length;
	// Let q1 = v/(2^32)^(k-1) * mu
	int l1 = DivMul(value, length, _mu, _muLength,
	_q1, _modulus.Length - 1);
	// Let q2 = q1/(2^32)^(k+1) * m
	int l2 = DivMul(_q1, l1, _modulus, _modulus.Length,
	_q2, _modulus.Length + 1);
	// Let v = (v - q2) % (2^32)^(k+1) - i*m
	return SubMod(value, length, _q2, l2,
	_modulus, _modulus.Length + 1);
	}

Here, we will think about upper bound of ActualLength(q1).
Since upper bound of V is base^(2k), upper bound of V / base^(k-1) is base^(k+1).
Since modulus should not be 0, r and mu can't be 0. So mu is less than base^(muLength).
It is easy to show that an upper bound on the values of both V and mu can be achieved.
From these facts, we can see that the upper bound of q1 is base^(muLength + k + 1).

Same as q1, we can see that upper bound of q2 is base^(muLength + k)

Now we have shown the upper bound of q1 and q2. From these, we can see that the upper bound of length of the buffer required for q1 and q2 is: ActualLength(q1) ≦ muLength + k + 1, ActualLength(q2) ≦ muLength + k.

[ghost]

Tagging subscribers to this area: @tannergooding, @pgovind
See info in area-owners.md if you want to be subscribed.

Issue Details

---

There is a bug in the current implementation of FastReducer that causes a runtime error when certain conditions described below are met.

• The value represented by modulus is less than or equal to base^{modulus.Length-1}.
• Length of value must be sufficiently large when call reduce (this length depending on the modulus).

When these conditions are met, the size of the buffers q1 and q2 used in the algorithm becomes insufficient, causing a assertion error on BigIntegerCalculator.FastReducer.cs#L77 in debug build, or IndexOutOfRangeException on BigIntegerCalculator.PowMod.cs#L401 in Release build.

FastReducer is used only in BigInteger.ModPow, but the specific condition cause above error. Following code is one of the example.

BigInteger.ModPow(BigInteger.One << 1008, 2, BigInteger.One << 992);

This PR fixes this bug by changing the allocated size of the buffer. The validity of the change is explained below.

Explanation

Unless otherwise stated, all variables are non-negative and division is defined as truncated division.

Before we begin, we will define a few constants and functions.

               base := 2^32
      Length(value) := Length of the current array representing the value
ActualLength(value) := Minimum array length that needs to represent value

First, we will look at the implementation of the constructor.

given m: uint array notation of modulus

let k := Length(m)
let r := base^(2k)

let mu := r / m
let muLength := ActualLength(mu)

Note: The variable k is implicitly defined and does not appear in the actual code.

The above corresponds to the following code.

runtime/src/libraries/System.Runtime.Numerics/src/System/Numerics/BigIntegerCalculator.FastReducer.cs

Lines 26 to 43 in 949fef4

	public FastReducer(uint[] modulus)
	{
	Debug.Assert(modulus != null);
	// Let r = (2^32)^(2k), with (2^32)^k > m
	uint[] r = new uint[modulus.Length * 2 + 1];
	r[r.Length - 1] = 1;
	// Let mu = r / m
	_mu = Divide(r, modulus);
	_modulus = modulus;
	// Allocate memory for quotients once
	_q1 = new uint[modulus.Length * 2 + 2];
	_q2 = new uint[modulus.Length * 2 + 1];
	_muLength = ActualLength(_mu);
	}

Next, we will look at the implementation of the Reduce function.

given V: value that ActualLength(V) is less than or equal 2k
      ⇔ V is less than base^(2k)

let q1 := V / base^(k-1) * mu
let q2 := q1 / base^(k+1) * m

The above corresponds to the following code.

runtime/src/libraries/System.Runtime.Numerics/src/System/Numerics/BigIntegerCalculator.FastReducer.cs

Lines 45 to 66 in 949fef4

	public int Reduce(uint[] value, int length)
	{
	Debug.Assert(value != null);
	Debug.Assert(length <= value.Length);
	Debug.Assert(value.Length <= _modulus.Length * 2);
	// Trivial: value is shorter
	if (length < _modulus.Length)
	return length;
	// Let q1 = v/(2^32)^(k-1) * mu
	int l1 = DivMul(value, length, _mu, _muLength,
	_q1, _modulus.Length - 1);
	// Let q2 = q1/(2^32)^(k+1) * m
	int l2 = DivMul(_q1, l1, _modulus, _modulus.Length,
	_q2, _modulus.Length + 1);
	// Let v = (v - q2) % (2^32)^(k+1) - i*m
	return SubMod(value, length, _q2, l2,
	_modulus, _modulus.Length + 1);
	}

Here, we will think about upper bound of ActualLength(q1).
Since upper bound of V is base^(2k), upper bound of V / base^(k-1) is base^(k+1).
Since modulus should not be 0, r and mu can't be 0. So mu is less than base^(muLength).
It is easy to show that an upper bound on the values of both V and mu can be achieved.
From these facts, we can see that the upper bound of q1 is base^(muLength + k + 1).

Same as q1, we can see that upper bound of q2 is base^(muLength + k)

Now we have shown the upper bound of q1 and q2. From these, we can see that the upper bound of length of the buffer required for q1 and q2 is: ActualLength(q1) < muLength + k + 1, ActualLength(q2) < muLength + k.

Author:	key-moon
Assignees:	-
Labels:	area-System.Numerics
Milestone:	-

[key-moon]

There is one additional bug in the current FastReducer implementation. This bug occurs when the lower half bit of the value is less than value%modulus.
This happens when you call BigInteger.ModPow using the FastReducer as follows:

BigInteger.ModPow(BigInteger.One << 1008, 2, (BigInteger.One << 992) + 1);

The above commit adds a test that covers this case and fixes the bug. Below is an explanation of the cause and fix for this bug.

Explanation

In this explanation, the symbols used in the above explanation are reused.
Barett Reduction is an algorithm that takes advantage of the fact that V mod m = V - ⌊v/m⌋m holds. In this algorithm, value of q2 is ⌊v/m⌋m. Now we will look at the algorithm of the last subtraction.

v := (v - q2) % base^(k+1)
while m ≦ v:
  v := v - m

The above corresponds to the following code.

runtime/src/libraries/System.Runtime.Numerics/src/System/Numerics/BigIntegerCalculator.FastReducer.cs

Lines 63 to 65 in 949fef4

	// Let v = (v - q2) % (2^32)^(k+1) - i*m
	return SubMod(value, length, _q2, l2,
	_modulus, _modulus.Length + 1);

runtime/src/libraries/System.Runtime.Numerics/src/System/Numerics/BigIntegerCalculator.FastReducer.cs

Lines 112 to 146 in 949fef4

	private static unsafe int SubMod(uint[] left, int leftLength,
	uint[] right, int rightLength,
	uint[] modulus, int k)
	{
	Debug.Assert(left != null);
	Debug.Assert(left.Length >= leftLength);
	Debug.Assert(right != null);
	Debug.Assert(right.Length >= rightLength);
	// Executes the subtraction algorithm for left and right,
	// but considers only the first k limbs, which is equivalent to
	// preceding reduction by 2^(32*k). Furthermore, if left is
	// still greater than modulus, further subtractions are used.
	if (leftLength > k)
	leftLength = k;
	if (rightLength > k)
	rightLength = k;
	fixed (uint* l = left, r = right, m = modulus)
	{
	SubtractSelf(l, leftLength, r, rightLength);
	leftLength = ActualLength(left, leftLength);
	while (Compare(l, leftLength, m, modulus.Length) >= 0)
	{
	SubtractSelf(l, leftLength, m, modulus.Length);
	leftLength = ActualLength(left, leftLength);
	}
	}
	Array.Clear(left, leftLength, left.Length - leftLength);
	return leftLength;
	}

The reason we divide by base^(k+1) first is that we are trying to speed up the process by using the fact that mod is greater than base^(k+1) by the definition of k. This algorithm itself is reasonable and there is no problem.
However, the actual implementation process is as follows:

given left, right, modulus, k

left := left % base^k
right := right % base^k
left := left - right
while modulus ≦ left:
  left := left - modulus

The SubtractSelf function, which is performing the subtraction in line 133, does not allow the subtrahend to be larger than minuend. Therefore, this code may result in an assertion error.
This can be avoided by using the SubtractSelf function which allows overflow. The reason this works is that the overflow in leftLength=k acts as operation on Z/(base^k)Z. Note that the overflow occurs only when left≧base^k, since left≧right holds.

Furthermore, since v%m<m<base^k holds, we know that the number to divide v can be base^k, not base^(k+1).

[key-moon]

I accidentally deleted the forked repository, so I'll close this one and repost it. I apologize for the inconvenience caused by this.