Initial commit.

[BoarSSL] / Crypto / ModInt.cs

/*
 * Copyright (c) 2017 Thomas Pornin <pornin@bolet.org>
 *
 * Permission is hereby granted, free of charge, to any person obtaining 
 * a copy of this software and associated documentation files (the
 * "Software"), to deal in the Software without restriction, including
 * without limitation the rights to use, copy, modify, merge, publish,
 * distribute, sublicense, and/or sell copies of the Software, and to
 * permit persons to whom the Software is furnished to do so, subject to
 * the following conditions:
 *
 * The above copyright notice and this permission notice shall be 
 * included in all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, 
 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND 
 * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
 * BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
 * ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
 * CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 */

using System;

namespace Crypto {

/*
 * Mutable container for a modular integer.
 *
 * Rules:
 *
 *   - Each instance is initialised over a given modulus, or by
 *     duplicating an existing instance.
 *
 *   - All operands for a given operation must share the same modulus,
 *     i.e. the instances must have been created with Dup() calls from
 *     a common ancestor.
 *
 *   - Modulus must be odd and greater than 1.
 */

public class ModInt {

        /*
         * Dedicated class for a modulus. It maintains some useful
         * values that depend only on the modulus.
         */
        class ZMod {

                /*
                 * val = modulus value, 31 bits per word, little-endian
                 * bitLen = modulus bit length (exact)
                 * n0i is such that n0i * val[0] == -1 mod 2^31
                 * R = 2^(31*val.Length) (modular)
                 * R2 = 2^(31*(val.Length*2)) (modular)
                 * Rx[i] = 2^(31*(val.Length+(2^i))) (modular)
                 */
                internal uint[] val;
                internal int bitLen;
                internal uint n0i;
                internal uint[] R;
                internal uint[] R2;
                internal uint[][] Rx;
                internal byte[] vm2;

                internal ZMod(byte[] bmod, int off, int len)
                {
                        bitLen = BigInt.BitLength(bmod, off, len);
                        if (bitLen <= 1 || (bmod[off + len - 1] & 1) == 0) {
                                throw new CryptoException("invalid modulus");
                        }
                        int n = (bitLen + 30) / 31;
                        val = new uint[n];
                        DecodeBE(bmod, off, len, val);
                        uint x = val[0];
                        uint y = 2 - x;
                        y *= 2 - y * x;
                        y *= 2 - y * x;
                        y *= 2 - y * x;
                        y *= 2 - y * x;
                        n0i = (1 + ~y) & 0x7FFFFFFF;

                        /*
                         * Compute the Rx[] values.
                         */
                        int zk = 0;
                        while ((n >> zk) != 0) {
                                zk ++;
                        }
                        R = new uint[n];
                        R[n - 1] = 1;
                        for (int i = 0; i < 31; i ++) {
                                ModMul2(R, val);
                        }
                        Rx = new uint[zk][];
                        Rx[0] = new uint[n];
                        Array.Copy(R, 0, Rx[0], 0, n);
                        for (int i = 0; i < 31; i ++) {
                                ModMul2(Rx[0], val);
                        }
                        for (int k = 1; k < zk; k ++) {
                                Rx[k] = new uint[n];
                                MontyMul(Rx[k - 1], Rx[k - 1], Rx[k], val, n0i);
                        }

                        /*
                         * Compute R2 by multiplying the relevant Rx[]
                         * values.
                         */
                        R2 = null;
                        uint[] tt = new uint[n];
                        for (int k = 0; k < zk; k ++) {
                                if (((n >> k) & 1) == 0) {
                                        continue;
                                }
                                if (R2 == null) {
                                        R2 = new uint[n];
                                        Array.Copy(Rx[k], 0, R2, 0, n);
                                } else {
                                        MontyMul(Rx[k], R2, tt, val, n0i);
                                        Array.Copy(tt, 0, R2, 0, n);
                                }
                        }

                        /*
                         * Compute N-2; used as an exponent for modular
                         * inversion. Since modulus N is odd and greater
                         * than 1, N-2 is necessarily positive.
                         */
                        vm2 = new byte[(bitLen + 7) >> 3];
                        int cc = 2;
                        for (int i = 0; i < vm2.Length; i ++) {
                                int b = bmod[off + len - 1 - i];
                                b -= cc;
                                vm2[vm2.Length - 1 - i] = (byte)b;
                                cc = (b >> 8) & 1;
                        }
                }
        }

