Ouroboros.Maths.cs

using System.Buffers;
using System.Numerics;
using System.Runtime.CompilerServices;
using System.Runtime.InteropServices;
using System.Security.Cryptography;

namespace Ouroboros.Maths;

public static class BinaryIntegerConstants<T> where T : IBinaryInteger<T>
{
    public static T Eleven { get; }
    public static T Five { get; }
    public static T Four { get; }
    public static T Log2Size { get; }
    public static T Seven { get; }
    public static T Seventeen { get; }
    public static T Six { get; }
    public static int Size { get; }
    public static T Ten { get; }
    public static T Thirteen { get; }
    public static T Three { get; }
    public static T Two { get; }

    static BinaryIntegerConstants() {
        var size = int.CreateChecked(value: T.PopCount(value: T.AllBitsSet));

        Eleven = T.CreateTruncating(value: 11U);
        Five = T.CreateTruncating(value: 5U);
        Four = T.CreateTruncating(value: 4U);
        Log2Size = T.Log2(value: T.CreateChecked(value: size));
        Six = T.CreateTruncating(value: 6U);
        Seven = T.CreateTruncating(value: 7U);
        Seventeen = T.CreateTruncating(value: 17U);
        Size = int.CreateChecked(value: size);
        Ten = T.CreateTruncating(value: 10U);
        Thirteen = T.CreateTruncating(value: 13U);
        Three = T.CreateTruncating(value: 3U);
        Two = T.CreateTruncating(value: 2U);
    }
}
public static class BinaryIntegerFunctions
{
    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    internal static T As<T>(this bool value) where T : IBinaryInteger<T> =>
        T.CreateTruncating(value: Unsafe.As<bool, byte>(source: ref value));
    /// <remarks>
    /// Should get translated to BLSR (or the equivalent) on supported platforms.
    /// </remarks>
    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    internal static T ClearLeastSignificantBit<T>(this T value) where T : IBinaryInteger<T> =>
        (value & (value - T.One));
    /// <remarks>
    /// Should get translated to BLSI (or the equivalent) on supported platforms.
    /// </remarks>
    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    internal static T ExtractLeastSignificantBit<T>(this T value) where T : IBinaryInteger<T> =>
        (value & (-value));
    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    internal static T IsGreaterThan<T>(this T value, T other) where T : IBinaryInteger<T> =>
        (value > other).As<T>();
    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    internal static T IsNonZero<T>(this T value) where T : IBinaryInteger<T> =>
        (T.Zero != value).As<T>();
    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    internal static T IsPowerOfTwo<T>(this T value) where T : IBinaryInteger<T> =>
        (value.ClearLeastSignificantBit().IsZero() - (T.Zero >= value).As<T>());
    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    internal static T IsZero<T>(this T value) where T : IBinaryInteger<T> =>
        (T.Zero == value).As<T>();
    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    internal static T Max<T>(this T value, T other) where T : IBinaryInteger<T> =>
        (value ^ ((value ^ other) & (-other.IsGreaterThan(other: value))));
    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    internal static T NthFermatMask<T>(this int value) where T : IBinaryInteger<T> =>
        (T.AllBitsSet / value.NthFermatNumber<T>());
    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    internal static T NthFermatNumber<T>(this int value) where T : IBinaryInteger<T> =>
        ((T.One << (1 << value)) + T.One);
    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    internal static T SetLeastSignificantBits<T>(this T value) where T : IBinaryInteger<T> =>
        (value | (value - T.One));

