/*
 * 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;
using System.Security.Cryptography;
using System.Text;

/*
 * A custom "big integer" implementation. It internally uses an array
 * of 32-bit integers, that encode the integer in little-endian convention,
 * using two's complement for negative integers.
 *
 * Apart from the array, a single 32-bit field is also present, which
 * encodes the sign. When the value is small (it fits on 32 bits), then
 * the array pointer is null, and the value is in the 32-bit field.
 * Since ZInt is a struct, this means that computations using ZInt do
 * not entail any dynamic (GC-based) memory allocation as long as the
 * value fits on 32 bits. This makes it substantially faster than usual
 * "big integer" implementations (including .NET's implementation, since
 * version 4.0) when values are small.
 *
 * Instances are immutable, and thus can be used as if they were
 * "plain integers".
 *
 * None of this code is "constant-time". As such, ZInt should be
 * considered unsuitable to implementations of cryptographic algorithms.
 */

public struct ZInt : IComparable, IComparable<ZInt>, IEquatable<ZInt> {

	/*
	 * CONVENTIONS:
	 *
	 * If varray == null, then "small" contains the integer value.
	 *
	 * If varray != null, then it contains the value, in little-endian
	 * convention (least significant word comes first) and of
	 * minimal encoded length (i.e. "trimmed"). Two's complement is
	 * used for negative values. "small" is then -1 or 0, depending
	 * on whether the value is negative or not.
	 *
	 * Note that the trimmed value does not always include the sign
	 * bit.
	 *
	 * If the integer value is in the range of values which can be
	 * represented in an "int", then magn is null. There is no allowed
	 * overlap between the two kinds of encodings.
	 *
	 * Default value thus encodes the integer zero.
	 */

	readonly int small;
	readonly uint[] varray;

	/*
	 * The value -1.
	 */
	public static ZInt MinusOne {
		get { return new ZInt(-1, null); }
	}

	/*
	 * The value 0.
	 */
	public static ZInt Zero {
		get { return new ZInt(0, null); }
	}

	/*
	 * The value 1.
	 */
	public static ZInt One {
		get { return new ZInt(1, null); }
	}

	/*
	 * The value 2.
	 */
	public static ZInt Two {
		get { return new ZInt(2, null); }
	}

	/*
	 * Internal constructor which assumes that the provided values are
	 * correct and normalized.
	 */
	private ZInt(int small, uint[] varray)
	{
		this.small = small;
		this.varray = varray;
#if DEBUG
		if (varray != null) {
			if (small != -1 && small != 0) {
				throw new Exception(
					"Bad sign encoding: " + small);
			}
			if (varray.Length == 0) {
				throw new Exception("Empty varray");
			}
			if (Length(small, varray) != varray.Length) {
				throw new Exception("Untrimmed varray");
			}
			if (varray.Length == 1) {
				/*
				 * If there was room for a sign bit, then
				 * the "small" encoding should have been used.
				 */
				if (((varray[0] ^ (uint)small) >> 31) == 0) {
					throw new Exception(
						"suboptimal encoding");
				}
			}
		}
#endif
	}

	/*
	 * Main internal build method. This method normalizes the encoding
	 * ("small" encoding is used when possible; otherwise, the varray
	 * is trimmed and the sign is normalized to -1 or 0).
	 *
	 * If "varray" is null, then the small value is used. Otherwise,
	 * only the sign bit (most significant bit) of "small" is used,
	 * and that value is normalized; the array is trimmed. If the
	 * small encoding then becomes applicable, then it is used.
	 */
	private static ZInt Make(int small, uint[] varray)
	{
		if (varray == null) {
			return new ZInt(small, null);
		}
		small >>= 31;
		int n = Length(small, varray);
		if (n == 1 && (((varray[0] ^ (uint)small) >> 31) == 0)) {
			small = (int)varray[0];
			varray = null;
		} else {
			/*
			 * Note: if n == 0 then the value is -1 or 0, and
			 * "small" already contains the correct value;
			 * Trim() will then return null, which is appropriate.
			 */
			varray = Trim(small, varray, n);
		}
		return new ZInt(small, varray);
	}

	/*
	 * Create an instance from a signed 32-bit integer.
	 */
	public ZInt(int val)
	{
		small = val;
		varray = null;
	}

	/*
	 * Create an instance from an unsigned 32-bit integer.
	 */
	public ZInt(uint val)
	{
		small = (int)val;
		if (small < 0) {
			small = 0;
			varray = new uint[1] { val };
		} else {
			varray = null;
		}
	}

	/*
	 * Create an instance from a signed 64-bit integer.
	 */
	public ZInt(long val)
	{
		small = (int)val;
		if ((long)small == val) {
			varray = null;
		} else {
			ulong uval = (ulong)val;
			uint w0 = (uint)uval;
			uint w1 = (uint)(uval >> 32);
			if (w1 == 0) {
				small = 0;
				varray = new uint[1] { w0 };
			} else if (w1 == 0xFFFFFFFF) {
				small = -1;
				varray = new uint[1] { w0 };
			} else {
				small = (int)w1 >> 31;
				varray = new uint[2] { w0, w1 };
			}
		}
	}

	/*
	 * Create an instance from an unsigned 64-bit integer.
	 */
	public ZInt(ulong val)
	{
		if (val <= 0x7FFFFFFF) {
			small = (int)val;
			varray = null;
		} else {
			small = 0;
			uint w0 = (uint)val;
			uint w1 = (uint)(val >> 32);
			if (w1 == 0) {
				varray = new uint[1] { w0 };
			} else {
				varray = new uint[2] { w0, w1 };
			}
		}
	}

	/*
	 * Create a ZInt instance for an integer expressed as an array
	 * of 32-bit integers (unsigned little-endian convention), with
	 * a specific number of value bits.
	 */
	public static ZInt Make(uint[] words, int numBits)
	{
		if (numBits == 0) {
			return Zero;
		}
		int n = (numBits + 31) >> 5;
		int kb = numBits & 31;
		uint[] m = new uint[n];
		Array.Copy(words, 0, m, 0, n);
		if (kb > 0) {
			m[n - 1] &= ~((uint)0xFFFFFFFF << kb);
		}
		return Make(0, m);
	}

	/*
	 * Create a ZInt instance for an integer expressed as an array
	 * of 64-bit integers (unsigned little-endian convention), with
	 * a specific number of value bits.
	 */
	public static ZInt Make(ulong[] words, int numBits)
	{
		if (numBits == 0) {
			return Zero;
		}
		int kw = (numBits + 63) >> 6;
		int kb = numBits & 63;
		int n = kw * 2;
		uint[] m = new uint[n];
		if (kb != 0) {
			ulong z = words[kw - 1]
				& ~((ulong)0xFFFFFFFFFFFFFFFF << kb);
			m[n - 1] = (uint)(z >> 32);
			m[n - 2] = (uint)z;
			kw --;
			n -= 2;
		}
		for (int i = kw - 1, j = n - 2; i >= 0; i --, j -= 2) {
			ulong z = words[i];
			m[j + 0] = (uint)z;
			m[j + 1] = (uint)(z >> 32);
		}
		return Make(0, m);
	}

	/*
	 * Test whether this value is 0.
	 */
	public bool IsZero {
		get {
			return small == 0 && varray == null;
		}
	}

	/*
	 * Test whether this value is 1.
	 */
	public bool IsOne {
		get {
			return small == 1;
		}
	}

	/*
	 * Test whether this value is even.
	 */
	public bool IsEven {
		get {
			uint w = (varray == null) ? (uint)small : varray[0];
			return (w & 1) == 0;
		}
	}

	/*
	 * Test whether this value is a power of 2. Note that 1 is a power
	 * of 2 (but 0 is not). For negative values, false is returned.
	 */
	public bool IsPowerOfTwo {
		get {
			if (small < 0) {
				return false;
			}
			if (varray == null) {
				return small != 0 && (small & -small) == small;
			}
			int n = varray.Length;
			int z = (int)varray[n - 1];
			if ((z & -z) != z) {
				return false;
			}
			for (int i = n - 2; i >= 0; i --) {
				if (varray[i] != 0) {
					return false;
				}
			}
			return true;
		}
	}

	/*
	 * Get the sign of this value as an integer (-1 for negative
	 * values, 1 for positive, 0 for zero).
	 */
	public int Sign {
		get {
			if (varray == null) {
				if (small < 0) {
					return -1;
				} else if (small == 0) {
					return 0;
				} else {
					return 1;
				}
			} else {
				return small | 1;
			}
		}
	}

	/*
	 * Test whether this value would fit in an 'int'.
	 */
	public bool IsInt {
		get {
			return varray == null;
		}
	}

	/*
	 * Test whether this value would fit in a "long".
	 */
	public bool IsLong {
		get {
			if (varray == null) {
				return true;
			}
			int n = varray.Length;
			if (n == 1) {
				return true;
			} else if (n == 2) {
				return ((int)varray[1] >> 31) == small;
			} else {
				return false;
			}
		}
	}

	/*
	 * Get the length of the value in bits. This is the minimal number
	 * of bits of the two's complement representation of the value,
	 * excluding the sign bit; thus, both 0 and -1 have bit length 0.
	 */
	public int BitLength {
		get {
			if (varray == null) {
				if (small < 0) {
					return 32 - LeadingZeros(~(uint)small);
				} else {
					return 32 - LeadingZeros((uint)small);
				}
			} else {
				int n = varray.Length;
				int bl = n << 5;
				uint w = varray[n - 1];
				if (small < 0) {
					return bl - LeadingZeros(~w);
				} else {
					return bl - LeadingZeros(w);
				}
			}
		}
	}

	/*
	 * Test whether the specified bit has value 1 or 0 ("true" is
	 * returned if the bit has value 1). Note that for negative values,
	 * two's complement representation is assumed.
	 */
	public bool TestBit(int n)
	{
		if (n < 0) {
			throw new ArgumentOutOfRangeException();
		}
		if (varray == null) {
			if (n >= 32) {
				return small < 0;
			} else {
				return (((uint)small >> n) & (uint)1) != 0;
			}
		} else {
			int nw = n >> 5;
			if (nw >= varray.Length) {
				return small != 0;
			}
			int nb = n & 31;
			return ((varray[nw] >> nb) & (uint)1) != 0;
		}
	}

	/*
	 * Copy some bits from this instance to the provided array. Bits
	 * are copied in little-endian order. First bit to be copied is
	 * bit at index "off", and exactly "num" bits are copied. This
	 * method modifies only the minimum number of destination words
	 * (i.e. the first "(num+31)/32" words, exactly). Remaining bits
	 * in the last touched word are set to 0.
	 */
	public void CopyBits(int off, int num, uint[] dest)
	{
		CopyBits(off, num, dest, 0);
	}