        ZMod mod;
        uint[] val;
        uint[] tmp1, tmp2;

        /*
         * Get the modulus bit length.
         */
        public int ModBitLength {
                get {
                        return mod.bitLen;
                }
        }

        /*
         * Test whether this value is zero.
         */
        public bool IsZero {
                get {
                        return IsZeroCT != 0;
                }
        }

        public uint IsZeroCT {
                get {
                        int z = (int)val[0];
                        for (int i = 1; i < val.Length; i ++) {
                                z |= (int)val[i];
                        }
                        return ~(uint)((z | -z) >> 31);
                }
        }

        /*
         * Test whether this value is one.
         */
        public bool IsOne {
                get {
                        return IsOneCT != 0;
                }
        }

        public uint IsOneCT {
                get {
                        int z = (int)val[0] ^ 1;
                        for (int i = 1; i < val.Length; i ++) {
                                z |= (int)val[i];
                        }
                        return ~(uint)((z | -z) >> 31);
                }
        }

        ModInt(ZMod m)
        {
                Init(m);
        }

        /*
         * Create a new instance by decoding the provided modulus
         * (unsigned big-endian). Value is zero.
         */
        public ModInt(byte[] modulus)
                : this(modulus, 0, modulus.Length)
        {
        }

        /*
         * Create a new instance by decoding the provided modulus
         * (unsigned big-endian). Value is zero.
         */
        public ModInt(byte[] modulus, int off, int len)
        {
                Init(new ZMod(modulus, off, len));
        }

        void Init(ZMod mod)
        {
                this.mod = mod;
                int n = mod.val.Length;
                val = new uint[n];
                tmp1 = new uint[n];
                tmp2 = new uint[n];
        }

        /*
         * Duplicate this instance. The new instance uses the same
         * modulus, and its value is initialized to the same value as
         * this instance.
         */
        public ModInt Dup()
        {
                ModInt m = new ModInt(mod);
                Array.Copy(val, 0, m.val, 0, val.Length);
                return m;
        }

        /*
         * Set the value in this instance to a copy of the value in
         * the provided instance. The 'm' instance may use a different
         * modulus, in which case the value may incur modular reduction.
         */
        public void Set(ModInt m)
        {
                /*
                 * If the instances use the same modulus array, then
                 * the value can simply be copied as is.
                 */
                if (mod == m.mod) {
                        Array.Copy(m.val, 0, val, 0, val.Length);
                        return;
                }
                Reduce(m.val);
        }

        /*
         * Set the value in this instance to a copy of the value in
         * the provided instance, but do so only conditionally. If
         * 'ctl' is -1, then the copy is done. If 'ctl' is 0, then the
         * copy is NOT done. This instance and the source instance
         * MUST use the same modulus.
         */
        public void CondCopy(ModInt m, uint ctl)
        {
                if (mod != m.mod) {
                        throw new CryptoException("Not same modulus");
                }
                for (int i = 0; i < val.Length; i ++) {
                        uint w = val[i];
                        val[i] = w ^ (ctl & (m.val[i] ^ w));
                }
        }

        /*
         * Set the value in this instance to a copy of either 'a'
         * (if ctl is -1) or 'b' (if ctl is 0).
         */
        public void CopyMux(uint ctl, ModInt a, ModInt b)
        {
                if (mod != a.mod || mod != b.mod) {
                        throw new CryptoException("Not same modulus");
                }
                for (int i = 0; i < val.Length; i ++) {
                        uint aw = a.val[i];
                        uint bw = b.val[i];
                        val[i] = bw ^ (ctl & (aw ^ bw));
                }
        }

        /*
         * Set the value in this instance to the provided integer value
         * (which MUST be nonnegative).
         */
        public void Set(int x)
        {
                val[0] = (uint)x;
                for (int i = 1; i < val.Length; i ++) {
                        val[i] = 0;
                }
        }