    public static T DigitalRoot<T>(this T value) where T : IBinaryInteger<T> {
        var x = value.IsNonZero();
        var y = T.Abs(value: value);
        var z = T.CreateChecked(value: 9U);

        return (x + ((y - x) % z));
    }
    public static IEnumerable<T> EnumerateDigits<T>(this T value) where T : IBinaryInteger<T> {
        var quotient = value;

        do {
            (quotient, var remainder) = T.DivRem(left: quotient, right: BinaryIntegerConstants<T>.Ten);

            yield return remainder;
        } while (T.Zero != quotient);
    }
    public static T GreatestCommonDivisor<T>(this T value, T other) where T : IBinaryInteger<T> {
        if (T.Zero == other) { return value; }
        else if (T.Zero == value) { return other; }

        other = T.Abs(value: other);
        value = T.Abs(value: value);

        var x = int.CreateTruncating(value: T.TrailingZeroCount(value: (other | value)));

        do {
            other >>= int.CreateTruncating(value: T.TrailingZeroCount(value: other));
            value >>= int.CreateTruncating(value: T.TrailingZeroCount(value: value));

            var y = ((other ^ value) & (-(value < other).As<T>()));

            other ^= y;
            value ^= y;
            value -= other;
        } while (T.Zero != value);

        return (other << x);
    }
    public static T LeastCommonMultiple<T>(this T value, T other) where T : IBinaryInteger<T> {
        var gcd = value.GreatestCommonDivisor(other: other);

        return ((value / gcd) * other);
    }
    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    public static T LeastSignificantBit<T>(this T value) where T : IBinaryInteger<T> =>
        (value.IsNonZero() * (T.TrailingZeroCount(value: value) + T.One));
    public static TResult MortonPair<TInput, TResult>(this TInput value, TInput other) where TInput : IBinaryInteger<TInput> where TResult : IBinaryInteger<TResult> {
        var x = TResult.CreateTruncating(value: value);
        var y = TResult.CreateTruncating(value: other);

        if (BinaryIntegerConstants<TResult>.Size > 128) {
            var i = (BitOperations.Log2(value: ((uint)(BinaryIntegerConstants<TResult>.Size))) - 7);

            do {
                x = ((x | (x << (1 << i))) & (i).NthFermatMask<TResult>());
                y = ((y | (y << (1 << i))) & (i).NthFermatMask<TResult>());
            } while (0 < --i);
        }

        if (BinaryIntegerConstants<TResult>.Size > 64) {
            x = ((x | (x << 64)) & (6).NthFermatMask<TResult>());
            y = ((y | (y << 64)) & (6).NthFermatMask<TResult>());
        }

        if (BinaryIntegerConstants<TResult>.Size > 32) {
            x = ((x | (x << 32)) & (5).NthFermatMask<TResult>());
            y = ((y | (y << 32)) & (5).NthFermatMask<TResult>());
        }

        if (BinaryIntegerConstants<TResult>.Size > 16) {
            x = ((x | (x << 16)) & (4).NthFermatMask<TResult>());
            y = ((y | (y << 16)) & (4).NthFermatMask<TResult>());
        }

        if (BinaryIntegerConstants<TResult>.Size > 8) {
            x = ((x | (x << 8)) & (3).NthFermatMask<TResult>());
            y = ((y | (y << 8)) & (3).NthFermatMask<TResult>());
        }

        if (BinaryIntegerConstants<TResult>.Size > 4) {
            x = ((x | (x << 4)) & (2).NthFermatMask<TResult>());
            y = ((y | (y << 4)) & (2).NthFermatMask<TResult>());
        }

        if (BinaryIntegerConstants<TResult>.Size > 2) {
            x = ((x | (x << 2)) & (1).NthFermatMask<TResult>());
            y = ((y | (y << 2)) & (1).NthFermatMask<TResult>());
        }

        x = ((x | (x << 1)) & (0).NthFermatMask<TResult>());
        y = ((y | (y << 1)) & (0).NthFermatMask<TResult>());

        return (x | (y << 1));
    }
    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    public static T MostSignificantBit<T>(this T value) where T : IBinaryInteger<T> =>
        (T.CreateTruncating(value: BinaryIntegerConstants<T>.Size) - T.LeadingZeroCount(value: value));
    public static T NthSquare<T>(this T value) where T : IBinaryInteger<T> =>
        (value * value);
    public static T PermuteBitsLexicographically<T>(this T value) where T : IBinaryInteger<T> {
        var x = value.SetLeastSignificantBits();
        var y = int.CreateTruncating(value: (T.TrailingZeroCount(value: value) + T.One));
        var z = (((~x).ExtractLeastSignificantBit() - T.One) >> y);

        return ((x + T.One) | z);
    }
    public static T PopulationParity<T>(this T value) where T : IBinaryInteger<T> =>
        (T.PopCount(value: value) & T.One);
    public static T ReflectedBinaryDecode<T>(this T value) where T : IBinaryInteger<T> {
        var x = (value ^ (value >> 1));

        if (BinaryIntegerConstants<T>.Size > 2) {
            x ^= (x >> 2);
        }

        if (BinaryIntegerConstants<T>.Size > 4) {
            x ^= (x >> 4);
        }

        if (BinaryIntegerConstants<T>.Size > 8) {
            x ^= (x >> 8);
        }

        if (BinaryIntegerConstants<T>.Size > 16) {
            x ^= (x >> 16);
        }

        if (BinaryIntegerConstants<T>.Size > 32) {
            x ^= (x >> 32);
        }

        if (BinaryIntegerConstants<T>.Size > 64) {
            x ^= (x >> 64);
        }

        if (BinaryIntegerConstants<T>.Size > 128) {
            var i = (BitOperations.Log2(value: ((uint)(BinaryIntegerConstants<T>.Size))) - 7);

            do {
                x ^= (x >> (int.CreateTruncating(value: BinaryIntegerConstants<T>.Size) >> i));
            } while (0 < --i);
        }

        return x;
    }
    public static T ReflectedBinaryEncode<T>(this T value) where T : IBinaryInteger<T> =>
        (value ^ (value >> 1));
}
public sealed class Pcg32XshRr
{
    private const ulong MAGIC_VALUE_DEFAULT = 6364136223846793005UL;
    private const string STREAM_VALUE_ERROR = "stream offset must be a positive integer less than 2^63";
    private const ulong STREAM_VALUE_MAX = ((1UL << 63) - 1UL);