	public void CopyBits(int off, int num, uint[] dest, int destOff)
	{
		if (off < 0 || num < 0) {
			throw new ArgumentOutOfRangeException();
		}
		if (num == 0) {
			return;
		}
		ZInt x = this;
		if (off > 0) {
			x >>= off;
		}
		int kw = num >> 5;
		int kb = num & 31;
		uint hmask = ~((uint)0xFFFFFFFF << kb);
		if (x.varray == null) {
			if (kw == 0) {
				dest[destOff] = (uint)x.small & hmask;
			} else {
				uint iw = (uint)(x.small >> 31);
				dest[destOff] = (uint)x.small;
				for (int i = 1; i < kw; i ++) {
					dest[destOff + i] = iw;
				}
				if (kb > 0) {
					dest[destOff + kw] = iw & hmask;
				}
			}
		} else {
			int n = x.varray.Length;
			if (kw <= n) {
				Array.Copy(x.varray, 0, dest, destOff, kw);
			} else {
				Array.Copy(x.varray, 0, dest, destOff, n);
				for (int i = n; i < kw; i ++) {
					dest[destOff + i] = (uint)x.small;
				}
			}
			if (kb > 0) {
				uint last;
				if (kw < n) {
					last = x.varray[kw] & hmask;
				} else {
					last = (uint)x.small & hmask;
				}
				dest[destOff + kw] = last;
			}
		}
	}

	/*
	 * Copy some bits from this instance to the provided array. Bits
	 * are copied in little-endian order. First bit to be copied is
	 * bit at index "off", and exactly "num" bits are copied. This
	 * method modifies only the minimum number of destination words
	 * (i.e. the first "(num+63)/64" words, exactly). Remaining bits
	 * in the last touched word are set to 0.
	 */
	public void CopyBits(int off, int num, ulong[] dest)
	{
		CopyBits(off, num, dest, 0);
	}

	public void CopyBits(int off, int num, ulong[] dest, int destOff)
	{
		if (off < 0 || num < 0) {
			throw new ArgumentOutOfRangeException();
		}
		if (num == 0) {
			return;
		}
		ZInt x = this;
		if (off > 0) {
			x >>= off;
		}
		int kw = num >> 6;
		int kb = num & 63;
		ulong hmask = ~((ulong)0xFFFFFFFFFFFFFFFF << kb);
		long xs = (long)x.small;
		if (x.varray == null) {
			if (kw == 0) {
				dest[destOff] = (ulong)xs & hmask;
			} else {
				ulong iw = (ulong)(xs >> 31);
				dest[destOff] = (ulong)xs;
				for (int i = 1; i < kw; i ++) {
					dest[destOff + i] = iw;
				}
				if (kb > 0) {
					dest[destOff + kw] = iw & hmask;
				}
			}
		} else {
			int n = x.varray.Length;
			uint iw = (uint)x.small;
			int j = 0;
			for (int i = 0; i < kw; i ++, j += 2) {
				uint w0 = (j < n) ? x.varray[j] : iw;
				uint w1 = ((j + 1) < n) ? x.varray[j + 1] : iw;
				dest[destOff + i] =
					(ulong)w0 | ((ulong)w1 << 32);
			}
			if (kb > 0) {
				uint w0 = (j < n) ? x.varray[j] : iw;
				uint w1 = ((j + 1) < n) ? x.varray[j + 1] : iw;
				ulong last = (ulong)w0 | ((ulong)w1 << 32);
				dest[destOff + kw] = last & hmask;
			}
		}
	}

	/*
	 * Extract a 32-bit word at a given offset (counted in bits).
	 * This function is equivalent to right-shifting the value by
	 * "off" bits, then returning the low 32 bits (however, this
	 * function may be more efficient).
	 */
	public uint GetWord(int off)
	{
		if (off < 0) {
			throw new ArgumentOutOfRangeException();
		}
		if (varray == null) {
			int x = small;
			if (off >= 32) {
				off = 31;
			}
			return (uint)(x >> off);
		}
		int n = varray.Length;
		int kw = off >> 5;
		if (kw >= n) {
			return (uint)small;
		}
		int kb = off & 31;
		if (kb == 0) {
			return varray[kw];
		} else {
			uint hi;
			if (kw == n - 1) {
				hi = (uint)small;
			} else {
				hi = varray[kw + 1];
			}
			uint lo = varray[kw];
			return (lo >> kb) | (hi << (32 - kb));
		}
	}

	/*
	 * Extract a 64-bit word at a given offset (counted in bits).
	 * This function is equivalent to right-shifting the value by
	 * "off" bits, then returning the low 64 bits (however, this
	 * function may be more efficient).
	 */
	public ulong GetWord64(int off)
	{
		if (off < 0) {
			throw new ArgumentOutOfRangeException();
		}
		if (varray == null) {
			int x = small;
			if (off >= 32) {
				off = 31;
			}
			return (ulong)(x >> off);
		}
		int n = varray.Length;
		int kw = off >> 5;
		if (kw >= n) {
			return (ulong)small;
		}
		int kb = off & 31;
		if (kb == 0) {
			if (kw == (n - 1)) {
				return (ulong)varray[kw]
					| ((ulong)small << 32);
			} else {
				return (ulong)varray[kw]
					| ((ulong)varray[kw + 1] << 32);
			}
		} else {
			uint v0, v1, v2;
			if (kw == (n - 1)) {
				v0 = varray[kw];
				v1 = (uint)small;
				v2 = (uint)small;
			} else if (kw == (n - 2)) {
				v0 = varray[kw];
				v1 = varray[kw + 1];
				v2 = (uint)small;
			} else {
				v0 = varray[kw];
				v1 = varray[kw + 1];
				v2 = varray[kw + 2];
			}
			uint lo = (v0 >> kb) | (v1 << (32 - kb));
			uint hi = (v1 >> kb) | (v2 << (32 - kb));
			return (ulong)lo | ((ulong)hi << 32);
		}
	}

	/*
	 * Convert this value to an 'int', using silent truncation if
	 * the value does not fit.
	 */
	public int ToInt {
		get {
			return (varray == null) ? small : (int)varray[0];
		}
	}

	/*
	 * Convert this value to an 'uint', using silent truncation if
	 * the value does not fit.
	 */
	public uint ToUInt {
		get {
			return (varray == null) ? (uint)small : varray[0];
		}
	}

	/*
	 * Convert this value to a 'long', using silent truncation if
	 * the value does not fit.
	 */
	public long ToLong {
		get {
			return (long)ToULong;
		}
	}

	/*
	 * Convert this value to an 'ulong', using silent truncation if
	 * the value does not fit.
	 */
	public ulong ToULong {
		get {
			if (varray == null) {
				return (ulong)small;
			} else if (varray.Length == 1) {
				uint iw = (uint)small;
				return (ulong)varray[0] | ((ulong)iw << 32);
			} else {
				return (ulong)varray[0]
					| ((ulong)varray[1] << 32);
			}
		}
	}

	/*
	 * Get the actual length of a varray encoding: this is the minimal
	 * length, in words, needed to encode the value. The value sign
	 * is provided as a negative or non-negative integer, and the
	 * encoding of minimal length does not necessarily include a
	 * sign bit. The value 0 is returned when the array encodes 0
	 * or -1 (depending on sign).
	 */
	static int Length(int sign, uint[] m)
	{
		if (m == null) {
			return 0;
		}
		uint sw = (uint)(sign >> 31);
		int n = m.Length;
		while (n > 0 && m[n - 1] == sw) {
			n --;
		}
		return n;
	}

	/*
	 * Trim an encoding to its minimal encoded length. If the provided
	 * array is already of minimal length, it is returned unchanged.
	 */
	static uint[] Trim(int sign, uint[] m)
	{
		int n = Length(sign, m);
		if (n == 0) {
			return null;
		} else if (n < m.Length) {
			uint[] mt = new uint[n];
			Array.Copy(m, 0, mt, 0, n);
			return mt;
		} else {
			return m;
		}
	}

	/*
	 * Trim or extend a value to the provided length. The returned
	 * array will have the length specified as "n" (if n == 0, then
	 * null is returned). If the source array already has the right
	 * size, then it is returned unchanged.
	 */
	static uint[] Trim(int sign, uint[] m, int n)
	{
		if (n == 0) {
			return null;
		} else if (m == null) {
			m = new uint[n];
			if (sign < 0) {
				Fill(0xFFFFFFFF, m);
			}
			return m;
		}
		int ct = m.Length;
		if (ct < n) {
			uint[] r = new uint[n];
			Array.Copy(m, 0, r, 0, ct);
			return r;
		} else if (ct == n) {
			return m;
		} else {
			uint[] r = new uint[n];
			Array.Copy(m, 0, r, 0, n);
			if (sign < 0) {
				Fill(0xFFFFFFFF, r, ct, n - ct);
			}
			return r;
		}
	}

	static void Fill(uint val, uint[] buf)
	{
		Fill(val, buf, 0, buf.Length);
	}

	static void Fill(uint val, uint[] buf, int off, int len)
	{
		while (len -- > 0) {
			buf[off ++] = val;
		}
	}

	// =================================================================
	/*
	 * Utility methods.
	 *
	 * The methods whose name begins with "Mutate" modify the array
	 * they are given as first parameter; the other methods instantiate
	 * a new array.
	 *
	 * As a rule, untrimmed arrays are accepted as input, and output
	 * may be untrimmed as well.
	 */

	/*
	 * Count the number of leading zeros for a 32-bit value (number of
	 * consecutive zeros, starting with the most significant bit). This
	 * value is between 0 (for a value equal to 2^31 or greater) and
	 * 32 (for zero).
	 */
	static int LeadingZeros(uint v)
	{
		if (v == 0) {
			return 32;
		}
		int n = 0;
		if (v > 0xFFFF) { v >>= 16; } else { n += 16; }
		if (v > 0x00FF) { v >>=  8; } else { n +=  8; }
		if (v > 0x000F) { v >>=  4; } else { n +=  4; }
		if (v > 0x0003) { v >>=  2; } else { n +=  2; }
		if (v <= 0x0001) { n ++; }
		return n;
	}

	/*
	 * Duplicate the provided magnitude array. No attempt is made at
	 * trimming. The source array MUST NOT be null.
	 */
	static uint[] Dup(uint[] m)
	{
		uint[] r = new uint[m.Length];
		Array.Copy(m, 0, r, 0, m.Length);
		return r;
	}

	/*
	 * Increment the provided array. If there is a resulting carry,
	 * then "true" is returned, "false" otherwise. The array MUST
	 * NOT be null.
	 */
	static bool MutateIncr(uint[] x)
	{
		int n = x.Length;
		for (int i = 0; i < n; i ++) {
			uint w = x[i] + 1;
			x[i] = w;
			if (w != 0) {
				return false;
			}
		}
		return true;
	}