        /*
         * Set this value to either 0 (if ctl is 0) or 1 (if ctl is -1),
         * in Montgomery representation.
         */
        public void SetMonty(uint ctl)
        {
                for (int i = 0; i < val.Length; i ++) {
                        val[i] = ctl & mod.R[i];
                }
        }

        /*
         * Set this value by decoding the provided integer (big-endian
         * unsigned encoding). If the source value does not fit (it is
         * not lower than the modulus), then this method return 0, and
         * this instance is set to 0. Otherwise, it returns -1.
         */
        public uint Decode(byte[] buf)
        {
                return Decode(buf, 0, buf.Length);
        }

        /*
         * Set this value by decoding the provided integer (big-endian
         * unsigned encoding). If the source value does not fit (it is
         * not lower than the modulus), then this method return 0, and
         * this instance is set to 0. Otherwise, it returns -1.
         */
        public uint Decode(byte[] buf, int off, int len)
        {
                /*
                 * Decode into val[]; if the truncation loses some
                 * non-zero bytes, then this returns 0.
                 */
                uint x = DecodeBE(buf, off, len, val);

                /*
                 * Compare with modulus. We want to get a 0 if the
                 * subtraction would not yield a carry, i.e. the
                 * source value is not lower than the modulus.
                 */
                x &= Sub(val, mod.val, 0);

                /*
                 * At that point, x is -1 if the value fits, 0
                 * otherwise.
                 */
                for (int i = 0; i < val.Length; i ++) {
                        val[i] &= x;
                }
                return x;
        }

        /*
         * Set this value by decoding the provided integer (big-endian
         * unsigned encoding). The source value is reduced if necessary.
         */
        public void DecodeReduce(byte[] buf)
        {
                DecodeReduce(buf, 0, buf.Length);
        }

        /*
         * Set this value by decoding the provided integer (big-endian
         * unsigned encoding). The source value is reduced if necessary.
         */
        public void DecodeReduce(byte[] buf, int off, int len)
        {
                uint[] x = new uint[((len << 3) + 30) / 31];
                DecodeBE(buf, off, len, x);
                Reduce(x);
        }

        /*
         * Set the value from the provided source array, with modular
         * reduction.
         */
        void Reduce(uint[] b)
        {
                /*
                 * If the modulus uses only one word then we must use
                 * a special code path.
                 */
                if (mod.val.Length == 1) {
                        ReduceSmallMod(b);
                        return;
                }

                /*
                 * Fast copy of words that do not incur modular
                 * reduction.
                 */
                int aLen = mod.val.Length;
                int bLen = b.Length;
                int cLen = Math.Min(aLen - 1, bLen);
                Array.Copy(b, bLen - cLen, val, 0, cLen);
                for (int i = cLen; i < aLen; i ++) {
                        val[i] = 0;
                }

                /*
                 * Inject extra words. We use the pre-computed
                 * Rx[] values to do shifts.
                 */
                for (int j = bLen - cLen; j > 0;) {
                        /*
                         * We can add power-of-2 words, but less
                         * than the modulus size. Note that the modulus
                         * uses at least two words, so this process works.
                         */
                        int k;
                        for (k = 0;; k ++) {
                                int nk = 1 << (k + 1);
                                if (nk >= aLen || nk > j) {
                                        break;
                                }
                        }
                        int num = 1 << k;
                        MontyMul(val, mod.Rx[k], tmp1, mod.val, mod.n0i);
                        j -= num;
                        Array.Copy(b, j, tmp2, 0, num);
                        for (int i = num; i < tmp2.Length; i ++) {
                                tmp2[i] = 0;
                        }
                        ModAdd(tmp1, tmp2, val, mod.val, 0);
                }
        }

