ModulusP.cs

using System;
using System.Numerics;
 
namespace Waher.Security.EllipticCurves
{
    public class ModulusP
    {
        protected readonly BigInteger p;
 
        public ModulusP(BigInteger Prime)
        {
            this.p = Prime;
        }
 
        public BigInteger Add(BigInteger a, BigInteger b)
        {
            BigInteger Sum = a + b;
 
            if (Sum >= this.p)
            {
                Sum -= this.p;
                if (Sum >= this.p)
                    Sum %= this.p;
            }
 
            return Sum;
        }
 
        public BigInteger Subtract(BigInteger a, BigInteger b)
        {
            BigInteger Diff = a - b;
            if (Diff.Sign < 0)
            {
                Diff += this.p;
                if (Diff.Sign < 0)
                {
                    Diff %= this.p;
                    if (Diff.Sign < 0)
                        Diff += this.p;
                }
            }
            else if (Diff >= this.p)
            {
                Diff -= this.p;
                if (Diff >= this.p)
                    Diff %= this.p;
            }
 
            return Diff;
        }
 
        public BigInteger Multiply(BigInteger a, BigInteger b)
        {
            return BigInteger.Remainder(a * b, this.p);
        }
 
        public BigInteger Divide(BigInteger a, BigInteger b)
        {
            b = this.Invert(b);
            return BigInteger.Remainder(a * b, this.p);
        }
 
        public BigInteger Negate(BigInteger x)
        {
            return this.p - x;
        }
 
        public BigInteger Invert(BigInteger x)
        {
            if (x.Sign < 0)
            {
                x = BigInteger.Remainder(x, p);
                if (x.Sign < 0)
                    x += p;
            }
            else if (x >= p)
                x = BigInteger.Remainder(x, p);
 
            BigInteger i = this.p;
            BigInteger j = x;
            BigInteger y1 = BigInteger.One;
            BigInteger y2 = BigInteger.Zero;
            BigInteger q, y;
 
            do
            {
                q = BigInteger.DivRem(i, j, out BigInteger r);
                y = y2 - y1 * q;
                i = j;
                j = r;
                y2 = y1;
                y1 = y;
            }
            while (!j.IsZero);
 
            if (!i.IsOne)
                throw new ArgumentException("Number not invertible.", nameof(x));
 
            BigInteger Result = BigInteger.Remainder(y2, this.p);
            if (Result.Sign < 0)
                Result += this.p;
 
            return Result;
        }
 
        public BigInteger Sqrt(BigInteger N)
        {
            return SqrtModP(N, this.p);
        }
 
        public static BigInteger SqrtModP(BigInteger N, BigInteger p)
        {
            // See: https://en.wikipedia.org/wiki/Tonelli–Shanks_algorithm
 
            if (N.Sign < 0)
            {
                N = BigInteger.Remainder(N, p);
                if (N.Sign < 0)
                    N += p;
            }
            else if (N >= p)
                N = BigInteger.Remainder(N, p);
 
            BigInteger pm1d2 = (p - 1) / 2;
            if (BigInteger.ModPow(N, pm1d2, p) != BigInteger.One)
                throw new ArgumentException("No root available.", nameof(N));
 
            BigInteger z = p - 1;
            BigInteger m = z;
            int s = 0;
 
            while (m.IsEven)
            {
                s++;
                m >>= 1;
            }
 
            while (BigInteger.ModPow(z, pm1d2, p) == BigInteger.One)
            {
                z--;
                if (z.IsZero)
                    throw new InvalidOperationException("Nonresidue not found.");
            }
 
            BigInteger c = BigInteger.ModPow(z, m, p);
            BigInteger c2 = BigInteger.Remainder(c * c, p);
            BigInteger u = BigInteger.ModPow(N, m, p);
            BigInteger v = BigInteger.ModPow(N, (m + 1) / 2, p);
 
            while (--s > 0)
            {
                if (BigInteger.ModPow(u, BigInteger.Pow(2, s - 1), p) != BigInteger.One)
                {
                    u = BigInteger.Remainder(u * c2, p);
                    v = BigInteger.Remainder(v * c, p);
                }
 
                c = c2;
                c2 = BigInteger.Remainder(c * c, p);
            }
 
            return v;
        }
 
        public static int CalcBits(BigInteger n)
        {
            if (n.IsZero)
                return 0;
 
            return CalcBits(n.ToByteArray());
        }
 
        public static int CalcBits(byte[] A)
        { 
            int c = A.Length - 1;
            int i = c << 3;
            byte b = A[c];
 
            while (b > 0)
            {
                i++;
                b >>= 1;
            }
 
            return i;
        }
 
    }
}