File: System\Numerics\BigIntegerCalculator.FastReducer.cs
Web Access
Project: src\src\libraries\System.Runtime.Numerics\src\System.Runtime.Numerics.csproj (System.Runtime.Numerics)
// Licensed to the .NET Foundation under one or more agreements.
// The .NET Foundation licenses this file to you under the MIT license.
 
using System.Diagnostics;
 
namespace System.Numerics
{
    internal static partial class BigIntegerCalculator
    {
        // If we need to reduce by a certain modulus again and again, it's much
        // more efficient to do this with multiplication operations. This is
        // possible, if we do some pre-computations first...
 
        // see https://en.wikipedia.org/wiki/Barrett_reduction
 
        private readonly ref struct FastReducer
        {
            private readonly ReadOnlySpan<uint> _modulus;
            private readonly ReadOnlySpan<uint> _mu;
            private readonly Span<uint> _q1;
            private readonly Span<uint> _q2;
 
            public FastReducer(ReadOnlySpan<uint> modulus, Span<uint> r, Span<uint> mu, Span<uint> q1, Span<uint> q2)
            {
                Debug.Assert(!modulus.IsEmpty);
                Debug.Assert(r.Length == modulus.Length * 2 + 1);
                Debug.Assert(mu.Length == r.Length - modulus.Length + 1);
                Debug.Assert(q1.Length == modulus.Length * 2 + 2);
                Debug.Assert(q2.Length == modulus.Length * 2 + 2);
 
                // Let r = 4^k, with 2^k > m
                r[r.Length - 1] = 1;
 
                // Let mu = 4^k / m
                Divide(r, modulus, mu);
                _modulus = modulus;
 
                _q1 = q1;
                _q2 = q2;
 
                _mu = mu.Slice(0, ActualLength(mu));
            }
 
            public int Reduce(Span<uint> value)
            {
                Debug.Assert(value.Length <= _modulus.Length * 2);
 
                // Trivial: value is shorter
                if (value.Length < _modulus.Length)
                    return value.Length;
 
                // Let q1 = v/2^(k-1) * mu
                _q1.Clear();
                int l1 = DivMul(value, _mu, _q1, _modulus.Length - 1);
 
                // Let q2 = q1/2^(k+1) * m
                _q2.Clear();
                int l2 = DivMul(_q1.Slice(0, l1), _modulus, _q2, _modulus.Length + 1);
 
                // Let v = (v - q2) % 2^(k+1) - i*m
                var length = SubMod(value, _q2.Slice(0, l2), _modulus, _modulus.Length + 1);
                value = value.Slice(length);
                value.Clear();
 
                return length;
            }
 
            private static int DivMul(ReadOnlySpan<uint> left, ReadOnlySpan<uint> right, Span<uint> bits, int k)
            {
                Debug.Assert(!right.IsEmpty);
                Debug.Assert(!bits.IsEmpty);
                Debug.Assert(bits.Length + k >= left.Length + right.Length);
 
                // Executes the multiplication algorithm for left and right,
                // but skips the first k limbs of left, which is equivalent to
                // preceding division by 2^(32*k). To spare memory allocations
                // we write the result to an already allocated memory.
 
                if (left.Length > k)
                {
                    left = left.Slice(k);
 
                    if (left.Length < right.Length)
                    {
                        Multiply(right,
                                 left,
                                 bits.Slice(0, left.Length + right.Length));
                    }
                    else
                    {
                        Multiply(left,
                                 right,
                                 bits.Slice(0, left.Length + right.Length));
                    }
 
                    return ActualLength(bits.Slice(0, left.Length + right.Length));
                }
 
                return 0;
            }
 
            private static int SubMod(Span<uint> left, ReadOnlySpan<uint> right, ReadOnlySpan<uint> modulus, int k)
            {
                // 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 (left.Length > k)
                    left = left.Slice(0, k);
                if (right.Length > k)
                    right = right.Slice(0, k);
 
                SubtractSelf(left, right);
                left = left.Slice(0, ActualLength(left));
 
                while (Compare(left, modulus) >= 0)
                {
                    SubtractSelf(left, modulus);
                    left = left.Slice(0, ActualLength(left));
                }
 
                return left.Length;
            }
        }
    }
}