        /*
         * Modular reduction in case the modulus fits on a single
         * word.
         */
        void ReduceSmallMod(uint[] b)
        {
                uint x = 0;
                uint n = mod.val[0];
                int nlen = mod.bitLen;
                uint n0i = mod.n0i;
                uint r2 = mod.R2[0];
                for (int i = b.Length - 1; i >= 0; i --) {
                        /*
                         * Multiply x by R (Montgomery multiplication by R^2).
                         */
                        ulong z = (ulong)x * (ulong)r2;
                        uint u = ((uint)z * n0i) & 0x7FFFFFFF;
                        z += (ulong)u * (ulong)n;
                        x = (uint)(z >> 31);

                        /*
                         * Ensure x fits on 31 bits (it may be up to twice
                         * the modulus at that point). If x >= 2^31 then,
                         * necessarily, x is greater than the modulus, and
                         * the subtraction is sound; moreover, in that case,
                         * subtracting the modulus brings back x to less
                         * than the modulus, hence fitting on 31 bits.
                         */
                        x -= (uint)((int)x >> 31) & n;

                        /*
                         * Add the next word, then reduce. The addition
                         * does not overflow since both operands fit on
                         * 31 bits.
                         *
                         * Since the modulus could be much smaller than
                         * 31 bits, we need a full remainder operation here.
                         */
                        x += b[i];

                        /*
                         * Constant-time modular reduction.
                         * We first perform two subtraction of the
                         * shifted modulus to ensure that the high bit
                         * is cleared. This allows the loop to work
                         * properly.
                         */
                        x -= (uint)((int)x >> 31) & (n << (31 - nlen));
                        x -= (uint)((int)x >> 31) & (n << (31 - nlen));
                        for (int j = 31 - nlen; j >= 0; j --) {
                                x -= (n << j);
                                x += (uint)((int)x >> 31) & (n << j);
                        }
                }
                val[0] = x;
        }

        /*
         * Encode into bytes. Big-endian unsigned encoding is used, the
         * returned array having the minimal length to encode the modulus.
         */
        public byte[] Encode()
        {
                return Encode(false);
        }

        /*
         * Encode into bytes. Big-endian encoding is used; if 'signed' is
         * true, then signed encoding is used: returned value will have a
         * leading bit set to 0. Returned array length is the minimal size
         * for encoding the modulus (with a sign bit if using signed
         * encoding).
         */
        public byte[] Encode(bool signed)
        {
                int x = mod.bitLen;
                if (signed) {
                        x ++;
                }
                byte[] buf = new byte[(x + 7) >> 3];
                Encode(buf, 0, buf.Length);
                return buf;
        }

        /*
         * Encode into bytes. The provided array is fully set; big-endian
         * encoding is used, and extra leading bytes of value 0 are added
         * if necessary. If the destination array is too small, then the
         * value is silently truncated.
         */
        public void Encode(byte[] buf)
        {
                Encode(buf, 0, buf.Length);
        }

        /*
         * Encode into bytes. The provided array chunk is fully set;
         * big-endian encoding is used, and extra leading bytes of value
         * 0 are added if necessary. If the destination array is too
         * small, then the value is silently truncated.
         */
        public void Encode(byte[] buf, int off, int len)
        {
                EncodeBE(val, buf, off, len);
        }

        /*
         * Get the least significant bit of the value (0 or 1).
         */
        public uint GetLSB()
        {
                return val[0] & (uint)1;
        }

        /*
         * Add a small integer to this instance. The small integer
         * 'x' MUST be lower than 2^31 and MUST be lower than the modulus.
         */
        public void Add(uint x)
        {
                tmp1[0] = x;
                for (int i = 1; i < tmp1.Length; i ++) {
                        tmp1[i] = 0;
                }
                ModAdd(val, tmp1, val, mod.val, 0);
        }

        /*
         * Add another value to this instance. The operand 'b' may
         * be the same as this instance.
         */
        public void Add(ModInt b)
        {
                if (mod != b.mod) {
                        throw new CryptoException("Not same modulus");
                }
                ModAdd(val, b.val, val, mod.val, 0);
        }

        /*
         * Subtract a small integer from this instance. The small integer
         * 'x' MUST be lower than 2^31 and MUST be lower than the modulus.
         */
        public void Sub(uint x)
        {
                tmp1[0] = x;
                for (int i = 1; i < tmp1.Length; i ++) {
                        tmp1[i] = 0;
                }
                ModSub(val, tmp1, val, mod.val);
        }

        /*
         * Subtract another value from this instance. The operand 'b'
         * may be the same as this instance.
         */
        public void Sub(ModInt b)
        {
                if (mod != b.mod) {
                        throw new CryptoException("Not same modulus");
                }
                ModSub(val, b.val, val, mod.val);
        }