	/*
	 * Decrement the provided array. If there is a resulting carry,
	 * then "true" is returned, "false" otherwise. The array MUST
	 * NOT be null.
	 */
	static bool MutateDecr(uint[] x)
	{
		int n = x.Length;
		for (int i = 0; i < n; i ++) {
			uint w = x[i];
			x[i] = w - 1;
			if (w != 0) {
				return false;
			}
		}
		return true;
	}

	/*
	 * Multiply a[] with b[] (unsigned multiplication).
	 */
	static uint[] Mul(uint[] a, uint[] b)
	{
		// TODO: use Karatsuba when operands are large.
		int na = Length(0, a);
		int nb = Length(0, b);
		if (na == 0 || nb == 0) {
			return null;
		}
		uint[] r = new uint[na + nb];
		for (int i = 0; i < na; i ++) {
			ulong ma = a[i];
			ulong carry = 0;
			for (int j = 0; j < nb; j ++) {
				ulong mb = (ulong)b[j];
				ulong mr = (ulong)r[i + j];
				ulong w = ma * mb + mr + carry;
				r[i + j] = (uint)w;
				carry = w >> 32;
			}
			r[i + nb] = (uint)carry;
		}
		return r;
	}

	/*
	 * Get the sign and magnitude of an integer. The sign is
	 * normalized to -1 (negative) or 0 (positive or zero). The
	 * magnitude is an array of length at least 1, containing the
	 * absolute value of this integer; if possible, the varray
	 * is reused (hence, the magnitude array MUST NOT be altered).
	 */
	static void ToAbs(ZInt x, out int sign, out uint[] magn)
	{
		if (x.small < 0) {
			sign = -1;
			x = -x;
		} else {
			sign = 0;
		}
		magn = x.varray;
		if (magn == null) {
			magn = new uint[1] { (uint)x.small }; 
		}
	}

	/*
	 * Compare two integers, yielding -1, 0 or 1.
	 */
	static int Compare(int a, int b)
	{
		if (a < b) {
			return -1;
		} else if (a == b) {
			return 0;
		} else {
			return 1;
		}
	}

	/*
	 * Compare a[] with b[] (unsigned). Returned value is 1 if a[]
	 * is greater than b[], 0 if they are equal, -1 otherwise.
	 */
	static int Compare(uint[] a, uint[] b)
	{
		int ka = Length(0, a);
		int kb = Length(0, b);
		if (ka < kb) {
			return -1;
		} else if (ka == kb) {
			while (ka > 0) {
				ka --;
				uint wa = a[ka];
				uint wb = b[ka];
				if (wa < wb) {
					return -1;
				} else if (wa > wb) {
					return 1;
				}
			}
			return 0;
		} else {
			return 1;
		}
	}

	/*
	 * Add b[] to a[] (unsigned). a[] is modified "in place". Only
	 * n words of a[] are modified. Moreover, the value of
	 * b[] which is added is left-shifted: words b[0]...b[n-1] are
	 * added to a[k]...a[k+n-1]. The final carry is returned ("true"
	 * for 1, "false" for 0). Neither a nor b may be null.
	 */
	static bool MutateAdd(uint[] a, int n, uint[] b, int k)
	{
		bool carry = false;
		for (int i = 0; i < n; i ++) {
			uint wa = a[i + k];
			uint wb = b[i];
			uint wc = wa + wb;
			if (carry) {
				wc ++;
				carry = wa >= wc;
			} else {
				carry = wa > wc;
			}
			a[i + k] = wc;
		}
		return carry;
	}

	/*
	 * Substract b[] from a[] (unsigned). a[] is modified "in
	 * place". Only n words of a[] are modified. Words
	 * b[0]...b[n-1] are subtracted from a[k]...a[k+n-1]. The final
	 * carry is returned ("true" for -1, "false" for 0). Neither a
	 * nor b may be null.
	 */
	static bool MutateSub(uint[] a, int n, uint[] b, int k)
	{
		bool carry = false;
		for (int i = 0; i < n; i ++) {
			uint wa = a[i + k];
			uint wb = b[i];
			uint wc = wa - wb;
			if (carry) {
				wc --;
				carry = wa <= wc;
			} else {
				carry = wa < wc;
			}
			a[i + k] = wc;
		}
		return carry;
	}

	/*
	 * Get the length (in words) of the result of a left-shift of
	 * the provided integer by k bits. If k < 0, then the value is
	 * computed for a right-shift by -k bits.
	 */
	static int GetLengthForLeftShift(ZInt x, int k)
	{
		if (k < 0) {
			if (k == Int32.MinValue) {
				return 0;
			}
			return GetLengthForRightShift(x, -k);
		}
		uint bl = (uint)x.BitLength + (uint)k;
		return (int)((bl + 31) >> 5);
	}

	/*
	 * Get the length (in words) of the result of a right-shift of
	 * the provided integer by k bits. If k < 0, then the value is
	 * computed for a left-shift by -k bits.
	 */
	static int GetLengthForRightShift(ZInt x, int k)
	{
		if (k < 0) {
			if (k == Int32.MinValue) {
				throw new OverflowException();
			}
			return GetLengthForLeftShift(x, -k);
		}
		uint bl = (uint)x.BitLength;
		if (bl <= (uint)k) {
			return 0;
		} else {
			return (int)((bl - k + 31) >> 5);
		}
	}

	/*
	 * Left-shift a[] (unsigned) by k bits. If k < 0, then this becomes
	 * a right-shift.
	 */
	static uint[] ShiftLeft(uint[] a, int k)
	{
		if (k < 0) {
			return ShiftRight(a, -k);
		} else if (k == 0) {
			return a;
		}
		int n = Length(0, a);
		if (n == 0) {
			return null;
		}

		/*
		 * Allocate the result array, with the exact proper size.
		 */
		int bl = ((n << 5) - LeadingZeros(a[n - 1])) + k;
		uint[] r = new uint[(bl + 31) >> 5];

		int kb = k & 31;
		int kw = k >> 5;

		/*
		 * Special case: shift by an integral amount of words.
		 */
		if (kb == 0) {
			Array.Copy(a, 0, r, kw, n);
			return r;
		}

		/*
		 * Copy the bits. This loop handles one source word at
		 * a time, and writes one destination word at a time.
		 * Some unhandled bits may remain at the end.
		 */
		uint bits = 0;
		int zkb = 32 - kb;
		for (int i = 0; i < n; i ++) {
			uint w = a[i];
			r[i + kw] = bits | (w << kb);
			bits = w >> zkb;
		}
		if (bits != 0) {
			r[n + kw] = bits;
		}
		return r;
	}

	/*
	 * Right-shift a[] by k bits. If k < 0, then this becomes
	 * a left-shift.
	 */
	static uint[] ShiftRight(uint[] a, int k)
	{
		if (k < 0) {
			return ShiftLeft(a, -k);
		} else if (k == 0) {
			return a;
		}
		int n = Length(0, a);
		if (n == 0) {
			return null;
		}
		int bl = (n << 5) - LeadingZeros(a[n - 1]) - k;
		if (bl <= 0) {
			return null;
		}
		uint[] r = new uint[(bl + 31) >> 5];

		int kb = k & 31;
		int kw = k >> 5;

		/*
		 * Special case: shift by an integral amount of words.
		 */
		if (kb == 0) {
			Array.Copy(a, kw, r, 0, r.Length);
			return r;
		}

		/*
		 * Copy the bits. This loop handles one source word at
		 * a time, and writes one destination word at a time.
		 * Some unhandled bits may remain at the end.
		 */
		uint bits = a[kw] >> kb;
		int zkb = 32 - kb;
		for (int i = kw + 1; i < n; i ++) {
			uint w = a[i];
			r[i - kw - 1] = bits | (w << zkb);
			bits = w >> kb;
		}
		if (bits != 0) {
			r[n - kw - 1] = bits;
		}
		return r;
	}

	/*
	 * Euclidian division of a[] (unsigned) by b (single word). This
	 * method assumes that b is not 0.
	 */
	static void DivRem(uint[] a, uint b, out uint[] q, out uint r)
	{
		int n = Length(0, a);
		if (n == 0) {
			q = null;
			r = 0;
			return;
		}
		q = new uint[n];
		ulong carry = 0;
		for (int i = n - 1; i >= 0; i --) {
			/*
			 * Performance: we strongly hope that the JIT
			 * compiler will notice that the "/" and "%"
			 * can be combined into a single operation.
			 * TODO: test whether we should replace the
			 * carry computation with:
			 *  carry = w - (ulong)q[i] * b;
			 */
			ulong w = (ulong)a[i] + (carry << 32);
			q[i] = (uint)(w / b);
			carry = w % b;
		}
		r = (uint)carry;
	}

	/*
	 * Euclidian division of a[] (unsigned) by b (single word). This
	 * method assumes that b is not 0. a[] is modified in place. The
	 * remainder (in the 0..b-1 range) is returned.
	 */
	static uint MutateDivRem(uint[] a, uint b)
	{
		int n = Length(0, a);
		if (n == 0) {
			return 0;
		}
		ulong carry = 0;
		for (int i = n - 1; i >= 0; i --) {
			/*
			 * Performance: we strongly hope that the JIT
			 * compiler will notice that the "/" and "%"
			 * can be combined into a single operation.
			 * TODO: test whether we should replace the
			 * carry computation with:
			 *  carry = w - (ulong)q[i] * b;
			 */
			ulong w = (ulong)a[i] + (carry << 32);
			a[i] = (uint)(w / b);
			carry = w % b;
		}
		return (uint)carry;
	}