    private static ulong Jump(ulong count, ulong magic, ulong state, ulong stream) {
        unchecked {
            var accMul = 1UL;
            var accAdd = 0UL;
            var curMul = magic;
            var curAdd = stream;

            while (count > 0UL) {
                if (0UL < (count & 1UL)) {
                    accMul *= curMul;
                    accAdd = ((accAdd * curMul) + curAdd);
                }

                curAdd = ((curMul + 1UL) * curAdd);
                curMul *= curMul;
                count >>= 1;
            }

            return ((accMul * state) + accAdd);
        }
    }
    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    private static uint NextUInt32(ulong magic, ref ulong state, ulong stream) {
        unchecked {
            state = ((state * magic) + stream);

            return uint.RotateRight(rotateAmount: ((int)(state >> 59)), value: ((uint)(((state >> 18) ^ state) >> 27)));
        }
    }
    private static uint Sample(uint exclusiveHigh, ulong magic, ref ulong state, ulong stream) {
        unchecked {
            var sample = ((ulong)NextUInt32(magic: magic, state: ref state, stream: stream));
            var result = (sample * exclusiveHigh);
            var leftover = ((uint)result);

            if (leftover < exclusiveHigh) {
                var threshold = ((((uint)(-exclusiveHigh)) % exclusiveHigh));

                while (leftover < threshold) {
                    sample = NextUInt32(magic: magic, state: ref state, stream: stream);
                    result = (sample * exclusiveHigh);
                    leftover = ((uint)result);
                }
            }

            return ((uint)(result >> 32));
        }
    }

    public static Pcg32XshRr New(ulong magic, ulong state, ulong stream) => new(magic: magic, state: state, stream: stream);
    public static Pcg32XshRr New(ulong state, ulong stream) => New(magic: MAGIC_VALUE_DEFAULT, state: state, stream: stream);
    public static Pcg32XshRr New() => New(state: SecureRandom.NextUInt64(), stream: SecureRandom.NextUInt64(maximum: STREAM_VALUE_MAX, minimum: 0UL));

    private readonly ulong m_magic;
    private readonly ulong m_stream;

    private ulong m_state;

    private Pcg32XshRr(ulong magic, ulong state, ulong stream) {
        if (stream > STREAM_VALUE_MAX) {
            throw new ArgumentOutOfRangeException(actualValue: stream, message: STREAM_VALUE_ERROR, paramName: nameof(stream));
        }

        m_magic = magic;
        m_state = state;
        m_stream = ((((~(stream & 1UL)) << 63) | stream) | 1UL);
    }

    public void Jump(ulong count) {
        m_state = Jump(count: count, magic: m_magic, state: m_state, stream: m_stream);
    }
    public uint NextPrime(uint maximum, uint minimum) {
        var result = (NextUInt32(maximum: maximum, minimum: minimum) | 1U);

        while (!result.IsPrime()) {
            result += 2U;
        }

        return result;
    }
    public uint NextUInt32() => NextUInt32(magic: m_magic, state: ref m_state, stream: m_stream);
    public uint NextUInt32(uint maximum, uint minimum) {
        var range = (maximum - minimum);

        return ((range != uint.MaxValue) ? (Sample(exclusiveHigh: (range + 1U), magic: m_magic, state: ref m_state, stream: m_stream) + minimum) : NextUInt32());
    }
}
public static class PrimeExtensions
{
    private static ReadOnlySpan<ulong> MillerRabinBases32 => new ulong[] { 2UL, 7UL, 61UL, };

    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    private static UInt128 GetModuloMultiplier(ulong divisor) =>
        ((UInt128.MaxValue / divisor) + UInt128.One);
    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    private static ulong Modulo(ulong divisor, UInt128 multiplier, ulong value) =>
        ((ulong)(((((multiplier * value) >> (sizeof(ulong) * 8)) + UInt128.One) * divisor) >> (sizeof(ulong) * 8)));