        /*
         * Negate this value.
         */
        public void Negate()
        {
                ModSub(null, val, val, mod.val);
        }

        /*
         * Convert this instance to Montgomery representation.
         */
        public void ToMonty()
        {
                MontyMul(val, mod.R2, tmp1, mod.val, mod.n0i);
                Array.Copy(tmp1, 0, val, 0, val.Length);
        }

        /*
         * Convert this instance back from Montgomery representation to
         * normal representation.
         */
        public void FromMonty()
        {
                tmp1[0] = 1;
                for (int i = 1; i < tmp1.Length; i ++) {
                        tmp1[i] = 0;
                }
                MontyMul(val, tmp1, tmp2, mod.val, mod.n0i);
                Array.Copy(tmp2, 0, val, 0, val.Length);
        }

        /*
         * Compute a Montgomery multiplication with the provided
         * value. The other operand may be this instance.
         */
        public void MontyMul(ModInt b)
        {
                if (mod != b.mod) {
                        throw new CryptoException("Not same modulus");
                }
                MontyMul(val, b.val, tmp1, mod.val, mod.n0i);
                Array.Copy(tmp1, 0, val, 0, val.Length);
        }

        /*
         * Montgomery-square this instance.
         */
        public void MontySquare()
        {
                MontyMul(val, val, tmp1, mod.val, mod.n0i);
                Array.Copy(tmp1, 0, val, 0, val.Length);
        }

        /*
         * Perform modular exponentiation. Exponent is in big-endian
         * unsigned encoding.
         */
        public void Pow(byte[] exp)
        {
                Pow(exp, 0, exp.Length);
        }

        /*
         * Perform modular exponentiation. Exponent is in big-endian
         * unsigned encoding.
         */
        public void Pow(byte[] exp, int off, int len)
        {
                MontyMul(val, mod.R2, tmp1, mod.val, mod.n0i);
                val[0] = 1;
                for (int i = 1; i < val.Length; i ++) {
                        val[i] = 0;
                }
                for (int i = 0; i < len; i ++) {
                        int x = exp[off + len - 1 - i];
                        for (int j = 0; j < 8; j ++) {
                                MontyMul(val, tmp1, tmp2, mod.val, mod.n0i);
                                uint ctl = (uint)-((x >> j) & 1);
                                for (int k = 0; k < val.Length; k ++) {
                                        val[k] = (tmp2[k] & ctl)
                                                | (val[k] & ~ctl);
                                }
                                MontyMul(tmp1, tmp1, tmp2, mod.val, mod.n0i);
                                Array.Copy(tmp2, 0, tmp1, 0, tmp2.Length);
                        }
                }
        }

        /*
         * Compute modular inverse of this value. If this instance is
         * zero, then it remains equal to zero. If the modulus is not
         * prime, then this function computes wrong values.
         */
        public void Invert()
        {
                Pow(mod.vm2);
        }

        /*
         * Compute the square root for this value. This method assumes
         * that the modulus is prime, greater than or equal to 7, and
         * equal to 3 modulo 4; if it is not, then the returned value
         * and the contents of this instance are indeterminate.
         *
         * Returned value is -1 if the value was indeed a square, 0
         * otherwise. In the latter case, array contents are the square
         * root of the opposite of the original value.
         */
        public uint SqrtBlum()
        {
                /*
                 * We suppose that p = 3 mod 4; we raise to the power
                 * (p-3)/4, then do an extra multiplication to go to
                 * power (p+1)/4.
                 *
                 * Since we know the modulus bit length, we can do
                 * the exponentiation from the high bits downwards.
                 */
                ToMonty();
                Array.Copy(val, 0, tmp1, 0, val.Length);
                int k = (mod.bitLen - 2) / 31;
                int j = mod.bitLen - 2 - k * 31;
                uint ew = mod.val[k];
                for (int i = mod.bitLen - 2; i >= 2; i --) {
                        uint ctl = ~(uint)-((int)(ew >> j) & 1);
                        MontyMul(tmp1, tmp1, tmp2, mod.val, mod.n0i);
                        MontyMul(val, tmp2, tmp1, mod.val, mod.n0i);
                        for (int m = 0; m < tmp1.Length; m ++) {
                                uint w = tmp1[m];
                                tmp1[m] = w ^ (ctl & (w ^ tmp2[m]));
                        }
                        if (-- j < 0) {
                                j = 30;
                                ew = mod.val[-- k];
                        }
                }

                /*
                 * The extra multiplication. Square root is written in
                 * tmp2 (in Montgomery representation).
                 */
                MontyMul(val, tmp1, tmp2, mod.val, mod.n0i);

                /*
                 * Square it back in tmp1, to see if it indeed yields
                 * val.
                 */
                MontyMul(tmp2, tmp2, tmp1, mod.val, mod.n0i);
                int z = 0;
                for (int i = 0; i < val.Length; i ++) {
                        z |= (int)(val[i] ^ tmp1[i]);
                }
                uint good = ~(uint)((z | -z) >> 31);

                /*
                 * Convert back the result to normal representation.
                 */
                Array.Copy(tmp2, 0, val, 0, val.Length);
                FromMonty();
                return good;
        }