	/*
	 * Euclidian division of a[] by b[] (unsigned). This method
	 * assumes that b[] is neither 0 or 1, and a[] is not smaller
	 * than b[] (this implies that the quotient won't be zero).
	 */
	static void DivRem(uint[] a, uint[] b, out uint[] q, out uint[] r)
	{
		int nb = Length(0, b);

		/*
		 * Special case when the divisor fits on one word.
		 */
		if (nb == 1) {
			r = new uint[1];
			DivRem(a, b[0], out q, out r[0]);
			return;
		}

		/*
		 * General case.
		 *
		 * We first normalize divisor and dividend such that the
		 * most significant bit of the most significant word of
		 * the divisor is set. We can then compute the quotient
		 * word by word. In details:
		 *
		 * Let:
		 *   w = 2^32 (one word)
		 *   a = (w*a0 + a1) * w^N + a2
		 *   b = b0 * w^N + b2
		 * such that:
		 *   0 <= a0 < w
		 *   0 <= a1 < w
		 *   0 <= a2 < w^N
		 *   w/2 <= b0 < w
		 *   0 <= b2 < w^N
		 *   a < w*b
		 * In other words, a0 and a1 are the two upper words of a[],
		 * b0 is the upper word of b[] and has length 32 bits exactly,
		 * and the quotient of a by b fits in one word.
		 *
		 * Under these conditions, define q and r such that:
		 *   a = b * q + r
		 *   q >= 0
		 *   0 <= r < b
		 * We can then compute a value u this way:
		 *   if a0 = b0, then let u = w-1
		 *   otherwise, let u be such that
		 *     (w*a0 + a1) = u*b0 + v, where 0 <= v < b0
		 * It can then be shown that all these inequations hold:
		 *   0 <= u < w
		 *   u-2 <= q <= u
		 *
		 * In plain words, this allows us to compute an almost-exact
		 * estimate of the upper word of the quotient, with only
		 * one 64-bit division.
		 */

		/*
		 * Normalize dividend and divisor. The normalized dividend
		 * will become the temporary array for the remainder, and
		 * we will modify it, so we make sure we have a copy.
		 */
		int norm = LeadingZeros(b[nb - 1]);
		r = ShiftLeft(a, norm);
		if (r == a) {
			r = new uint[a.Length];
			Array.Copy(a, 0, r, 0, a.Length);
		}
		uint[] b2 = ShiftLeft(b, norm);
		int nr = Length(0, r);
#if DEBUG
		if (Length(0, b2) != nb) {
			throw new Exception("normalize error 1");
		}
		if (b2[nb - 1] < 0x80000000) {
			throw new Exception("normalize error 2");
		}
		{
			uint[] ta = ShiftRight(r, norm);
			if (Compare(a, ta) != 0) {
				throw new Exception("normalize error 3");
			}
			uint[] tb = ShiftRight(b2, norm);
			if (Compare(b, tb) != 0) {
				throw new Exception("normalize error 4");
			}
		}
#endif
		b = b2;

		/*
		 * Length of the quotient will be (at most) k words. This
		 * is the number of iterations in the loop below.
		 */
		int k = (nr - nb) + 1;
#if DEBUG
		if (k <= 0) {
			throw new Exception("wrong iteration count: " + k);
		}
#endif
		q = new uint[k];

		/*
		 * The upper word of a[] (the one we modify, i.e. currently
		 * stored in r[]) is in a0; it is carried over from the
		 * previous loop iteration. Initially, it is zero.
		 */
		uint a0 = 0;
		uint b0 = b[nb - 1];
		int j = nr;
		while (k -- > 0) {
			uint a1 = r[-- j];
			uint u;
			if (a0 == b0) {
				u = 0xFFFFFFFF;
			} else {
				ulong ah = ((ulong)a0 << 32) | (ulong)a1;
				u = (uint)(ah / b0);
			}

			/*
			 * Candidate word for the quotient:
			 * -- if u = 0 then qw is necessarily 0, and the
			 *    rest of this iteration is trivial;
			 * -- if u = 1 then we try qw = 1;
			 * -- otherwise, we try qw = u-1.
			 */
			uint qw;
			if (u == 0) {
				q[k] = 0;
				a0 = a1;
				continue;
			} else if (u == 1) {
				qw = 1;
			} else {
				qw = u - 1;
			}

			/*
			 * "Trying" a candidate word means subtracting from
			 * r[] the product qw*b (b[] being shifted by k words).
			 * The result may be negative, in which case we
			 * overestimated qw; it may be greater than b or
			 * equal to b, in which case we underestimated qw;
			 * or it may be just fine.
			 */
			ulong carry = 0;
			bool tooBig = true;
			for (int i = 0; i < nb; i ++) {
				uint wb = b[i];
				ulong z = (ulong)wb * (ulong)qw + carry;
				carry = z >> 32;
				uint wa = r[i + k];
				uint wc = wa - (uint)z;
				if (wc > wa) {
					carry ++;
				}
				r[i + k] = wc;
				if (wc != wb) {
					tooBig = wc > wb;
				}
			}

			/*
			 * Once we have adjusted everything, the upper word
			 * of r[] will be nullified; wo do it now. Note that
			 * for the first loop iteration, that upper word
			 * may be absent (so already zero, but also "virtual").
			 */
			if (nb + k < nr) {
				r[nb + k] = 0;
			}

			/*
			 * At that point, "carry" should be equal to a0
			 * if we estimated right.
			 */
			if (carry < (ulong)a0) {
				/*
				 * Underestimate.
				 */
				qw ++;
#if DEBUG
				if (carry + 1 != (ulong)a0) {
					throw new Exception("div error 1");
				}
				if (!MutateSub(r, nb, b, k)) {
					throw new Exception("div error 2");
				}
#else
				MutateSub(r, nb, b, k);
#endif
			} else if (carry > (ulong)a0) {
				/*
				 * Overestimate.
				 */
				qw --;
#if DEBUG
				if (carry - 1 != (ulong)a0) {
					throw new Exception("div error 3");
				}
				if (!MutateAdd(r, nb, b, k)) {
					throw new Exception("div error 4");
				}
#else
				MutateAdd(r, nb, b, k);
#endif
			} else if (tooBig) {
				/*
				 * Underestimate, but no expected carry.
				 */
				qw ++;
#if DEBUG
				if (MutateSub(r, nb, b, k)) {
					throw new Exception("div error 5");
				}
#else
				MutateSub(r, nb, b, k);
#endif
			}

			q[k] = qw;
			a0 = r[j];
		}

		/*
		 * At that point, r[] contains the remainder but needs
		 * to be shifted back, to account for the normalization
		 * performed before. q[] is correct (but possibly
		 * untrimmed).
		 */
		r = ShiftRight(r, norm);
	}

	// =================================================================

	/*
	 * Conversion to sbyte; an OverflowException is thrown if the
	 * value does not fit. Use ToInt to get a truncating conversion.
	 */
	public static explicit operator sbyte(ZInt val)
	{
		int x = (int)val;
		if (x < SByte.MinValue || x > SByte.MaxValue) {
			throw new OverflowException();
		}
		return (sbyte)x;
	}

	/*
	 * Conversion to byte; an OverflowException is thrown if the
	 * value does not fit. Use ToInt to get a truncating conversion.
	 */
	public static explicit operator byte(ZInt val)
	{
		int x = (int)val;
		if (x > Byte.MaxValue) {
			throw new OverflowException();
		}
		return (byte)x;
	}

	/*
	 * Conversion to short; an OverflowException is thrown if the
	 * value does not fit. Use ToInt to get a truncating conversion.
	 */
	public static explicit operator short(ZInt val)
	{
		int x = (int)val;
		if (x < Int16.MinValue || x > Int16.MaxValue) {
			throw new OverflowException();
		}
		return (short)x;
	}

	/*
	 * Conversion to ushort; an OverflowException is thrown if the
	 * value does not fit. Use ToInt to get a truncating conversion.
	 */
	public static explicit operator ushort(ZInt val)
	{
		int x = (int)val;
		if (x > UInt16.MaxValue) {
			throw new OverflowException();
		}
		return (ushort)x;
	}

	/*
	 * Conversion to int; an OverflowException is thrown if the
	 * value does not fit. Use ToInt to get a truncating conversion.
	 */
	public static explicit operator int(ZInt val)
	{
		if (val.varray != null) {
			throw new OverflowException();
		}
		return val.small;
	}

	/*
	 * Conversion to uint; an OverflowException is thrown if the
	 * value does not fit. Use ToUInt to get a truncating conversion.
	 */
	public static explicit operator uint(ZInt val)
	{
		int s = val.small;
		if (s < 0) {
			throw new OverflowException();
		}
		uint[] m = val.varray;
		if (m == null) {
			return (uint)s;
		} else if (m.Length > 1) {
			throw new OverflowException();
		} else {
			return m[0];
		}
	}

	/*
	 * Conversion to long; an OverflowException is thrown if the
	 * value does not fit. Use ToLong to get a truncating conversion.
	 */
	public static explicit operator long(ZInt val)
	{
		int s = val.small;
		uint[] m = val.varray;
		if (m == null) {
			return (long)s;
		} else if (m.Length == 1) {
			return (long)m[0] | ((long)s << 32);
		} else if (m.Length == 2) {
			uint w0 = m[0];
			uint w1 = m[1];
			if (((w1 ^ (uint)s) >> 31) != 0) {
				throw new OverflowException();
			}
			return (long)w0 | ((long)w1 << 32);
		} else {
			throw new OverflowException();
		}
	}

	/*
	 * Conversion to ulong; an OverflowException is thrown if the
	 * value does not fit. Use ToULong to get a truncating conversion.
	 */
	public static explicit operator ulong(ZInt val)
	{
		int s = val.small;
		if (s < 0) {
			throw new OverflowException();
		}
		uint[] m = val.varray;
		if (m == null) {
			return (ulong)s;
		} else if (m.Length == 1) {
			return (ulong)m[0];
		} else if (m.Length == 2) {
			return (ulong)m[0] | ((ulong)m[1] << 32);
		} else {
			throw new OverflowException();
		}
	}

	/*
	 * By definition, conversion from sbyte conserves the value.
	 */
	public static implicit operator ZInt(sbyte val)
	{
		return new ZInt((int)val);
	}

	/*
	 * By definition, conversion from byte conserves the value.
	 */
	public static implicit operator ZInt(byte val)
	{
		return new ZInt((uint)val);
	}

	/*
	 * By definition, conversion from short conserves the value.
	 */
	public static implicit operator ZInt(short val)
	{
		return new ZInt((int)val);
	}

	/*
	 * By definition, conversion from ushort conserves the value.
	 */
	public static implicit operator ZInt(ushort val)
	{
		return new ZInt((uint)val);
	}

	/*
	 * By definition, conversion from int conserves the value.
	 */
	public static implicit operator ZInt(int val)
	{
		return new ZInt(val);
	}

	/*
	 * By definition, conversion from uint conserves the value.
	 */
	public static implicit operator ZInt(uint val)
	{
		return new ZInt(val);
	}

	/*
	 * By definition, conversion from long conserves the value.
	 */
	public static implicit operator ZInt(long val)
	{
		return new ZInt(val);
	}

	/*
	 * By definition, conversion from ulong conserves the value.
	 */
	public static implicit operator ZInt(ulong val)
	{
		return new ZInt(val);
	}

	/*
	 * Unary '+' operator is a no-operation.
	 */
	public static ZInt operator +(ZInt a)
	{
		return a;
	}