    public static bool IsPrime(this uint value) {
        if (0U == (value & 01U)) { return (02U == value); }
        if (0U == (value % 03U)) { return (03U == value); }
        if (0U == (value % 05U)) { return (05U == value); }
        if (0U == (value % 07U)) { return (07U == value); }
        if (0U == (value % 11U)) { return (11U == value); }
        if (0U == (value % 13U)) { return (13U == value); }
        if (0U == (value % 17U)) { return (17U == value); }
        if (0U == (value % 19U)) { return (19U == value); }
        if (0U == (value % 23U)) { return (23U == value); }
        if (0U == (value % 29U)) { return (29U == value); }
        if (0U == (value % 31U)) { return (31U == value); }
        if (0U == (value % 37U)) { return (37U == value); }
        if (0U == (value % 41U)) { return (41U == value); }
        if (0U == (value % 43U)) { return (43U == value); }
        if (0U == (value % 47U)) { return (47U == value); }
        if (0U == (value % 53U)) { return (53U == value); }
        if (0U == (value % 59U)) { return (59U == value); }
        if (0U == (value % 61U)) { return (61U == value); }
        if (0U == (value % 67U)) { return (67U == value); }
        if (0U == (value % 71U)) { return (71U == value); }
        if (0U == (value % 73U)) { return (73U == value); }
        if (0U == (value % 79U)) { return (79U == value); }
        if (0U == (value % 83U)) { return (83U == value); }
        if (0U == (value % 89U)) { return (89U == value); }
        if (0U == (value % 97U)) { return (97U == value); }

        var exponent = (value - 1U);
        var multiplier = GetModuloMultiplier(divisor: value);
        var shift = ((int)uint.TrailingZeroCount(value: exponent));
        var valueMinusOne = exponent;

        exponent >>= shift;

        for (var i = 0; (i < MillerRabinBases32.Length); ++i) { // IsStrongPseudoprime
            var @base = MillerRabinBases32[i];
            var index = exponent;
            var witness = 1UL;

            do { // ModularExponentiation
                if (0U != (index & 1U)) {
                    witness = Modulo(
                        divisor: value,
                        multiplier: multiplier,
                        value: (@base * witness)
                    );
                }

                @base = Modulo(
                    divisor: value,
                    multiplier: multiplier,
                    value: (@base * @base)
                );
            } while (0U != (index >>= 1));

            if ((1U != witness) && (valueMinusOne != witness)) {
                while (++index < shift) {
                    witness = Modulo(
                        divisor: value,
                        multiplier: multiplier,
                        value: (witness * witness)
                    );

                    if ((1U == witness)) { return false; }
                    if (valueMinusOne == witness) { break; }
                }

                if (index == shift) { return false; }
            }
        }

        return (1U != value);
    }
    public static uint NthPrime(this uint value) {
        if (203280220U < value) { return 0U; }
        if (3U > value) { return ((0U == value) ? 2U : ((1U == value) ? 3U : 5U)); }

        var factor = (((uint)((value * 1.10445d) * Math.Log(d: value))) | 1U);
        var index = ((factor.IsPrime() ? (factor - 1U) : factor).Totient() - 1U);

        do {
            if (factor.IsPrime() && (++index == value)) { break; }

            factor += 2U;
        } while (index <= value);

        return factor;
    }
    /// <remarks>
    /// Based on the source code:
    ///     https://github.com/favre49/ThayirSadamLibrary/blob/main/number-theory/PrimeCount.hpp
    ///
    /// One believes that that maximum memory usage of this function is 512 kilobytes.
    /// Variables w, x, y, and z are used to indicate that one has no idea WTF is going on.
    /// Variables i, j, and k are used as indices to various arrays that one also doesn't really understand.
    /// </remarks>
    public static uint Totient(this uint value) {
        if (value.IsPrime()) { return --value; }
        if (value < 9U) { return ((value < 2U) ? 0U : ((value < 3U) ? 1U : ((value < 5U) ? 2U : ((value < 7U) ? 3U : 4U)))); }

        var squareRoot = ((uint)Math.Sqrt(d: value));
        var squareRootHalved = ((squareRoot + 1U) >> 1);
        var larges = ArrayPool<uint>.Shared.Rent(minimumLength: ((int)squareRootHalved));
        var roughs = ArrayPool<uint>.Shared.Rent(minimumLength: ((int)squareRootHalved));
        var smalls = ArrayPool<uint>.Shared.Rent(minimumLength: ((int)squareRootHalved));

        for (var index = 0U; (index < squareRootHalved); ++index) {
            larges[index] = (((value / ((index << 1) + 1U)) - 1U) >> 1);
            roughs[index] = ((index << 1) + 1U);
            smalls[index] = index;
        }

        var counter = 0U;
        var factor = 3U;

        if (factor <= squareRoot) {
            var ignore = ArrayPool<bool>.Shared.Rent(minimumLength: ((int)(squareRoot + 1U)));

            Array.Clear(array: ignore);

            do {
                if (!ignore[factor]) {
                    var factorSquared = (factor * factor);

                    if ((((ulong)factorSquared) * factorSquared) > value) { break; }

                    ignore[factor] = true;

                    var i = factorSquared;
                    var j = 0U;
                    var k = 0U;

                    while (i <= squareRoot) {
                        ignore[i] = true;
                        i += (factor << 1);
                    }
                    while (j < squareRootHalved) {
                        i = roughs[j];

                        if (!ignore[i]) {
                            var x = (i * factor);
                            var y = ((x > squareRoot) ? smalls[(((value / x) - 1U) >> 1)] : larges[(smalls[(x >> 1)] - counter)]);
                            var z = (larges[j] - y);

                            larges[k] = (z + counter);
                            roughs[k++] = i;
                        }

                        ++j;
                    }

                    i = ((squareRoot - 1U) >> 1);
                    j = (((squareRoot / factor) - 1U) | 1U);
                    squareRootHalved = k;

                    while (j >= factor) {
                        var x = (smalls[(j >> 1)] - counter);

                        for (k = ((j * factor) >> 1); (i >= k); --i) {
                            smalls[i] -= x;
                        }

                        j -= 2U;
                    }

                    ++counter;
                }

                factor += 2U;
            } while (factor <= squareRoot);

            ArrayPool<bool>.Shared.Return(array: ignore);
        }

        larges[0] += (((squareRootHalved + ((counter - 1U) << 1)) * (squareRootHalved - 1U)) >> 1);

        for (var i = 1U; (i < squareRootHalved); ++i) { larges[0] -= larges[i]; }
        for (var i = 1U; (i < squareRootHalved); ++i) {
            var w = roughs[i];
            var x = (value / w);
            var y = (smalls[(((x / w) - 1U) >> 1)] - counter);

            if (y < (i + 1U)) { break; }

            var z = 0U;

            for (var j = (i + 1U); (j <= y); ++j) {
                z += smalls[(((x / roughs[j]) - 1U) >> 1)];
            }

            larges[0] += (z - ((y - i) * ((counter + i) - 1U)));
        }

        value = (larges[0] + 1U);

        ArrayPool<uint>.Shared.Return(array: smalls);
        ArrayPool<uint>.Shared.Return(array: roughs);
        ArrayPool<uint>.Shared.Return(array: larges);

        return value;
    }
}
public static class SecureRandom
{
    private static uint NextUInt32(uint exclusiveHigh) {
        var range = (uint.MaxValue - (((uint.MaxValue % exclusiveHigh) + 1U) % exclusiveHigh));
        uint result;

        do {
            result = NextUInt32();
        } while (result > range);

        return (result % exclusiveHigh);
    }
    private static ulong NextUInt64(ulong exclusiveHigh) {
        var range = (ulong.MaxValue - (((ulong.MaxValue % exclusiveHigh) + 1UL) % exclusiveHigh));
        ulong result;

        do {
            result = NextUInt64();
        } while (result > range);

        return (result % exclusiveHigh);
    }