        /*
         * Conditionally swap this instance with the one provided in
         * parameter. The swap is performed if ctl is -1, not performed
         * if ctl is 0.
         */
        public void CondSwap(ModInt b, uint ctl)
        {
                if (mod != b.mod) {
                        throw new CryptoException("Not same modulus");
                }
                for (int i = 0; i < val.Length; i ++) {
                        uint x = val[i];
                        uint y = b.val[i];
                        uint m = ctl & (x ^ y);
                        val[i] = x ^ m;
                        b.val[i] = y ^ m;
                }
        }

        /*
         * Compare for equality this value with another. Comparison still
         * works if the two values use distinct moduli.
         */
        public bool Eq(ModInt b)
        {
                return EqCT(b) != 0;
        }

        /*
         * Compare for equality this value with another. Comparison still
         * works if the two values use distinct moduli. Returned value is
         * -1 on equality, 0 otherwise.
         */
        public uint EqCT(ModInt b)
        {
                uint z = 0;
                if (b.val.Length > val.Length) {
                        for (int i = 0; i < val.Length; i ++) {
                                z |= val[i] ^ b.val[i];
                        }
                        for (int i = val.Length; i < b.val.Length; i ++) {
                                z |= b.val[i];
                        }
                } else {
                        for (int i = 0; i < b.val.Length; i ++) {
                                z |= val[i] ^ b.val[i];
                        }
                        for (int i = b.val.Length; i < val.Length; i ++) {
                                z |= b.val[i];
                        }
                }
                int x = (int)z;
                return ~(uint)((x | -x) >> 31);
        }

        /* ============================================================== */

        /*
         * Decode value (unsigned big-endian) into dst[] (little-endian,
         * 31 bits per word). All words of dst[] are initialised.
         * Returned value is -1 on success, 0 if some non-zero source
         * bits had to be ignored.
         */
        static uint DecodeBE(byte[] buf, int off, int len, uint[] dst)
        {
                int i = 0;
                uint acc = 0;
                int accLen = 0;
                off += len;
                while (i < dst.Length) {
                        uint b;
                        if (len > 0) {
                                b = buf[-- off];
                                len --;
                        } else {
                                b = 0;
                        }
                        acc |= (b << accLen);
                        accLen += 8;
                        if (accLen >= 31) {
                                dst[i ++] = acc & 0x7FFFFFFF;
                                accLen -= 31;
                                acc = b >> (8 - accLen);
                        }
                }
                while (len -- > 0) {
                        acc |= buf[-- off];
                }
                int x = (int)acc;
                return ~(uint)((x | -x) >> 31);
        }

        /*
         * Encode an integer (array of words, little-endian, 31 bits per
         * word) into big-endian encoding, with the provided length (in
         * bytes).
         */
        static byte[] EncodeBE(uint[] x, int len)
        {
                byte[] val = new byte[len];
                EncodeBE(x, val, 0, len);
                return val;
        }

        /*
         * Encode an integer (array of words, little-endian, 31 bits per
         * word) into big-endian encoding, with the provided length (in
         * bytes).
         */
        static void EncodeBE(uint[] x, byte[] val)
        {
                EncodeBE(x, val, 0, val.Length);
        }