	/*
	 * Unary negation.
	 */
	public static ZInt operator -(ZInt a)
	{
		int s = a.small;
		uint[] m = a.varray;
		if (m == null) {
			if (s == Int32.MinValue) {
				return new ZInt(0, new uint[1] { 0x80000000 });
			} else {
				return new ZInt(-s, null);
			}
		}

		/*
		 * Two's complement: invert all bits, then add 1. The
		 * result array will usually have the same size, but may
		 * be one word longer or one word shorter.
		 */
		int n = Length(s, m);
		uint[] bm = new uint[n];
		for (int i = 0; i < n; i ++) {
			bm[i] = ~m[i];
		}
		if (MutateIncr(bm)) {
			/*
			 * Extra carry. This may happen only if the source
			 * array contained only zeros, which may happen only
			 * for a source value -2^(32*k) for some integer k
			 * (k > 0).
			 */
			bm = new uint[n + 1];
			bm[n] = 1;
			return new ZInt(0, bm);
		} else {
			/*
			 * The resulting array might be too big by one word,
			 * so we must not assume that it is trimmed.
			 */
			return Make((int)~(uint)s, bm);
		}
	}

	/*
	 * Addition operator.
	 */
	public static ZInt operator +(ZInt a, ZInt b)
	{
		int sa = a.small;
		int sb = b.small;
		uint[] ma = a.varray;
		uint[] mb = b.varray;
		if (ma == null) {
			if (sa == 0) {
				return b;
			}
			if (mb == null) {
				if (sb == 0) {
					return a;
				}
				return new ZInt((long)sa + (long)sb);
			}
			ma = new uint[1] { (uint)sa };
			sa >>= 31;
		} else if (mb == null) {
			if (sb == 0) {
				return a;
			}
			mb = new uint[1] { (uint)sb };
			sb >>= 31;
		}
		int na = ma.Length;
		int nb = mb.Length;
		int n = Math.Max(na, nb) + 1;
		uint[] mc = new uint[n];
		ulong carry = 0;
		for (int i = 0; i < n; i ++) {
			uint wa = i < na ? ma[i] : (uint)sa;
			uint wb = i < nb ? mb[i] : (uint)sb;
			ulong z = (ulong)wa + (ulong)wb + carry;
			mc[i] = (uint)z;
			carry = z >> 32;
		}
		return Make((-(int)carry) ^ sa ^ sb, mc);
	}

	/*
	 * Subtraction operator.
	 */
	public static ZInt operator -(ZInt a, ZInt b)
	{
		int sa = a.small;
		int sb = b.small;
		uint[] ma = a.varray;
		uint[] mb = b.varray;
		if (ma == null) {
			if (sa == 0) {
				return -b;
			}
			if (mb == null) {
				if (sb == 0) {
					return a;
				}
				return new ZInt((long)sa - (long)sb);
			}
			ma = new uint[1] { (uint)sa };
			sa >>= 31;
		} else if (mb == null) {
			if (sb == 0) {
				return a;
			}
			mb = new uint[1] { (uint)sb };
			sb >>= 31;
		}
		int na = ma.Length;
		int nb = mb.Length;
		int n = Math.Max(na, nb) + 1;
		uint[] mc = new uint[n];
		long carry = 0;
		for (int i = 0; i < n; i ++) {
			uint wa = i < na ? ma[i] : (uint)sa;
			uint wb = i < nb ? mb[i] : (uint)sb;
			long z = (long)wa - (long)wb + carry;
			mc[i] = (uint)z;
			carry = z >> 32;
		}
		return Make((int)carry ^ sa ^ sb, mc);
	}

	/*
	 * Increment operator.
	 */
	public static ZInt operator ++(ZInt a)
	{
		int s = a.small;
		uint[] ma = a.varray;
		if (ma == null) {
			return new ZInt((long)s + 1);
		}
		uint[] mb = Dup(ma);
		if (MutateIncr(mb)) {
			int n = ma.Length;
			mb = new uint[n + 1];
			mb[n] = 1;
			return new ZInt(0, mb);
		} else {
			return Make(s, mb);
		}
	}

	/*
	 * Decrement operator.
	 */
	public static ZInt operator --(ZInt a)
	{
		int s = a.small;
		uint[] ma = a.varray;
		if (ma == null) {
			return new ZInt((long)s - 1);
		}

		/*
		 * MutateDecr() will report a carry only if the varray
		 * contained only zeros; since this value was not small,
		 * then it must have been negative.
		 */
		uint[] mb = Dup(ma);
		if (MutateDecr(mb)) {
			int n = mb.Length;
			uint[] mc = new uint[n + 1];
			Array.Copy(mb, 0, mc, 0, n);
			mc[n] = 0xFFFFFFFE;
			return new ZInt(-1, mc);
		} else {
			return Make(s, mb);
		}
	}

	/*
	 * Multiplication operator.
	 */
	public static ZInt operator *(ZInt a, ZInt b)
	{
		int sa = a.small;
		int sb = b.small;
		uint[] ma = a.varray;
		uint[] mb = b.varray;

		/*
		 * Special cases:
		 * -- one of the operands is zero
		 * -- both operands are small
		 */
		if (ma == null) {
			if (sa == 0) {
				return Zero;
			}
			if (mb == null) {
				if (sb == 0) {
					return Zero;
				}
				return new ZInt((long)sa * (long)sb);
			}
		} else if (mb == null) {
			if (sb == 0) {
				return Zero;
			}
		}

		/*
		 * Get both values in sign+magnitude representation.
		 */
		ToAbs(a, out sa, out ma);
		ToAbs(b, out sb, out mb);

		/*
		 * Compute the product. Set the sign.
		 */
		ZInt r = Make(0, Mul(ma, mb));
		if ((sa ^ sb) < 0) {
			r = -r;
		}
		return r;
	}

	/*
	 * Integer division: the quotient is returned, and the remainder
	 * is written in 'r'.
	 *
	 * This operation follows the C# rules for integer division:
	 * -- rounding is towards 0
	 * -- quotient is positive if dividend and divisor have the same
	 *    sign, negative otherwise
	 * -- remainder has the sign of the dividend
	 *
	 * Attempt at dividing by zero triggers a DivideByZeroException.
	 */
	public static ZInt DivRem(ZInt a, ZInt b, out ZInt r)
	{
		ZInt q;
		DivRem(a, b, out q, out r);
		return q;
	}

	static void DivRem(ZInt a, ZInt b,
		out ZInt q, out ZInt r)
	{
		int sa = a.small;
		int sb = b.small;
		uint[] ma = a.varray;
		uint[] mb = b.varray;

		/*
		 * Division by zero triggers an exception.
		 */
		if (mb == null && sb == 0) {
			throw new DivideByZeroException();
		}

		/*
		 * If the dividend is zero, then both quotient and
		 * remainder are zero.
		 */
		if (ma == null && sa == 0) {
			q = Zero;
			r = Zero;
			return;
		}

		/*
		 * If both dividend and divisor are small, then we
		 * use the native 64-bit integer types. If only the
		 * divisor is small, then we have a special fast case
		 * for division by 1 or -1.
		 */
		if (ma == null && mb == null) {
			q = new ZInt((long)sa / (long)sb);
			r = new ZInt((long)sa % (long)sb);
			return;
		}
		if (mb == null) {
			if (sb == 1) {
				q = a;
				r = Zero;
				return;
			} else if (sb == -1) {
				q = -a;
				r = Zero;
				return;
			}
		}

		/*
		 * We know that the dividend is not 0, and the divisor
		 * is not -1, 0 or 1. We now want the sign+magnitude
		 * representations of both operands.
		 */
		ToAbs(a, out sa, out ma);
		ToAbs(b, out sb, out mb);

		/*
		 * If the divisor is greater (in magnitude) than the
		 * dividend, then the quotient is zero and the remainder
		 * is equal to the dividend. If the divisor and dividend
		 * are equal in magnitude, then the remainder is zero and
		 * the quotient is 1 if divisor and dividend have the same
		 * sign, -1 otherwise.
		 */
		int cc = Compare(ma, mb);
		if (cc < 0) {
			q = Zero;
			r = a;
			return;
		} else if (cc == 0) {
			q = (sa == sb) ? One : MinusOne;
			r = Zero;
			return;
		}

		/*
		 * At that point, we know that the divisor is not -1, 0
		 * or 1, and that the quotient will not be 0. We perform
		 * the unsigned division (with the magnitudes), then
		 * we adjust the signs.
		 */

		uint[] mq, mr;
		DivRem(ma, mb, out mq, out mr);

		/*
		 * Quotient is positive if divisor and dividend have the
		 * same sign, negative otherwise. Remainder always has
		 * the sign of the dividend, but it may be zero.
		 */
		q = Make(0, mq);
		if (sa != sb) {
			q = -q;
		}
		r = Make(0, mr);
		if (sa < 0) {
			r = -r;
		}
#if DEBUG
		if (q * b + r != a) {
			throw new Exception("division error");
		}
#endif
	}

	/*
	 * Division operator: see DivRem() for details.
	 */
	public static ZInt operator /(ZInt a, ZInt b)
	{
		ZInt q, r;
		DivRem(a, b, out q, out r);
		return q;
	}

	/*
	 * Remainder operator: see DivRem() for details.
	 */
	public static ZInt operator %(ZInt a, ZInt b)
	{
		ZInt q, r;
		DivRem(a, b, out q, out r);
		return r;
	}

	/*
	 * Reduce this value modulo the provided m. This differs from
	 * '%' in that the returned value is always in the 0 to abs(m)-1
	 * range.
	 */
	public ZInt Mod(ZInt m)
	{
		ZInt r = this % m;
		if (r.small < 0) {
			if (m.small < 0) {
				m = -m;
			}
			r += m;
		}
		return r;
	}

	/*
	 * Left-shift operator.
	 */
	public static ZInt operator <<(ZInt a, int n)
	{
		if (n < 0) {
			if (n == Int32.MinValue) {
				return Zero;
			}
			return a >> (-n);
		}
		int sa = a.small;
		uint[] ma = a.varray;
		if (ma == null) {
			if (n <= 32) {
				return new ZInt((long)sa << n);
			}
			if (sa == 0) {
				return Zero;
			}
		}
		uint[] mr = new uint[GetLengthForLeftShift(a, n)];

		/*
		 * Special case when the shift is a multiple of 32.
		 */
		int kw = n >> 5;
		int kb = n & 31;
		if (kb == 0) {
			if (ma == null) {
				mr[kw] = (uint)sa;
				return new ZInt(sa >> 31, mr);
			} else {
				Array.Copy(ma, 0, mr, kw, ma.Length);
				return new ZInt(sa, mr);
			}
		}

		/*
		 * At that point, we know that the source integer does
		 * not fit in a signed "int", or is shifted by more than
		 * 32 bits, or both. Either way, the result will not fit
		 * in an "int".
		 *
		 * We process all input words one by one.
		 */
		uint rem = 0;
		int ikb = 32 - kb;
		int j;
		if (ma == null) {
			j = 1;
			uint wa = (uint)sa;
			mr[kw] = wa << kb;
			rem = wa >> ikb;
		} else {
			j = ma.Length;
			for (int i = 0; i < j; i ++) {
				uint wa = ma[i];
				mr[i + kw] = rem | (wa << kb);
				rem = wa >> ikb;
			}
		}
		sa >>= 31;
#if DEBUG
		if ((j + kw) == mr.Length - 1) {
			if (rem == ((uint)sa >> ikb)) {
				throw new Exception(
					"Wrong left-shift: untrimmed");
			}
		} else if ((j + kw) == mr.Length) {
			if (rem != ((uint)sa >> ikb)) {
				throw new Exception(
					"Wrong left-shift: dropped bits");
			}
		} else {
			throw new Exception(
				"Wrong left-shift: oversized output length");
		}
#endif
		if ((j + kw) < mr.Length) {
			mr[j + kw] = rem | ((uint)sa << kb);
		}
		return new ZInt(sa, mr);
	}

