Add BigInteger.ModInverse and BigInteger.LeastCommonMultiple

[su-senka]

Add BigInteger.ModInverse and BigInteger.LeastCommonMultiple

Implements the two APIs approved in #130021:

public static BigInteger ModInverse(BigInteger value, BigInteger modulus);
public static BigInteger LeastCommonMultiple(BigInteger left, BigInteger right);

Both match the shape signed off by @bartonjs in API review. Nothing outside the approved
surface is added — in particular there is no ExtendedEuclidian and no TryModInverse.

Semantics

ModInverse(value, modulus)

• Returns the least non-negative x in [0, modulus) with value * x == 1 (mod modulus).
• value may be negative; it is reduced to its canonical representative first.
• Returns Zero when modulus is 1, for any value including 0.
• Throws ArgumentOutOfRangeException when modulus <= 0, with ParamName "modulus".
• Throws ArithmeticException when value and modulus are not coprime, which includes
  the case where value is congruent to zero modulo modulus.

LeastCommonMultiple(left, right)

The proposal did not pin down the edge cases, so these were posted on the issue and drew
no objection before implementation started:

• The result is always non-negative. The sign of either operand does not change the set of
  their common multiples.
• The result is Zero when either operand is zero, including (0, 0). Zero is a multiple
  of every integer, so it is the least common multiple in that case rather than an error.
• The method never throws.

It is computed as Abs(left / gcd) * Abs(right), dividing before multiplying so the
intermediate value never exceeds the result.

Documentation wording

The <returns> text on ModInverse states the contract as a congruence:

the least non-negative integer x in the range [0, modulus) that satisfies
value * x == 1 (mod modulus)

rather than the (value * x) % modulus == 1 phrasing used in the original proposal text.
That phrasing is wrong. C#'s % truncates toward zero, so for negative value the
remainder is negative and the stated equality does not hold even when the returned x is
correct. ModInverse(-3, 7) is 2, and (-3 * 2) % 7 is -6, not 1. The congruence
form is true for every input the method accepts.

Implementation

The work landed in five reviewable steps, one commit each.

Schoolbook baseline first. ModInverse was first implemented as a plain extended
Euclidean algorithm over BigInteger operations, purely to establish correctness. It was
then retained temporarily as a debug-only differential oracle — a Debug.Assert comparing
every result computed by the real implementation against the schoolbook one — while the
optimized paths were developed. The final commit removes both the oracle and the second
implementation, so this is development scaffolding that does not ship.

Lehmer-based extended GCD. The general path is BigIntegerCalculator.ModInverse in
BigIntegerCalculator.GcdInv.cs, reusing the existing Lehmer guessing step and Jebelean
termination condition from Gcd. The Bézout cofactor is carried alongside the remainder
sequence using the recurrence

t[0] = 0, t[1] = 1, t[i+1] = t[i-1] + q[i] * t[i]

which satisfies r[i] == (-1)^(i+1) * t[i] * value (mod modulus). Every cofactor is
non-negative and the recurrence only ever adds, so the cofactors stay unsigned like the
rest of the calculator and the alternating sign is carried as the parity of the remainder
index alone.

This falls out especially cleanly for the Lehmer step. LehmerCore applies
x' = a*x - b*y, y' = d*y - c*x to the remainders. Cofactors take the same linear map,
and substituting the alternating signs cancels both subtractions exactly, leaving
u' = a*u + b*v, v' = c*u + d*v on magnitudes — the same matrix with the signs dropped,
not its transpose. The existing odd-iteration swap that restores remainder ordering
doubles as the parity flip, so no separate sign bookkeeping is needed.

Buffer sizing comes from the identity t[i] * r[i-1] + t[i-1] * r[i] == modulus, which
bounds every cofactor by the modulus, including while one is being accumulated. Cofactor
buffers therefore carry exactly one spare limb.

Hensel lifting for power-of-two moduli. Inverting modulo 2^k does not need the
Lehmer path. ModInversePowerOfTwo runs Newton iteration on the 2-adic inverse: if
value * x == 1 (mod 2^m) then x * (2 - value * x) is correct modulo 2^2m, because
writing value * x = 1 + e * 2^m gives (1 + e*2^m) * (1 - e*2^m) = 1 - e^2 * 2^2m. Each
step doubles the number of correct bits, so the lift terminates after
ceil(log2(limbCount)) steps.

The seed is the inverse modulo a single limb, taken from the existing
ComputeMontgomeryInverse in BigIntegerCalculator.PowMod.cs. That helper returns the
negated inverse, because that is what Montgomery reduction consumes, so the sign is
corrected before use. The lift works in whole limbs and the top limb is then masked back
to the bits the modulus spans, which is sound because 2^k divides the limb radix raised
to the limb count.

Note for future contributors to this file

BigIntegerCalculator.Multiply requires a destination of exactly left.Length + right.Length limbs, not merely one that is large enough. The documented assert only says
bits.Length >= left.Length + right.Length, but the Karatsuba path splits the destination
as bits[..2n] / bits[2n..] and asserts bitsLow.Length >= bitsHigh.Length, which fails
for an oversized destination. Reusing one generously sized scratch buffer across several
multiplies of differing widths — the obvious thing to do — trips this. It is worth knowing
before writing new code against that API.

32-bit

One branch in this change is architecture-conditional and was never executed locally,
because this was developed and tested on arm64:

• BigIntegerCalculator.GcdInv.cs:321 — int lehmerMinimum = nint.Size == 4 ? 3 : 2;

This is the minimum operand length before Lehmer's guess can be attempted, mirroring
ExtractDigits' own asserts (>= 3 limbs on 32-bit, >= 2 on 64-bit). Too low and
ExtractDigits reads past the operand; too high only costs performance. It deserves
attention from 32-bit CI legs.