        /*
         * Encode an integer (array of words, little-endian, 31 bits per
         * word) into big-endian encoding, with the provided length (in
         * bytes).
         */
        static void EncodeBE(uint[] x, byte[] val, int off, int len)
        {
                uint acc = 0;
                int accLen = 0;
                int j = 0;
                for (int i = len - 1; i >= 0; i --) {
                        uint b;
                        if (accLen < 8) {
                                uint z = (j < x.Length) ? x[j ++] : 0;
                                b = acc | (z << accLen);
                                acc = z >> (8 - accLen);
                                accLen += 23;
                        } else {
                                b = acc;
                                accLen -= 8;
                                acc >>= 8;
                        }
                        val[off + i] = (byte)b;
                }
        }

        /*
         * Subtract b from a; the carry is returned (-1 if carry, 0
         * otherwise). The operation is done only if ctl is -1; if
         * ctl is 0, then a[] is unmodified, but the carry is still
         * computed and returned.
         *
         * The two operand arrays MUST have the same size.
         */
        static uint Sub(uint[] a, uint[] b, uint ctl)
        {
                int n = a.Length;
                int cc = 0;
                ctl >>= 1;
                for (int i = 0; i < n; i ++) {
                        uint aw = a[i];
                        uint bw = b[i];
                        uint cw = (uint)cc + aw - bw;
                        cc = (int)cw >> 31;
                        a[i] = (aw & ~ctl) | (cw & ctl);
                }
                return (uint)cc;
        }

        /*
         * Left-shift value by one bit; the extra bit (carry) is
         * return as -1 or 0.
         */
        static uint LShift(uint[] a)
        {
                int n = a.Length;
                uint cc = 0;
                for (int i = 0; i < n; i ++) {
                        uint aw = a[i];
                        a[i] = (cc | (aw << 1)) & 0x7FFFFFFF;
                        cc = aw >> 30;
                }
                return (uint)-(int)cc;
        }

        /*
         * Modular left-shift value by one bit. Value and modulus MUST
         * have the same length. Value MUST be lower than modulus.
         */
        static void ModMul2(uint[] a, uint[] mod)
        {
                int n = a.Length;

                /*
                 * First pass: compute 2*a-mod, but don't keep the
                 * result, only the final carry (0 or -1).
                 */
                uint cc1 = 0;
                int cc2 = 0;
                for (int i = 0; i < n; i ++) {
                        uint aw = a[i];
                        uint aws = ((aw << 1) | cc1) & 0x7FFFFFFF;
                        cc1 = aw >> 30;
                        uint z = aws - mod[i] + (uint)cc2;
                        cc2 = (int)z >> 31;
                }
                cc2 += (int)cc1;

                /*
                 * If cc2 is 0, then the subtraction yields no carry and
                 * must be done. Otherwise, cc2 is -1, and the subtraction
                 * must not be done.
                 */
                uint ctl = ~(uint)cc2;
                cc1 = 0;
                cc2 = 0;
                for (int i = 0; i < n; i ++) {
                        uint aw = a[i];
                        uint aws = ((aw << 1) | cc1) & 0x7FFFFFFF;
                        cc1 = aw >> 30;
                        uint z = aws - (mod[i] & ctl) + (uint)cc2;
                        cc2 = (int)z >> 31;
                        a[i] = z & 0x7FFFFFFF;
                }
        }