	/*
	 * Right-shift operator.
	 */
	public static ZInt operator >>(ZInt a, int n)
	{
		if (n < 0) {
			if (n == Int32.MinValue) {
				throw new OverflowException();
			}
			return a << (-n);
		}
		int sa = a.small;
		uint[] ma = a.varray;
		if (ma == null) {
			/*
			 * If right-shifting a "small" value, then we can
			 * do the computation with the native ">>" operator
			 * on "int" values, unless the shift count is 32
			 * or more, in which case we get either 0 or -1,
			 * depending on the source value sign.
			 */
			if (n < 32) {
				return new ZInt(sa >> n, null);
			} else {
				return new ZInt(sa >> 31, null);
			}
		}

		/*
		 * At that point, we know that the source value uses
		 * a non-null varray. We compute the bit length of the
		 * result. If the result would fit in an "int" (bit length
		 * of 31 or less) then we handle it as a special case.
		 */
		int kw = n >> 5;
		int kb = n & 31;
		int bl = a.BitLength - n;
		if (bl <= 0) {
			return new ZInt(sa, null);
		} else if (bl <= 31) {
			if (kb == 0) {
				return new ZInt((int)ma[kw], null);
			} else {
				int p = ma.Length;
				uint w0 = ma[kw];
				uint w1 = (kw + 1) < p ? ma[kw + 1] : (uint)sa;
				return new ZInt((int)((w0 >> kb)
					| (w1 << (32 - kb))), null);
			}
		}

		/*
		 * Result will require an array. Let's allocate it.
		 */
		uint[] mr = new uint[(bl + 31) >> 5];

		/*
		 * Special case when the shift is a multiple of 32.
		 */
		if (kb == 0) {
#if DEBUG
			if (mr.Length != (ma.Length - kw)) {
				throw new Exception(
					"Wrong right-shift: output length");
			}
#endif
			Array.Copy(ma, kw, mr, 0, ma.Length - kw);
			return new ZInt(sa, mr);
		}

		/*
		 * We process all input words one by one.
		 */
		int ikb = 32 - kb;
		uint rem = ma[kw] >> kb;
		int j = ma.Length;
		for (int i = kw + 1; i < j; i ++) {
			uint wa = ma[i];
			mr[i - kw - 1] = rem | (wa << ikb);
			rem = wa >> kb;
		}
#if DEBUG
		if ((j - kw - 1) == mr.Length - 1) {
			if (rem == ((uint)sa >> kb)) {
				throw new Exception(
					"Wrong right-shift: untrimmed");
			}
		} else if ((j - kw - 1) == mr.Length) {
			if (rem != ((uint)sa >> kb)) {
				throw new Exception(
					"Wrong right-shift: dropped bits");
			}
		} else {
			throw new Exception(
				"Wrong right-shift: oversized output length");
		}
#endif
		if ((j - kw - 1) < mr.Length) {
			mr[j - kw - 1] = rem | ((uint)sa << ikb);
		}
		return new ZInt(sa, mr);
	}

	/*
	 * NOTES ON BITWISE BOOLEAN OPERATIONS
	 *
	 * When both operands are "small" then the result is necessarily
	 * small: in "small" encoding, all bits beyond bit 31 of the
	 * two's complement encoding are equal to bit 31, so the result
	 * of computing the operation on bits 31 will also be valid for
	 * all subsequent bits. Therefore, when the two operands are
	 * small, we can just do the operation on the "int" values and the
	 * result is guaranteed "small" as well.
	 */

	/*
	 * Compute the bitwise AND between a "big" and a "small" values.
	 */
	static ZInt AndSmall(int s, uint[] m, int x)
	{
		if (x < 0) {
			uint[] r = Dup(m);
			r[0] &= (uint)x;
			return new ZInt(s, r);
		} else {
			return new ZInt(x & (int)m[0], null);
		}
	}

	/*
	 * Bitwise AND operator.
	 */
	public static ZInt operator &(ZInt a, ZInt b)
	{
		int sa = a.small;
		int sb = b.small;
		uint[] ma = a.varray;
		uint[] mb = b.varray;
		if (ma == null) {
			if (mb == null) {
				return new ZInt(sa & sb, null);
			} else {
				return AndSmall(sb, mb, sa);
			}
		} else if (mb == null) {
			return AndSmall(sa, ma, sb);
		}

		/*
		 * Both values are big. Since upper zero bits force the
		 * result to contain zeros, the result size is that of the
		 * positive operand (the smallest of the two if both are
		 * positive). If both operands are negative, then the
		 * result magnitude may be as large as the largest of the
		 * two source magnitudes.
		 *
		 * Result is negative if and only if both operands are
		 * negative.
		 */
		int na = ma.Length;
		int nb = mb.Length;
		int nr;
		if (sa >= 0) {
			if (sb >= 0) {
				nr = Math.Min(na, nb);
			} else {
				nr = na;
			}
		} else {
			if (sb >= 0) {
				nr = nb;
			} else {
				nr = Math.Max(na, nb);
			}
		}
		uint[] mr = new uint[nr];
		for (int i = 0; i < nr; i ++) {
			uint wa = i < na ? ma[i] : (uint)sa;
			uint wb = i < nb ? mb[i] : (uint)sb;
			mr[i] = wa & wb;
		}
		return Make(sa & sb, mr);
	}

	/*
	 * Compute the bitwise OR between a "big" value (sign and
	 * magnitude, already normalized/trimmed), and a small value
	 * (which can be positive or negative).
	 */
	static ZInt OrSmall(int s, uint[] m, int x)
	{
		if (x < 0) {
			return new ZInt(x | (int)m[0], null);
		} else {
			uint[] r = Dup(m);
			r[0] |= (uint)x;
			return new ZInt(s, r);
		}
	}

	/*
	 * Bitwise OR operator.
	 */
	public static ZInt operator |(ZInt a, ZInt b)
	{
		int sa = a.small;
		int sb = b.small;
		uint[] ma = a.varray;
		uint[] mb = b.varray;
		if (ma == null) {
			if (mb == null) {
				return new ZInt(sa | sb, null);
			} else {
				return OrSmall(sb, mb, sa);
			}
		} else if (mb == null) {
			return OrSmall(sa, ma, sb);
		}

		/*
		 * Both values are big. Since upper one bits force the
		 * result to contain ones, the result size is that of
		 * the negative operand (the greater, i.e. "smallest" of
		 * the two if both are negative). If both operands are
		 * positive, then the result magnitude may be as large
		 * as the largest of the two source magnitudes.
		 *
		 * Result is positive if and only if both operands are
		 * positive.
		 */
		int na = ma.Length;
		int nb = mb.Length;
		int nr;
		if (sa >= 0) {
			if (sb >= 0) {
				nr = Math.Max(na, nb);
			} else {
				nr = nb;
			}
		} else {
			if (sb >= 0) {
				nr = na;
			} else {
				nr = Math.Min(na, nb);
			}
		}
		uint[] mr = new uint[nr];
		for (int i = 0; i < nr; i ++) {
			uint wa = i < na ? ma[i] : (uint)sa;
			uint wb = i < nb ? mb[i] : (uint)sb;
			mr[i] = wa | wb;
		}
		return Make(sa | sb, mr);
	}

	/*
	 * Bitwise XOR operator.
	 */
	public static ZInt operator ^(ZInt a, ZInt b)
	{
		int sa = a.small;
		int sb = b.small;
		uint[] ma = a.varray;
		uint[] mb = b.varray;
		if (ma == null && mb == null) {
			return new ZInt(sa ^ sb, null);
		}
		if (ma == null) {
			int st = sa;
			sa = sb;
			sb = st;
			uint[] mt = ma;
			ma = mb;
			mb = mt;
		}
		if (mb == null) {
			int nx = ma.Length;
			uint[] mx = new uint[nx];
			mx[0] = ma[0] ^ (uint)sb;
			if (nx > 1) {
				if (sb < 0) {
					for (int i = 1; i < nx; i ++) {
						mx[i] = ~ma[i];
					}
				} else {
					Array.Copy(ma, 1, mx, 1, nx - 1);
				}
			}
			return Make(sa ^ (sb >> 31), mx);
		}

		/*
		 * Both operands use varrays.
		 * Result can be as big as the bigger of the two operands
		 * (it _will_ be that big, necessarily, if the two operands
		 * have distinct sizes).
		 */
		int na = ma.Length;
		int nb = mb.Length;
		int nr = Math.Max(na, nb);
		uint[] mr = new uint[nr];
		for (int i = 0; i < nr; i ++) {
			uint wa = i < na ? ma[i] : (uint)sa;
			uint wb = i < nb ? mb[i] : (uint)sb;
			mr[i] = wa ^ wb;
		}
		return Make((sa ^ sb) >> 31, mr);
	}

	/*
	 * Bitwise inversion operator.
	 */
	public static ZInt operator ~(ZInt a)
	{
		int s = a.small;
		uint[] m = a.varray;
		if (m == null) {
			return new ZInt(~s, null);
		} else {
			int n = m.Length;
			uint[] r = new uint[n];
			for (int i = 0; i < n; i ++) {
				r[i] = ~m[i];
			}
			return new ZInt(~s, r);
		}
	}

	/*
	 * Basic comparison; returned value is -1, 0 or 1 depending on
	 * whether this instance is to be considered lower then, equal
	 * to, or greater than the provided object.
	 *
	 * All ZInt instances are considered greater than 'null', and
	 * lower than any non-null object that is not a ZInt.
	 */
	public int CompareTo(object obj)
	{
		if (obj == null) {
			return 1;
		}
		if (!(obj is ZInt)) {
			return -1;
		}
		return CompareTo((ZInt)obj);
	}