    public static uint NextPrime() {
        var result = (NextUInt32() | 1U);

        while (!result.IsPrime()) {
            result += 2U;
        }

        return result;
    }
    public static uint NextUInt32() {
        var result = 0U;

        RandomNumberGenerator.Fill(data: MemoryMarshal.AsBytes(span: new Span<uint>(reference: ref result)));

        return result;
    }
    public static uint NextUInt32(uint maximum, uint minimum) {
        var range = (maximum - minimum);

        return ((range != uint.MaxValue) ? (NextUInt32(exclusiveHigh: (range + 1U)) + minimum) : NextUInt32());
    }
    public static ulong NextUInt64() {
        var result = 0UL;

        RandomNumberGenerator.Fill(data: MemoryMarshal.AsBytes(span: new Span<ulong>(reference: ref result)));

        return result;
    }
    public static ulong NextUInt64(ulong maximum, ulong minimum) {
        var range = (maximum - minimum);

        return ((range != ulong.MaxValue) ? (NextUInt64(exclusiveHigh: (range + 1UL)) + minimum) : NextUInt64());
    }
}
public static class UnsignedNumberFunctions
{
    public static TResult ElegantPair<TInput, TResult>(this TInput value, TInput other) where TInput : IBinaryInteger<TInput>, IUnsignedNumber<TInput> where TResult : IBinaryInteger<TResult>, IUnsignedNumber<TResult> {
        var x = value.Max(other: other);
        var y = ((value ^ other) * (x & TInput.One));
        var z = TResult.CreateTruncating(value: x);

        return ((z * (z + TResult.One)) + (TResult.CreateTruncating(value: (y ^ other)) - TResult.CreateTruncating(value: (y ^ value))));
    }
    public static (TResult x, TResult y) ElegantUnpair<TInput, TResult>(this TInput value) where TInput : IBinaryInteger<TInput>, IUnsignedNumber<TInput> where TResult : IBinaryInteger<TResult>, IUnsignedNumber<TResult> {
        var x = value.SquareRoot();
        var y = TResult.CreateTruncating(value: (value - (x * x)));
        var z = TResult.CreateTruncating(value: x);

        if (y < z) {
            (y, z) = (z, y);
        }
        else {
            y = ((z << 1) - y);
        }

        if (TResult.IsOddInteger(value: TResult.Max(x: y, y: z))) {
            (y, z) = (z, y);
        }

        return (y, z);
    }
    public static IEnumerable<T> EnumeratePrimeFactors<T>(this T value) where T : IBinaryInteger<T>, IUnsignedNumber<T> {
        if (BinaryIntegerConstants<T>.Four > value) { yield break; }
        if (BinaryIntegerConstants<T>.Five == value) { yield break; }
        if (BinaryIntegerConstants<T>.Seven == value) { yield break; }
        if (BinaryIntegerConstants<T>.Eleven == value) { yield break; }
        if (BinaryIntegerConstants<T>.Thirteen == value) { yield break; }

        var index = value;

        while (T.Zero == (index & T.One)/* enumerate factors of 2 */) {
            yield return BinaryIntegerConstants<T>.Two;

            index >>= 1;
        }
        while (T.Zero == (index % BinaryIntegerConstants<T>.Three)) { // enumerate factors of 3
            yield return BinaryIntegerConstants<T>.Three;

            index /= BinaryIntegerConstants<T>.Three;
        }
        while (T.Zero == (index % BinaryIntegerConstants<T>.Five)/* enumerate factors of 5 */) {
            yield return BinaryIntegerConstants<T>.Five;

            index /= BinaryIntegerConstants<T>.Five;
        }
        while (T.Zero == (index % BinaryIntegerConstants<T>.Seven)/* enumerate factors of 7 */) {
            yield return BinaryIntegerConstants<T>.Seven;

            index /= BinaryIntegerConstants<T>.Seven;
        }
        while (T.Zero == (index % BinaryIntegerConstants<T>.Eleven)/* enumerate factors of 11 */) {
            yield return BinaryIntegerConstants<T>.Eleven;

            index /= BinaryIntegerConstants<T>.Eleven;
        }
        while (T.Zero == (index % BinaryIntegerConstants<T>.Thirteen)/* enumerate factors of 13 */) {
            yield return BinaryIntegerConstants<T>.Thirteen;

            index /= BinaryIntegerConstants<T>.Thirteen;
        }

        var factor = BinaryIntegerConstants<T>.Seventeen;
        var limit = index.SquareRoot();

        if (factor <= limit) {
            do {
                while (T.Zero == (index % factor)/* enumerate factors of (30k - 13) */) {
                    yield return factor;

                    index /= factor;
                }

                factor += BinaryIntegerConstants<T>.Two;

                while (T.Zero == (index % factor)/* enumerate factors of (30k - 11) */) {
                    yield return factor;

                    index /= factor;
                }

                factor += BinaryIntegerConstants<T>.Four;

                while (T.Zero == (index % factor)/* enumerate factors of (30k - 7) */) {
                    yield return factor;

                    index /= factor;
                }

                factor += BinaryIntegerConstants<T>.Six;

                while (T.Zero == (index % factor)/* enumerate factors of (30k - 1) */) {
                    yield return factor;

                    index /= factor;
                }

                factor += BinaryIntegerConstants<T>.Two;

                while (T.Zero == (index % factor)/* enumerate factors of (30k + 1) */) {
                    yield return factor;

                    index /= factor;
                }

                factor += BinaryIntegerConstants<T>.Six;

                while (T.Zero == (index % factor)/* enumerate factors of (30k + 7) */) {
                    yield return factor;

                    index /= factor;
                }

                factor += BinaryIntegerConstants<T>.Four;

                while (T.Zero == (index % factor)/* enumerate factors of (30k + 11) */) {
                    yield return factor;

                    index /= factor;
                }

                factor += BinaryIntegerConstants<T>.Two;

                while (T.Zero == (index % factor)/* enumerate factors of (30k + 13) */) {
                    yield return factor;

                    index /= factor;
                }

                factor += BinaryIntegerConstants<T>.Four;
                limit = index.SquareRoot();
            } while (factor <= limit);
        }

        if ((index != T.One) && (index != value)) {
            yield return index;
        }
    }
    public static T Log10<T>(this T value) where T : IBinaryInteger<T>, IUnsignedNumber<T> {
        return BinaryIntegerConstants<T>.Size switch {
#if !FORCE_SOFTWARE_LOG10
            8 => (T.CreateTruncating(value: ((uint)MathF.Log10(x: uint.CreateTruncating(value: value)))) + T.One),
            16 => (T.CreateTruncating(value: ((uint)MathF.Log10(x: uint.CreateTruncating(value: value)))) + T.One),
            32 => (T.CreateTruncating(value: ((uint)Math.Log10(d: uint.CreateTruncating(value: value)))) + T.One),
#endif
            _ => SoftwareImplementation(value: value),
        };

        static T SoftwareImplementation(T value) {
            var quotient = value;
            var result = T.Zero;

            do {
                quotient /= BinaryIntegerConstants<T>.Ten;
                ++result;
            } while (T.Zero != quotient);

            return result;
        }
    }
    /// <remarks>
    /// Based on the paper:
    ///     An Improved Integer Multiplicative Inverse(modulo 2w)
    ///     Jeffrey Hurchalla, April 2022
    ///     https://arxiv.org/ftp/arxiv/papers/2204/2204.04342.pdf
    /// </remarks>
    public static T ModularInverse<T>(this T value) where T : IBinaryInteger<T>, IUnsignedNumber<T> {
        var x = ((T.CreateChecked(value: 3) * value) ^ T.CreateChecked(value: 2));
        var y = (T.One - (value * x));

        x *= (y + T.One);

        if (BinaryIntegerConstants<T>.Size > 8) {
            y *= y;
            x *= (y + T.One);
        }

        if (BinaryIntegerConstants<T>.Size > 16) {
            y *= y;
            x *= (y + T.One);
        }

        if (BinaryIntegerConstants<T>.Size > 32) {
            y *= y;
            x *= (y + T.One);
        }

        if (BinaryIntegerConstants<T>.Size > 64) {
            y *= y;
            x *= (y + T.One);
        }

        if (BinaryIntegerConstants<T>.Size > 128) {
            var i = (BitOperations.Log2(value: (((uint)BinaryIntegerConstants<T>.Size) / 4)) - 5);

            do {
                y *= y;
                x *= (y + T.One);
            } while (0 < --i);
        }

        return x;
    }
    public static T NextPowerOfTwo<T>(this T value) where T : IBinaryInteger<T>, IUnsignedNumber<T> {
        var x = (BinaryIntegerConstants<T>.Size - int.CreateTruncating(value: T.LeadingZeroCount(value: (value - T.One))));
        var y = int.CreateTruncating(value: BinaryIntegerConstants<T>.Log2Size);

        return ((T.One ^ T.CreateTruncating(value: (((uint)x) >> y))) << x);
    }
    public static T NextSquare<T>(this T value) where T : IBinaryInteger<T>, IUnsignedNumber<T> {
        var squareRootPlusOne = (value.SquareRoot() + T.One);

        return (squareRootPlusOne * squareRootPlusOne);
    }
    /// <remarks>
    /// Based on the paper:
    ///     Square Rooting Algorithms for Integer and Floating-Point Numbers
    ///     IEEE Transactions on Computers, August 1990
    ///     Volume: 39, Issue: 8, Pages: 1025 - 1029
    /// </remarks>
    public static T SquareRoot<T>(this T value) where T : IBinaryInteger<T>, IUnsignedNumber<T> {
        return BinaryIntegerConstants<T>.Size switch {
#if !FORCE_SOFTWARE_SQRT
            8 => T.CreateTruncating(value: ((uint)MathF.Sqrt(x: uint.CreateTruncating(value: value)))),
            16 => T.CreateTruncating(value: ((uint)MathF.Sqrt(x: uint.CreateTruncating(value: value)))),
            32 => T.CreateTruncating(value: ((uint)Math.Sqrt(d: uint.CreateTruncating(value: value)))),
            64 => T.CreateTruncating(value: Sqrt(value: ulong.CreateTruncating(value: value))),
#endif
            _ => SoftwareImplementation(value: value),
        };

        /*
             Credit goes to njuffa for providing a reference implementation:
                 https://stackoverflow.com/a/31149161/1186165

             Notes:
                 - This implementation of the algorithm runs in constant time, based on the size of T.
                 - Ignoring the loop that is entered when the size of T exceeds 64, all branches get eliminated during JIT compilation.
         */
        static T SoftwareImplementation(T value) {
            var msb = int.CreateTruncating(value: value.MostSignificantBit());
            var msbIsOdd = (msb & 1);
            var m = ((msb + 1) >> 1);
            var mMinusOne = (m - 1);
            var mPlusOne = (m + 1);
            var x = (T.One << mMinusOne);
            var y = (x - (value >> (mPlusOne - msbIsOdd)));
            var z = y;

            x += x;

            if (BinaryIntegerConstants<T>.Size > 8) {
                y = (((y * y) >> mPlusOne) + z);
                y = (((y * y) >> mPlusOne) + z);
            }

            if (BinaryIntegerConstants<T>.Size > 16) {
                y = (((y * y) >> mPlusOne) + z);
                y = (((y * y) >> mPlusOne) + z);
                y = (((y * y) >> mPlusOne) + z);
                y = (((y * y) >> mPlusOne) + z);
            }

            if (BinaryIntegerConstants<T>.Size > 32) {
                y = (((y * y) >> mPlusOne) + z);
                y = (((y * y) >> mPlusOne) + z);
                y = (((y * y) >> mPlusOne) + z);
                y = (((y * y) >> mPlusOne) + z);
                y = (((y * y) >> mPlusOne) + z);
                y = (((y * y) >> mPlusOne) + z);
                y = (((y * y) >> mPlusOne) + z);
                y = (((y * y) >> mPlusOne) + z);
            }

            if (BinaryIntegerConstants<T>.Size > 64) {
                var i = T.CreateTruncating(value: (BinaryIntegerConstants<T>.Size >> 3));

                do {
                    i -= (T.One << 3);
                    y = (((y * y) >> mPlusOne) + z);
                    y = (((y * y) >> mPlusOne) + z);
                    y = (((y * y) >> mPlusOne) + z);
                    y = (((y * y) >> mPlusOne) + z);
                    y = (((y * y) >> mPlusOne) + z);
                    y = (((y * y) >> mPlusOne) + z);
                    y = (((y * y) >> mPlusOne) + z);
                    y = (((y * y) >> mPlusOne) + z);
                } while (i != T.Zero);
            }

            y = (x - y);
            x = T.CreateTruncating(value: msbIsOdd);
            y -= BinaryIntegerConstants<T>.Size switch {
                8 => (x * ((y * T.CreateChecked(value: 5UL)) >> 4)),
                16 => (x * ((y * T.CreateChecked(value: 75UL)) >> 8)),
                32 => (x * ((y * T.CreateChecked(value: 19195UL)) >> 16)),
                64 => (x * ((y * T.CreateChecked(value: 1257966796UL)) >> 32)),
                128 => (x * ((y * T.CreateChecked(value: 5402926248376769403UL)) >> 64)),
                _ => throw new NotSupportedException(), // TODO: Research a way to calculate the proper constant at runtime.
            };
            x = (T.One << (BinaryIntegerConstants<T>.Size - 1));
            y -= (value - (y * y)).IsGreaterThan(other: x);

            if (BinaryIntegerConstants<T>.Size > 8) {
                y -= (value - (y * y)).IsGreaterThan(other: x);
                y -= (value - (y * y)).IsGreaterThan(other: x);
            }

            if (BinaryIntegerConstants<T>.Size > 32) {
                y -= (value - (y * y)).IsGreaterThan(other: x);
                y -= (value - (y * y)).IsGreaterThan(other: x);
                y -= (value - (y * y)).IsGreaterThan(other: x);
            }

            return (y & (T.AllBitsSet >> 1));
        }
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        static uint Sqrt(ulong value) {
            var x = ((uint)Math.Sqrt(d: unchecked((long)value)));
            var y = (unchecked(((ulong)x) * x) > value).As<uint>(); // ((x * x) > value) ? 1 : 0
            var z = ((uint)(value >> 63)); // (64 == value.MostSignificantBit()) ? 1 : 0

            return unchecked(x - (y | z));
        }
    }
}