        /*
         * Modular addition.
         *
         * If 'hi' is zero, then this computes a+b-n if a+b >= n,
         * a+b otherwise.
         *
         * If 'hi' is non-zero, then this computes 2^k+a+b-n, where
         * k = n.Length*31.
         *
         * Result is written in d[]. The same array may be used in
         * several operands.
         */
        static void ModAdd(uint[] a, uint[] b, uint[] d, uint[] n, uint hi)
        {
                /*
                 * Set ctl to -1 if hi is non-zero, 0 otherwise.
                 * 'ctl' computes whether the subtraction with n[]
                 * is needed in the second pass.
                 */
                int x = (int)hi;
                uint ctl = (uint)((x | -x) >> 31);

                for (int pass = 0; pass < 2; pass ++) {
                        /*
                         * cc1 is the carry for a+b (0 or 1)
                         *
                         * cc2 is the carry for the modulus
                         * subtraction (0 or -1)
                         */
                        uint cc1 = 0;
                        int cc2 = 0;
                        for (int i = 0; i < n.Length; i ++) {
                                uint aw = a[i];
                                uint bw = b[i];
                                uint nw = n[i];
                                uint sw = aw + bw + cc1;
                                cc1 = sw >> 31;
                                sw &= 0x7FFFFFFF;
                                uint dw = sw - nw + (uint)cc2;
                                cc2 = (int)dw >> 31;
                                if (pass == 1) {
                                        dw &= 0x7FFFFFFF;
                                        d[i] = (dw & ctl) | (sw & ~ctl);
                                }
                        }

                        /*
                         * Compute aggregate subtraction carry. This should
                         * not be 1 if the operands are correct (it would
                         * mean that a+b-n overflows, so both a and b are
                         * larger than n). If it is 0, then a+b-n >= 0,
                         * so the subtraction must be done; otherwise, the
                         * aggregate carry will be -1, and the subtraction
                         * should not be done (unless forced by a non-zero
                         * 'hi' value).
                         */
                        cc2 += (int)cc1;
                        ctl |= ~(uint)(cc2 >> 31);
                }
        }

        /*
         * Modular subtraction.
         *
         * This computes a-b, then adds n if the result is negative.
         * If a is null, then it is assumed to be all-zeros.
         *
         * Result is written in d[]. The same array may be used in
         * several operands.
         */
        static void ModSub(uint[] a, uint[] b, uint[] d, uint[] n)
        {
                uint ctl = 0;
                for (int pass = 0; pass < 2; pass ++) {
                        /*
                         * cc1 = carry for a-b (0 or -1)
                         * cc2 = carry for modulus addition (0 or 1)
                         */
                        int cc1 = 0;
                        uint cc2 = 0;
                        for (int i = 0; i < n.Length; i ++) {
                                uint aw = (a == null) ? 0 : a[i];
                                uint bw = b[i];
                                uint nw = n[i];
                                uint sw = aw - bw + (uint)cc1;
                                cc1 = (int)sw >> 31;
                                sw &= 0x7FFFFFFF;
                                uint dw = sw + nw + cc2;
                                cc2 = dw >> 31;
                                if (pass == 1) {
                                        dw &= 0x7FFFFFFF;
                                        d[i] = (dw & ctl) | (sw & ~ctl);
                                }
                        }

                        /*
                         * Modulus addition must be done if and only if
                         * a-b had a final carry.
                         */
                        ctl = (uint)cc1;
                }
        }

        /*
         * Compute Montgomery multiplication of a[] by b[], result in
         * d[], with modulus n[]. All arrays must have the same length.
         * d[] must be distinct from a[], b[] and n[]. Modulus must be
         * odd. n0i must be such that n[0]*n0i = -1 mod 2^31. Values a[]
         * and b[] must be lower than n[] (so, in particular, a[] and
         * b[] cannot be the same array as n[]).
         */
        static void MontyMul(uint[] a, uint[] b, uint[] d, uint[] n, uint n0i)
        {
                int len = n.Length;
                for (int i = 0; i < len; i ++) {
                        d[i] = 0;
                }
                ulong dh = 0;
                for (int i = 0; i < len; i ++) {
                        uint ai = a[i];
                        uint u = ((d[0] + ai * b[0]) * n0i) & 0x7FFFFFFF;
                        ulong cc = 0;
                        for (int j = 0; j < len; j ++) {
                                ulong z = (ulong)d[j]
                                        + (ulong)ai * (ulong)b[j]
                                        + (ulong)u * (ulong)n[j] + cc;
                                cc = z >> 31;
                                if (j > 0) {
                                        d[j - 1] = (uint)z & 0x7FFFFFFF;
                                } else {
                                        // DEBUG
                                        if (((uint)z & 0x7FFFFFFF) != 0) {
                                                throw new Exception("BAD!");
                                        }
                                }
                        }
                        dh += cc;
                        d[len - 1] = (uint)dh & 0x7FFFFFFF;
                        dh >>= 31;
                }
                int x = (int)dh;
                uint ctl = (uint)((x | -x) >> 31);
                Sub(d, n, ctl | ~Sub(d, n, 0));
        }
}

}