	/*
	 * Basic comparison; returned value is -1, 0 or 1 depending on
	 * whether this instance is to be considered lower then, equal
	 * to, or greater than the provided value.
	 */
	public int CompareTo(ZInt v)
	{
		int sv = v.small;
		uint[] mv = v.varray;
		int sign1 = small >> 31;
		int sign2 = sv >> 31;
		if (sign1 != sign2) {
			/*
			 * One of the sign* values is -1, the other is 0.
			 */
			return sign1 - sign2;
		}

		/*
		 * Both values have the same sign. Since the varrays are
		 * trimmed, we can use their presence and length to
		 * quickly resolve most cases.
		 */
		if (small < 0) {
			if (varray == null) {
				if (mv == null) {
					return Compare(small, sv);
				} else {
					return 1;
				}
			} else {
				if (mv == null) {
					return -1;
				} else {
					int n1 = varray.Length;
					int n2 = mv.Length;
					if (n1 < n2) {
						return 1;
					} else if (n1 == n2) {
						return Compare(varray, mv);
					} else {
						return -1;
					}
				}
			}
		} else {
			if (varray == null) {
				if (mv == null) {
					return Compare(small, sv);
				} else {
					return -1;
				}
			} else {
				if (mv == null) {
					return 1;
				} else {
					return Compare(varray, mv);
				}
			}
		}
	}

	/*
	 * Equality comparison: a ZInt instance is equal only to another
	 * ZInt instance that encodes the same integer value.
	 */
	public override bool Equals(object obj)
	{
		if (obj == null) {
			return false;
		}
		if (!(obj is ZInt)) {
			return false;
		}
		return CompareTo((ZInt)obj) == 0;
	}

	/*
	 * Equality comparison: a ZInt instance is equal only to another
	 * ZInt instance that encodes the same integer value.
	 */
	public bool Equals(ZInt v)
	{
		return CompareTo(v) == 0;
	}

	/*
	 * The hash code for a ZInt is equal to its lower 32 bits.
	 */
	public override int GetHashCode()
	{
		if (varray == null) {
			return small;
		} else {
			return (int)varray[0];
		}
	}

	/*
	 * Equality operator.
	 */
	public static bool operator ==(ZInt a, ZInt b)
	{
		return a.CompareTo(b) == 0;
	}

	/*
	 * Inequality operator.
	 */
	public static bool operator !=(ZInt a, ZInt b)
	{
		return a.CompareTo(b) != 0;
	}

	/*
	 * Lower-than operator.
	 */
	public static bool operator <(ZInt a, ZInt b)
	{
		return a.CompareTo(b) < 0;
	}

	/*
	 * Lower-or-equal operator.
	 */
	public static bool operator <=(ZInt a, ZInt b)
	{
		return a.CompareTo(b) <= 0;
	}

	/*
	 * Greater-than operator.
	 */
	public static bool operator >(ZInt a, ZInt b)
	{
		return a.CompareTo(b) > 0;
	}

	/*
	 * Greater-or-equal operator.
	 */
	public static bool operator >=(ZInt a, ZInt b)
	{
		return a.CompareTo(b) >= 0;
	}

	/*
	 * Power function: this raises x to the power e. The exponent e
	 * MUST NOT be negative. If x and e are both zero, then 1 is
	 * returned.
	 */
	public static ZInt Pow(ZInt x, int e)
	{
		if (e < 0) {
			throw new ArgumentOutOfRangeException();
		}
		if (e == 0) {
			return One;
		}
		if (e == 1 || x.IsZero || x.IsOne) {
			return x;
		}
		bool neg = false;
		if (x.Sign < 0) {
			x = -x;
			neg = (e & 1) != 0;
		}
		if (x.IsPowerOfTwo) {
			int t = x.BitLength - 1;
			long u = (long)t * (long)e;
			if (u > (long)Int32.MaxValue) {
				throw new OverflowException();
			}
			x = One << (int)u;
		} else {
			ZInt y = One;
			for (;;) {
				if ((e & 1) != 0) {
					y *= x;
				}
				e >>= 1;
				if (e == 0) {
					break;
				}
				x *= x;
			}
			x = y;
		}
		return neg ? -x : x;
	}

	/*
	 * Modular exponentation: this function raises v to the power e
	 * modulo m. The returned value is reduced modulo m: it will be
	 * in the 0 to abs(m)-1 range.
	 *
	 * The modulus m must be positive. If m is 1, then the result is
	 * 0 (regardless of the values of v and e).
	 *
	 * The exponent e must be nonnegative (this function does not
	 * compute modular inverses). If e is zero, then the result is 1
	 * (except if m is 1).
	 */
	public static ZInt ModPow(ZInt v, ZInt e, ZInt m)
	{
		int se = e.Sign;
		if (se < 0) {
			throw new ArgumentOutOfRangeException();
		}
		int sm = m.Sign;
		if (sm < 0) {
			m = -m;
		} else if (sm == 0) {
			throw new DivideByZeroException();
		}
		if (m.varray == null && m.small == 1) {
			return Zero;
		}
		if (se == 0) {
			return One;
		}
		if (v.IsZero) {
			return Zero;
		}

		// TODO: use Montgomery's multiplication when the exponent
		// is large.
		ZInt x = v.Mod(m);
		for (int n = e.BitLength - 2; n >= 0; n --) {
			x = (x * x).Mod(m);
			if (e.TestBit(n)) {
				x = (x * v).Mod(m);
			}
		}
		return x;
	}

	/*
	 * Get the absolute value of a ZInt.
	 */
	public static ZInt Abs(ZInt x)
	{
		return (x.Sign < 0) ? -x : x;
	}

	private static void AppendHex(StringBuilder sb, uint w, bool trim)
	{
		int i = 28;
		if (trim) {
			for (; i >= 0 && (w >> i) == 0; i -= 4);
			if (i < 0) {
				sb.Append("0");
				return;
			}
		}
		for (; i >= 0; i -= 4) {
			sb.Append("0123456789ABCDEF"[(int)((w >> i) & 0x0F)]);
		}
	}

	/*
	 * Convert this value to hexadecimal. If this instance is zero,
	 * then "0" is returned. Otherwise, the number of digits is
	 * minimal (no leading '0'). A leading '-' is used for negative
	 * values. Hexadecimal digits are uppercase.
	 */
	public string ToHexString()
	{
		if (varray == null && small == 0) {
			return "0";
		}
		StringBuilder sb = new StringBuilder();
		ZInt x = this;
		if (x.small < 0) {
			sb.Append('-');
			x = -x;
		}
		if (x.varray == null) {
			AppendHex(sb, (uint)x.small, true);
		} else {
			int n = x.varray.Length;
			AppendHex(sb, x.varray[n - 1], true);
			for (int j = n - 2; j >= 0; j --) {
				AppendHex(sb, x.varray[j], false);
			}
		}
		return sb.ToString();
	}

	/*
	 * Convert this value to decimal. A leading '-' sign is used for
	 * negative value. The number of digits is minimal (no leading '0',
	 * except for zero, which is returned as "0").
	 */
	public override string ToString()
	{
		return ToString(10);
	}

	private static int[] NDIGITS32 = {
		0, 0, 31, 20, 15, 13, 12, 11, 10, 10, 9, 9, 8, 8, 8, 8, 7, 7,
		7, 7, 7, 7, 7, 7, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6
	};

	private static uint[] RADIXPP32 = {
		0, 0, 2147483648, 3486784401, 1073741824, 1220703125,
		2176782336, 1977326743, 1073741824, 3486784401, 1000000000,
		2357947691, 429981696, 815730721, 1475789056, 2562890625,
		268435456, 410338673, 612220032, 893871739, 1280000000,
		1801088541, 2494357888, 3404825447, 191102976, 244140625,
		308915776, 387420489, 481890304, 594823321, 729000000,
		887503681, 1073741824, 1291467969, 1544804416, 1838265625,
		2176782336
	};

	private static char[] DCHAR =
		"0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ".ToCharArray();

	static void AppendDigits(StringBuilder sb, uint v, int radix, int num)
	{
		while (num -- > 0) {
			sb.Append(DCHAR[v % (uint)radix]);
			v = v / (uint)radix;
		}
	}

	static void AppendDigits(StringBuilder sb, uint v, int radix)
	{
		while (v > 0) {
			sb.Append(DCHAR[v % (uint)radix]);
			v = v / (uint)radix;
		}
	}

	/*
	 * Convert this value to a string, in the provided radix. The
	 * radix must be in the 2 to 36 range (inclusive); uppercase
	 * letters 'A' to 'Z' are used for digits of value 10 to 35.
	 * If the value is zero, then "0" is returned; otherwise, leading
	 * '0' digits are removed. A leading '-' sign is added for negative
	 * values.
	 */
	public string ToString(int radix)
	{
		if (radix < 2 || radix > 36) {
			throw new ArgumentOutOfRangeException();
		}

		/*
		 * Special optimized case for base 16.
		 */
		if (radix == 16) {
			return ToHexString();
		}

		if (IsZero) {
			return "0";
		}

		ZInt x = this;
		if (x.Sign < 0) {
			x = -x;
		}
		StringBuilder sb = new StringBuilder();
		if (x.varray == null) {
			AppendDigits(sb, (uint)x.small, radix);
		} else {
			uint[] m = new uint[x.varray.Length];
			Array.Copy(x.varray, 0, m, 0, m.Length);
			uint prad = RADIXPP32[radix];
			int dnum = NDIGITS32[radix];
			for (;;) {
				uint v = MutateDivRem(m, prad);
				bool qz = Length(0, m) == 0;
				if (qz) {
					AppendDigits(sb, v, radix);
					break;
				} else {
					AppendDigits(sb, v, radix, dnum);
				}
			}
		}
		if (Sign < 0) {
			sb.Append('-');
		}
		return Reverse(sb);
	}

	static string Reverse(StringBuilder sb)
	{
		int n = sb.Length;
		char[] tc = new char[n];
		sb.CopyTo(0, tc, 0, n);
		for (int i = 0, j = n - 1; i < j; i ++, j --) {
			char c = tc[i];
			tc[i] = tc[j];
			tc[j] = c;
		}
		return new string(tc);
	}

	static uint DigitValue(char c, int radix)
	{
		int d;
		if (c >= '0' && c <= '9') {
			d = c - '0';
		} else if (c >= 'a' && c <= 'z') {
			d = (c - 'a') + 10;
		} else if (c >= 'A' && c <= 'Z') {
			d = (c - 'A') + 10;
		} else {
			d = -1;
		}
		if (d < 0 || d >= radix) {
			throw new ArgumentException();
		}
		return (uint)d;
	}