The surface was deliberately kept this small. The cofactor combination helpers use a
UInt128 accumulator with BitsPerLimb shifts and so are architecture-neutral, and the
Euclid step was written as a fully general fallback used whenever operands are too short
for Lehmer, rather than as a separate terminal path that only 32-bit would reach. The
power-of-two fast path introduces no architecture-conditional branch at all.

Testing

Two new test files, modinverse.cs and leastcommonmultiple.cs, add 84 tests. The full
System.Runtime.Numerics suite is 8718 tests, 0 failures, up from a baseline of 8634.

The tests verify the defining mathematical identities directly rather than diffing against
a reference implementation: the congruence value * x == 1 (mod modulus), canonicality of
the representative, ModInverse being an involution, agreement with ModPow via Fermat's
little theorem for Mersenne prime moduli, lcm(a,b) * gcd(a,b) == |a*b|, coprimality of
the LCM cofactors, associativity, commutativity, and identity. This means they remain
meaningful independently of how the implementation is written.

Coverage includes exhaustive power-of-two bit lengths — every length through 300 plus the
larger limb boundaries — operand widths from 1 to 512 bytes, asymmetric operand sizes,
Mersenne and shifted-limb structured operands, and the non-coprime and non-positive-modulus
error paths.

Mutation testing was used to check the tests are actually sensitive to the errors most
likely to occur here, since a green suite alone says little about a sign convention:

• Inverting the final parity: 10+ failures.
• Removing the parity flip only on the Lehmer odd-iteration path, keeping the swap — a bug
  reachable only on large operands: 6 failures, including LargeOperands and
  InverseOfInverseIsReducedValue.
• Dropping the ComputeMontgomeryInverse sign correction: caught.

Separately, during development an out-of-tree stress harness exercised roughly 480,000
exhaustive small cases, 300 consecutive Fibonacci pairs (the worst case for Euclid step
counts, and so for parity churn), and randomized operands across 26 widths, with every call
cross-checked against the schoolbook oracle described above.

Performance

Measured on a Release build of clr+libs. The harness prints the resolved
System.Runtime.Numerics location and its DebuggableAttribute optimizer status before
running, and refuses to present numbers from a Debug assembly — this code carries per-limb
Debug.Asserts, including an O(n) buffer scan, that would otherwise dominate.

BenchmarkDotNet v0.15.2, macOS 26.6 (25G72) [Darwin 25.6.0]
Apple M4, 1 CPU, 10 logical and 10 physical cores
  [Host] : .NET 11.0.0 (42.42.42.42424), Arm64 RyuJIT AdvSIMD
Toolchain=InProcessEmitToolchain

Method	Bits	Mean	Error	StdDev	Ratio	Gen0	Allocated
GreatestCommonDivisor	256	809.2 ns	4.64 ns	4.34 ns	1.00	-	-
ModInverse	256	1,679.4 ns	2.06 ns	1.93 ns	2.08	0.0134	112 B
ModInverse (2^k modulus)	256	122.1 ns	1.26 ns	1.17 ns	0.15	0.0067	56 B
LeastCommonMultiple	256	955.0 ns	1.62 ns	1.52 ns	1.18	0.0172	144 B
GreatestCommonDivisor	2048	8,083.4 ns	13.68 ns	12.13 ns	1.00	-	-
ModInverse	2048	16,547.0 ns	18.55 ns	17.35 ns	2.05	0.0610	560 B
ModInverse (2^k modulus)	2048	1,435.1 ns	4.60 ns	4.08 ns	0.18	0.0324	280 B
LeastCommonMultiple	2048	9,199.4 ns	24.11 ns	20.13 ns	1.14	0.0916	816 B
GreatestCommonDivisor	4096	18,974.7 ns	30.60 ns	28.62 ns	1.00	-	-
ModInverse	4096	50,605.2 ns	26.36 ns	24.66 ns	2.67	0.1221	1,072 B
ModInverse (2^k modulus)	4096	4,237.6 ns	5.56 ns	5.20 ns	0.22	0.0610	536 B
LeastCommonMultiple	4096	22,211.7 ns	20.16 ns	15.74 ns	1.17	0.1831	1,584 B

Takeaways:

• Cofactor tracking costs about 2x plain GCD — 2.08x, 2.05x, 2.67x at 256, 2048 and
  4096 bits. GreatestCommonDivisor is the baseline precisely so this overhead is visible;
  it is the price of the extended algorithm over the plain one on the shared Lehmer path.
  The rise to 2.67x at 4096 bits is outside the noise: cofactors grow toward modulus width
  while remainders shrink, so cofactor updates become a larger share of each Lehmer block
  as operands grow.
• The power-of-two fast path is 12-14x faster than the general path, at 0.15-0.22x the
  cost of a GCD. That is Hensel lifting's O(M(n) log n) against Lehmer's O(n^2).
• LeastCommonMultiple costs GCD plus about 15%, consistently — one exact division and
  one multiplication on top of the GCD it has to compute anyway.

The GreatestCommonDivisor baseline shows zero allocation because its operands are coprime
by construction, as they must be for an inverse to exist, so it returns 1, which
BigInteger represents inline. The work is a full-length Lehmer run either way; only the
allocation column is asymmetric.

[azure-pipelines]

Azure Pipelines:
Successfully started running 3 pipeline(s).
13 pipeline(s) were filtered out due to trigger conditions.
There may be pipelines that require an authorized user to comment /azp run to run.

[dotnet-policy-service]

Tagging subscribers to this area: @dotnet/area-system-numerics
See info in area-owners.md if you want to be subscribed.