	static ZInt ParseUnsigned(string s, int radix)
	{
		if (s.Length == 0) {
			throw new ArgumentException();
		}
		ZInt x = Zero;
		uint acc = 0;
		int accNum = 0;
		uint prad = RADIXPP32[radix];
		int dnum = NDIGITS32[radix];
		foreach (char c in s) {
			uint d = DigitValue(c, radix);
			acc = acc * (uint)radix + d;
			if (++ accNum == dnum) {
				x = x * (ZInt)prad + (ZInt)acc;
				acc = 0;
				accNum = 0;
			}
		}
		if (accNum > 0) {
			uint p = 1;
			while (accNum -- > 0) {
				p *= (uint)radix;
			}
			x = x * (ZInt)p + (ZInt)acc;
		}
		return x;
	}

	/*
	 * Parse a string:
	 * -- A leading '-' is allowed, to denote a negative value.
	 * -- If there is a "0b" or "0B" header (after any '-' sign),
	 *    then the value is interpreted in base 2.
	 * -- If there is a "0x" or "0X" header (after any '-' sign),
	 *    then the value is interpreted in base 16 (hexadecimal).
	 *    Both uppercase and lowercase letters are accepted.
	 * -- If there is no header, then decimal interpretation is used.
	 *
	 * Unexpected characters (including spaces) trigger exceptions.
	 * There must be at least one digit.
	 */
	public static ZInt Parse(string s)
	{
		s = s.Trim();
		bool neg = false;
		if (s.StartsWith("-")) {
			neg = true;
			s = s.Substring(1);
		}
		int radix;
		if (s.StartsWith("0b") || s.StartsWith("0B")) {
			radix = 2;
			s = s.Substring(2);
		} else if (s.StartsWith("0x") || s.StartsWith("0X")) {
			radix = 16;
			s = s.Substring(2);
		} else {
			radix = 10;
		}
		ZInt x = ParseUnsigned(s, radix);
		return neg ? -x : x;
	}

	/*
	 * Parse a string in the specified radix. The radix must be in
	 * the 2 to 36 range (inclusive). Uppercase and lowercase letters
	 * are accepted for digits in the 10 to 35 range.
	 *
	 * A leading '-' sign is allowed, to denote a negative value.
	 * Otherwise, only digits (acceptable with regards to the radix)
	 * may appear. There must be at least one digit.
	 */
	public static ZInt Parse(string s, int radix)
	{
		if (radix < 2 || radix > 36) {
			throw new ArgumentOutOfRangeException();
		}
		s = s.Trim();
		bool neg = false;
		if (s.StartsWith("-")) {
			neg = true;
			s = s.Substring(1);
		}
		ZInt x = ParseUnsigned(s, radix);
		return neg ? -x : x;
	}

	static uint DecU32BE(byte[] buf, int off, int len)
	{
		switch (len) {
		case 0:
			return 0;
		case 1:
			return buf[off];
		case 2:
			return ((uint)buf[off] << 8)
				| (uint)buf[off + 1];
		case 3:
			return ((uint)buf[off] << 16)
				| ((uint)buf[off + 1] << 8)
				| (uint)buf[off + 2];
		default:
			return ((uint)buf[off] << 24)
				| ((uint)buf[off + 1] << 16)
				| ((uint)buf[off + 2] << 8)
				| (uint)buf[off + 3];
		}
	}

	/*
	 * Decode an integer, assuming unsigned big-endian encoding.
	 * An empty array is decoded as 0.
	 */
	public static ZInt DecodeUnsignedBE(byte[] buf)
	{
		return DecodeUnsignedBE(buf, 0, buf.Length);
	}

	/*
	 * Decode an integer, assuming unsigned big-endian encoding.
	 * An empty array is decoded as 0.
	 */
	public static ZInt DecodeUnsignedBE(byte[] buf, int off, int len)
	{
		while (len > 0 && buf[off] == 0) {
			off ++;
			len --;
		}
		if (len == 0) {
			return Zero;
		} else if (len <= 4) {
			return new ZInt(DecU32BE(buf, off, len));
		}
		uint[] m = new uint[(len + 3) >> 2];
		int i = 0;
		for (int j = len; j > 0; j -= 4) {
			int k = j - 4;
			uint w;
			if (k < 0) {
				w = DecU32BE(buf, off, j);
			} else {
				w = DecU32BE(buf, off + k, 4);
			}
			m[i ++] = w;
		}
		return new ZInt(0, m);
	}

	static RNGCryptoServiceProvider RNG = new RNGCryptoServiceProvider();

	/*
	 * Create a random integer of the provided size. Returned value
	 * is in the 0 (inclusive) to 2^size (exclusive) range. A
	 * cryptographically strong RNG is used to ensure uniform selection.
	 */
	public static ZInt MakeRand(int size)
	{
		if (size <= 0) {
			throw new ArgumentOutOfRangeException();
		}
		byte[] buf = new byte[(size + 7) >> 3];
		RNG.GetBytes(buf);
		int kb = size & 7;
		if (kb != 0) {
			buf[0] &= (byte)(0xFF >> (8 - kb));
		}
		return DecodeUnsignedBE(buf);
	}

	/*
	 * Create a random integer in the 0 (inclusive) to max (exclusive)
	 * range. 'max' must be positive. A cryptographically strong RNG
	 * is used to ensure uniform selection.
	 */
	public static ZInt MakeRand(ZInt max)
	{
		if (max.Sign <= 0) {
			throw new ArgumentOutOfRangeException();
		}
		int bl = max.BitLength;
		for (;;) {
			ZInt x = MakeRand(bl);
			if (x < max) {
				return x;
			}
		}
	}

	/*
	 * Create a random integer in the min (inclusive) to max (exclusive)
	 * range. 'max' must be greater than min. A cryptographically
	 * strong RNG is used to ensure uniform selection.
	 */
	public static ZInt MakeRand(ZInt min, ZInt max)
	{
		if (max <= min) {
			throw new ArgumentOutOfRangeException();
		}
		return min + MakeRand(max - min);
	}

	/*
	 * Check whether this integer is prime. A probabilistic algorithm
	 * is used, that theoretically ensures that a non-prime won't be
	 * declared prime with probability greater than 2^(-100). Note that
	 * this holds regardless of how the integer was generated (this
	 * method does not assume that uniform random selection was used).
	 *
	 * (Realistically, the probability of a computer hardware
	 * malfunction is way greater than 2^(-100), so this property
	 * returns the primality status with as much certainty as can be
	 * achieved with a computer.)
	 */
	public bool IsPrime {
		get {
			return IsProbablePrime(50);
		}
	}

	static uint[] PRIMES_BF = new uint[] {
		0xA08A28AC, 0x28208A20, 0x02088288, 0x800228A2,
		0x20A00A08, 0x80282088, 0x800800A2, 0x08028228,
		0x0A20A082, 0x22880020, 0x28020800, 0x88208082,
		0x02022020, 0x08828028, 0x8008A202, 0x20880880
	};

	private bool IsProbablePrime(int rounds)
	{
		ZInt x = this;
		int cc = x.Sign;
		if (cc == 0) {
			return false;
		} else if (cc < 0) {
			x = -x;
		}
		if (x.varray == null) {
			if (x.small < (PRIMES_BF.Length << 5)) {
				return (PRIMES_BF[x.small >> 5]
					& ((uint)1 << (x.small & 31))) != 0;
			}
		}
		if (!x.TestBit(0)) {
			return false;
		}

		ZInt xm1 = x;
		xm1 --;
		ZInt m = xm1;
		int a;
		for (a = 0; !m.TestBit(a); a ++);
		m >>= a;
		while (rounds -- > 0) {
			ZInt b = MakeRand(Two, x);
			ZInt z = ModPow(b, m, x);
			for (int j = 0; j < a; j ++) {
				if (z == One) {
					if (j > 0) {
						return false;
					}
					break;
				}
				if (z == xm1) {
					break;
				}
				if ((j + 1) < a) {
					z = (z * z) % x;
				} else {
					return false;
				}
			}
		}
		return true;
	}

	/*
	 * Encode this integer as bytes (signed big-endian convention).
	 * Encoding is of minimal length that still contains a sign bit
	 * (compatible with ASN.1 DER encoding).
	 */
	public byte[] ToBytesBE()
	{
		byte[] r = new byte[(BitLength + 8) >> 3];
		ToBytesBE(r, 0, r.Length);
		return r;
	}

	/*
	 * Encode this integer as bytes (signed little-endian convention).
	 * Encoding is of minimal length that still contains a sign bit.
	 */
	public byte[] ToBytesLE()
	{
		byte[] r = new byte[(BitLength + 8) >> 3];
		ToBytesLE(r, 0, r.Length);
		return r;
	}

	/*
	 * Encode this integer as bytes (signed big-endian convention).
	 * Output length is provided; exactly that many bytes will be
	 * written. The value is sign-extended or truncated if needed.
	 */
	public void ToBytesBE(byte[] buf, int off, int len)
	{
		ToBytes(true, buf, off, len);
	}

	/*
	 * Encode this integer as bytes (signed little-endian convention).
	 * Output length is provided; exactly that many bytes will be
	 * written. The value is sign-extended or truncated if needed.
	 */
	public void ToBytesLE(byte[] buf, int off, int len)
	{
		ToBytes(false, buf, off, len);
	}

	/*
	 * Encode this integer as bytes (unsigned big-endian convention).
	 * Encoding is of minimal length, possibly without a sign bit. If
	 * this value is zero, then an empty array is returned. If this
	 * value is negative, then an ArgumentOutOfRangeException is thrown.
	 */
	public byte[] ToBytesUnsignedBE()
	{
		if (Sign < 0) {
			throw new ArgumentOutOfRangeException();
		}
		byte[] r = new byte[(BitLength + 7) >> 3];
		ToBytesBE(r, 0, r.Length);
		return r;
	}

	/*
	 * Encode this integer as bytes (unsigned little-endian convention).
	 * Encoding is of minimal length, possibly without a sign bit. If
	 * this value is zero, then an empty array is returned. If this
	 * value is negative, then an ArgumentOutOfRangeException is thrown.
	 */
	public byte[] ToBytesUnsignedLE()
	{
		if (Sign < 0) {
			throw new ArgumentOutOfRangeException();
		}
		byte[] r = new byte[(BitLength + 7) >> 3];
		ToBytesLE(r, 0, r.Length);
		return r;
	}

	void ToBytes(bool be, byte[] buf, int off, int len)
	{
		uint iw = (uint)small >> 31;
		for (int i = 0; i < len; i ++) {
			int j = i >> 2;
			uint w;
			if (varray == null) {
				w = (j == 0) ? (uint)small : iw;
			} else {
				w = (j < varray.Length) ? varray[j] : iw;
			}
			byte v = (byte)(w >> ((i & 3) << 3));
			if (be) {
				buf[off + len - 1 - i] = v;
			} else {
				buf[off + i] = v;
			}
		}
	}
}