2 * Copyright (c) 2017 Thomas Pornin <pornin@bolet.org>
4 * Permission is hereby granted, free of charge, to any person obtaining
5 * a copy of this software and associated documentation files (the
6 * "Software"), to deal in the Software without restriction, including
7 * without limitation the rights to use, copy, modify, merge, publish,
8 * distribute, sublicense, and/or sell copies of the Software, and to
9 * permit persons to whom the Software is furnished to do so, subject to
10 * the following conditions:
12 * The above copyright notice and this permission notice shall be
13 * included in all copies or substantial portions of the Software.
15 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
16 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
17 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
18 * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
19 * BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
20 * ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
21 * CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
26 using System.Security.Cryptography;
30 * A custom "big integer" implementation. It internally uses an array
31 * of 32-bit integers, that encode the integer in little-endian convention,
32 * using two's complement for negative integers.
34 * Apart from the array, a single 32-bit field is also present, which
35 * encodes the sign. When the value is small (it fits on 32 bits), then
36 * the array pointer is null, and the value is in the 32-bit field.
37 * Since ZInt is a struct, this means that computations using ZInt do
38 * not entail any dynamic (GC-based) memory allocation as long as the
39 * value fits on 32 bits. This makes it substantially faster than usual
40 * "big integer" implementations (including .NET's implementation, since
41 * version 4.0) when values are small.
43 * Instances are immutable, and thus can be used as if they were
46 * None of this code is "constant-time". As such, ZInt should be
47 * considered unsuitable to implementations of cryptographic algorithms.
50 public struct ZInt : IComparable, IComparable<ZInt>, IEquatable<ZInt> {
55 * If varray == null, then "small" contains the integer value.
57 * If varray != null, then it contains the value, in little-endian
58 * convention (least significant word comes first) and of
59 * minimal encoded length (i.e. "trimmed"). Two's complement is
60 * used for negative values. "small" is then -1 or 0, depending
61 * on whether the value is negative or not.
63 * Note that the trimmed value does not always include the sign
66 * If the integer value is in the range of values which can be
67 * represented in an "int", then magn is null. There is no allowed
68 * overlap between the two kinds of encodings.
70 * Default value thus encodes the integer zero.
74 readonly uint[] varray;
79 public static ZInt MinusOne {
80 get { return new ZInt(-1, null); }
86 public static ZInt Zero {
87 get { return new ZInt(0, null); }
93 public static ZInt One {
94 get { return new ZInt(1, null); }
100 public static ZInt Two {
101 get { return new ZInt(2, null); }
105 * Internal constructor which assumes that the provided values are
106 * correct and normalized.
108 private ZInt(int small, uint[] varray)
111 this.varray = varray;
113 if (varray != null) {
114 if (small != -1 && small != 0) {
116 "Bad sign encoding: " + small);
118 if (varray.Length == 0) {
119 throw new Exception("Empty varray");
121 if (Length(small, varray) != varray.Length) {
122 throw new Exception("Untrimmed varray");
124 if (varray.Length == 1) {
126 * If there was room for a sign bit, then
127 * the "small" encoding should have been used.
129 if (((varray[0] ^ (uint)small) >> 31) == 0) {
131 "suboptimal encoding");
139 * Main internal build method. This method normalizes the encoding
140 * ("small" encoding is used when possible; otherwise, the varray
141 * is trimmed and the sign is normalized to -1 or 0).
143 * If "varray" is null, then the small value is used. Otherwise,
144 * only the sign bit (most significant bit) of "small" is used,
145 * and that value is normalized; the array is trimmed. If the
146 * small encoding then becomes applicable, then it is used.
148 private static ZInt Make(int small, uint[] varray)
150 if (varray == null) {
151 return new ZInt(small, null);
154 int n = Length(small, varray);
155 if (n == 1 && (((varray[0] ^ (uint)small) >> 31) == 0)) {
156 small = (int)varray[0];
160 * Note: if n == 0 then the value is -1 or 0, and
161 * "small" already contains the correct value;
162 * Trim() will then return null, which is appropriate.
164 varray = Trim(small, varray, n);
166 return new ZInt(small, varray);
170 * Create an instance from a signed 32-bit integer.
179 * Create an instance from an unsigned 32-bit integer.
181 public ZInt(uint val)
186 varray = new uint[1] { val };
193 * Create an instance from a signed 64-bit integer.
195 public ZInt(long val)
198 if ((long)small == val) {
201 ulong uval = (ulong)val;
202 uint w0 = (uint)uval;
203 uint w1 = (uint)(uval >> 32);
206 varray = new uint[1] { w0 };
207 } else if (w1 == 0xFFFFFFFF) {
209 varray = new uint[1] { w0 };
211 small = (int)w1 >> 31;
212 varray = new uint[2] { w0, w1 };
218 * Create an instance from an unsigned 64-bit integer.
220 public ZInt(ulong val)
222 if (val <= 0x7FFFFFFF) {
228 uint w1 = (uint)(val >> 32);
230 varray = new uint[1] { w0 };
232 varray = new uint[2] { w0, w1 };
238 * Create a ZInt instance for an integer expressed as an array
239 * of 32-bit integers (unsigned little-endian convention), with
240 * a specific number of value bits.
242 public static ZInt Make(uint[] words, int numBits)
247 int n = (numBits + 31) >> 5;
248 int kb = numBits & 31;
249 uint[] m = new uint[n];
250 Array.Copy(words, 0, m, 0, n);
252 m[n - 1] &= ~((uint)0xFFFFFFFF << kb);
258 * Create a ZInt instance for an integer expressed as an array
259 * of 64-bit integers (unsigned little-endian convention), with
260 * a specific number of value bits.
262 public static ZInt Make(ulong[] words, int numBits)
267 int kw = (numBits + 63) >> 6;
268 int kb = numBits & 63;
270 uint[] m = new uint[n];
272 ulong z = words[kw - 1]
273 & ~((ulong)0xFFFFFFFFFFFFFFFF << kb);
274 m[n - 1] = (uint)(z >> 32);
279 for (int i = kw - 1, j = n - 2; i >= 0; i --, j -= 2) {
282 m[j + 1] = (uint)(z >> 32);
288 * Test whether this value is 0.
292 return small == 0 && varray == null;
297 * Test whether this value is 1.
306 * Test whether this value is even.
310 uint w = (varray == null) ? (uint)small : varray[0];
316 * Test whether this value is a power of 2. Note that 1 is a power
317 * of 2 (but 0 is not). For negative values, false is returned.
319 public bool IsPowerOfTwo {
324 if (varray == null) {
325 return small != 0 && (small & -small) == small;
327 int n = varray.Length;
328 int z = (int)varray[n - 1];
332 for (int i = n - 2; i >= 0; i --) {
333 if (varray[i] != 0) {
342 * Get the sign of this value as an integer (-1 for negative
343 * values, 1 for positive, 0 for zero).
347 if (varray == null) {
350 } else if (small == 0) {
362 * Test whether this value would fit in an 'int'.
366 return varray == null;
371 * Test whether this value would fit in a "long".
375 if (varray == null) {
378 int n = varray.Length;
382 return ((int)varray[1] >> 31) == small;
390 * Get the length of the value in bits. This is the minimal number
391 * of bits of the two's complement representation of the value,
392 * excluding the sign bit; thus, both 0 and -1 have bit length 0.
394 public int BitLength {
396 if (varray == null) {
398 return 32 - LeadingZeros(~(uint)small);
400 return 32 - LeadingZeros((uint)small);
403 int n = varray.Length;
405 uint w = varray[n - 1];
407 return bl - LeadingZeros(~w);
409 return bl - LeadingZeros(w);
416 * Test whether the specified bit has value 1 or 0 ("true" is
417 * returned if the bit has value 1). Note that for negative values,
418 * two's complement representation is assumed.
420 public bool TestBit(int n)
423 throw new ArgumentOutOfRangeException();
425 if (varray == null) {
429 return (((uint)small >> n) & (uint)1) != 0;
433 if (nw >= varray.Length) {
437 return ((varray[nw] >> nb) & (uint)1) != 0;
442 * Copy some bits from this instance to the provided array. Bits
443 * are copied in little-endian order. First bit to be copied is
444 * bit at index "off", and exactly "num" bits are copied. This
445 * method modifies only the minimum number of destination words
446 * (i.e. the first "(num+31)/32" words, exactly). Remaining bits
447 * in the last touched word are set to 0.
449 public void CopyBits(int off, int num, uint[] dest)
451 CopyBits(off, num, dest, 0);
454 public void CopyBits(int off, int num, uint[] dest, int destOff)
456 if (off < 0 || num < 0) {
457 throw new ArgumentOutOfRangeException();
468 uint hmask = ~((uint)0xFFFFFFFF << kb);
469 if (x.varray == null) {
471 dest[destOff] = (uint)x.small & hmask;
473 uint iw = (uint)(x.small >> 31);
474 dest[destOff] = (uint)x.small;
475 for (int i = 1; i < kw; i ++) {
476 dest[destOff + i] = iw;
479 dest[destOff + kw] = iw & hmask;
483 int n = x.varray.Length;
485 Array.Copy(x.varray, 0, dest, destOff, kw);
487 Array.Copy(x.varray, 0, dest, destOff, n);
488 for (int i = n; i < kw; i ++) {
489 dest[destOff + i] = (uint)x.small;
495 last = x.varray[kw] & hmask;
497 last = (uint)x.small & hmask;
499 dest[destOff + kw] = last;
505 * Copy some bits from this instance to the provided array. Bits
506 * are copied in little-endian order. First bit to be copied is
507 * bit at index "off", and exactly "num" bits are copied. This
508 * method modifies only the minimum number of destination words
509 * (i.e. the first "(num+63)/64" words, exactly). Remaining bits
510 * in the last touched word are set to 0.
512 public void CopyBits(int off, int num, ulong[] dest)
514 CopyBits(off, num, dest, 0);
517 public void CopyBits(int off, int num, ulong[] dest, int destOff)
519 if (off < 0 || num < 0) {
520 throw new ArgumentOutOfRangeException();
531 ulong hmask = ~((ulong)0xFFFFFFFFFFFFFFFF << kb);
532 long xs = (long)x.small;
533 if (x.varray == null) {
535 dest[destOff] = (ulong)xs & hmask;
537 ulong iw = (ulong)(xs >> 31);
538 dest[destOff] = (ulong)xs;
539 for (int i = 1; i < kw; i ++) {
540 dest[destOff + i] = iw;
543 dest[destOff + kw] = iw & hmask;
547 int n = x.varray.Length;
548 uint iw = (uint)x.small;
550 for (int i = 0; i < kw; i ++, j += 2) {
551 uint w0 = (j < n) ? x.varray[j] : iw;
552 uint w1 = ((j + 1) < n) ? x.varray[j + 1] : iw;
554 (ulong)w0 | ((ulong)w1 << 32);
557 uint w0 = (j < n) ? x.varray[j] : iw;
558 uint w1 = ((j + 1) < n) ? x.varray[j + 1] : iw;
559 ulong last = (ulong)w0 | ((ulong)w1 << 32);
560 dest[destOff + kw] = last & hmask;
566 * Extract a 32-bit word at a given offset (counted in bits).
567 * This function is equivalent to right-shifting the value by
568 * "off" bits, then returning the low 32 bits (however, this
569 * function may be more efficient).
571 public uint GetWord(int off)
574 throw new ArgumentOutOfRangeException();
576 if (varray == null) {
581 return (uint)(x >> off);
583 int n = varray.Length;
598 uint lo = varray[kw];
599 return (lo >> kb) | (hi << (32 - kb));
604 * Extract a 64-bit word at a given offset (counted in bits).
605 * This function is equivalent to right-shifting the value by
606 * "off" bits, then returning the low 64 bits (however, this
607 * function may be more efficient).
609 public ulong GetWord64(int off)
612 throw new ArgumentOutOfRangeException();
614 if (varray == null) {
619 return (ulong)(x >> off);
621 int n = varray.Length;
629 return (ulong)varray[kw]
630 | ((ulong)small << 32);
632 return (ulong)varray[kw]
633 | ((ulong)varray[kw + 1] << 32);
641 } else if (kw == (n - 2)) {
650 uint lo = (v0 >> kb) | (v1 << (32 - kb));
651 uint hi = (v1 >> kb) | (v2 << (32 - kb));
652 return (ulong)lo | ((ulong)hi << 32);
657 * Convert this value to an 'int', using silent truncation if
658 * the value does not fit.
662 return (varray == null) ? small : (int)varray[0];
667 * Convert this value to an 'uint', using silent truncation if
668 * the value does not fit.
672 return (varray == null) ? (uint)small : varray[0];
677 * Convert this value to a 'long', using silent truncation if
678 * the value does not fit.
682 return (long)ToULong;
687 * Convert this value to an 'ulong', using silent truncation if
688 * the value does not fit.
690 public ulong ToULong {
692 if (varray == null) {
694 } else if (varray.Length == 1) {
695 uint iw = (uint)small;
696 return (ulong)varray[0] | ((ulong)iw << 32);
698 return (ulong)varray[0]
699 | ((ulong)varray[1] << 32);
705 * Get the actual length of a varray encoding: this is the minimal
706 * length, in words, needed to encode the value. The value sign
707 * is provided as a negative or non-negative integer, and the
708 * encoding of minimal length does not necessarily include a
709 * sign bit. The value 0 is returned when the array encodes 0
710 * or -1 (depending on sign).
712 static int Length(int sign, uint[] m)
717 uint sw = (uint)(sign >> 31);
719 while (n > 0 && m[n - 1] == sw) {
726 * Trim an encoding to its minimal encoded length. If the provided
727 * array is already of minimal length, it is returned unchanged.
729 static uint[] Trim(int sign, uint[] m)
731 int n = Length(sign, m);
734 } else if (n < m.Length) {
735 uint[] mt = new uint[n];
736 Array.Copy(m, 0, mt, 0, n);
744 * Trim or extend a value to the provided length. The returned
745 * array will have the length specified as "n" (if n == 0, then
746 * null is returned). If the source array already has the right
747 * size, then it is returned unchanged.
749 static uint[] Trim(int sign, uint[] m, int n)
753 } else if (m == null) {
762 uint[] r = new uint[n];
763 Array.Copy(m, 0, r, 0, ct);
765 } else if (ct == n) {
768 uint[] r = new uint[n];
769 Array.Copy(m, 0, r, 0, n);
771 Fill(0xFFFFFFFF, r, ct, n - ct);
777 static void Fill(uint val, uint[] buf)
779 Fill(val, buf, 0, buf.Length);
782 static void Fill(uint val, uint[] buf, int off, int len)
789 // =================================================================
793 * The methods whose name begins with "Mutate" modify the array
794 * they are given as first parameter; the other methods instantiate
797 * As a rule, untrimmed arrays are accepted as input, and output
798 * may be untrimmed as well.
802 * Count the number of leading zeros for a 32-bit value (number of
803 * consecutive zeros, starting with the most significant bit). This
804 * value is between 0 (for a value equal to 2^31 or greater) and
807 static int LeadingZeros(uint v)
813 if (v > 0xFFFF) { v >>= 16; } else { n += 16; }
814 if (v > 0x00FF) { v >>= 8; } else { n += 8; }
815 if (v > 0x000F) { v >>= 4; } else { n += 4; }
816 if (v > 0x0003) { v >>= 2; } else { n += 2; }
817 if (v <= 0x0001) { n ++; }
822 * Duplicate the provided magnitude array. No attempt is made at
823 * trimming. The source array MUST NOT be null.
825 static uint[] Dup(uint[] m)
827 uint[] r = new uint[m.Length];
828 Array.Copy(m, 0, r, 0, m.Length);
833 * Increment the provided array. If there is a resulting carry,
834 * then "true" is returned, "false" otherwise. The array MUST
837 static bool MutateIncr(uint[] x)
840 for (int i = 0; i < n; i ++) {
851 * Decrement the provided array. If there is a resulting carry,
852 * then "true" is returned, "false" otherwise. The array MUST
855 static bool MutateDecr(uint[] x)
858 for (int i = 0; i < n; i ++) {
869 * Multiply a[] with b[] (unsigned multiplication).
871 static uint[] Mul(uint[] a, uint[] b)
873 // TODO: use Karatsuba when operands are large.
874 int na = Length(0, a);
875 int nb = Length(0, b);
876 if (na == 0 || nb == 0) {
879 uint[] r = new uint[na + nb];
880 for (int i = 0; i < na; i ++) {
883 for (int j = 0; j < nb; j ++) {
884 ulong mb = (ulong)b[j];
885 ulong mr = (ulong)r[i + j];
886 ulong w = ma * mb + mr + carry;
890 r[i + nb] = (uint)carry;
896 * Get the sign and magnitude of an integer. The sign is
897 * normalized to -1 (negative) or 0 (positive or zero). The
898 * magnitude is an array of length at least 1, containing the
899 * absolute value of this integer; if possible, the varray
900 * is reused (hence, the magnitude array MUST NOT be altered).
902 static void ToAbs(ZInt x, out int sign, out uint[] magn)
912 magn = new uint[1] { (uint)x.small };
917 * Compare two integers, yielding -1, 0 or 1.
919 static int Compare(int a, int b)
931 * Compare a[] with b[] (unsigned). Returned value is 1 if a[]
932 * is greater than b[], 0 if they are equal, -1 otherwise.
934 static int Compare(uint[] a, uint[] b)
936 int ka = Length(0, a);
937 int kb = Length(0, b);
940 } else if (ka == kb) {
947 } else if (wa > wb) {
958 * Add b[] to a[] (unsigned). a[] is modified "in place". Only
959 * n words of a[] are modified. Moreover, the value of
960 * b[] which is added is left-shifted: words b[0]...b[n-1] are
961 * added to a[k]...a[k+n-1]. The final carry is returned ("true"
962 * for 1, "false" for 0). Neither a nor b may be null.
964 static bool MutateAdd(uint[] a, int n, uint[] b, int k)
967 for (int i = 0; i < n; i ++) {
983 * Substract b[] from a[] (unsigned). a[] is modified "in
984 * place". Only n words of a[] are modified. Words
985 * b[0]...b[n-1] are subtracted from a[k]...a[k+n-1]. The final
986 * carry is returned ("true" for -1, "false" for 0). Neither a
989 static bool MutateSub(uint[] a, int n, uint[] b, int k)
992 for (int i = 0; i < n; i ++) {
1008 * Get the length (in words) of the result of a left-shift of
1009 * the provided integer by k bits. If k < 0, then the value is
1010 * computed for a right-shift by -k bits.
1012 static int GetLengthForLeftShift(ZInt x, int k)
1015 if (k == Int32.MinValue) {
1018 return GetLengthForRightShift(x, -k);
1020 uint bl = (uint)x.BitLength + (uint)k;
1021 return (int)((bl + 31) >> 5);
1025 * Get the length (in words) of the result of a right-shift of
1026 * the provided integer by k bits. If k < 0, then the value is
1027 * computed for a left-shift by -k bits.
1029 static int GetLengthForRightShift(ZInt x, int k)
1032 if (k == Int32.MinValue) {
1033 throw new OverflowException();
1035 return GetLengthForLeftShift(x, -k);
1037 uint bl = (uint)x.BitLength;
1038 if (bl <= (uint)k) {
1041 return (int)((bl - k + 31) >> 5);
1046 * Left-shift a[] (unsigned) by k bits. If k < 0, then this becomes
1049 static uint[] ShiftLeft(uint[] a, int k)
1052 return ShiftRight(a, -k);
1053 } else if (k == 0) {
1056 int n = Length(0, a);
1062 * Allocate the result array, with the exact proper size.
1064 int bl = ((n << 5) - LeadingZeros(a[n - 1])) + k;
1065 uint[] r = new uint[(bl + 31) >> 5];
1071 * Special case: shift by an integral amount of words.
1074 Array.Copy(a, 0, r, kw, n);
1079 * Copy the bits. This loop handles one source word at
1080 * a time, and writes one destination word at a time.
1081 * Some unhandled bits may remain at the end.
1085 for (int i = 0; i < n; i ++) {
1087 r[i + kw] = bits | (w << kb);
1097 * Right-shift a[] by k bits. If k < 0, then this becomes
1100 static uint[] ShiftRight(uint[] a, int k)
1103 return ShiftLeft(a, -k);
1104 } else if (k == 0) {
1107 int n = Length(0, a);
1111 int bl = (n << 5) - LeadingZeros(a[n - 1]) - k;
1115 uint[] r = new uint[(bl + 31) >> 5];
1121 * Special case: shift by an integral amount of words.
1124 Array.Copy(a, kw, r, 0, r.Length);
1129 * Copy the bits. This loop handles one source word at
1130 * a time, and writes one destination word at a time.
1131 * Some unhandled bits may remain at the end.
1133 uint bits = a[kw] >> kb;
1135 for (int i = kw + 1; i < n; i ++) {
1137 r[i - kw - 1] = bits | (w << zkb);
1141 r[n - kw - 1] = bits;
1147 * Euclidian division of a[] (unsigned) by b (single word). This
1148 * method assumes that b is not 0.
1150 static void DivRem(uint[] a, uint b, out uint[] q, out uint r)
1152 int n = Length(0, a);
1160 for (int i = n - 1; i >= 0; i --) {
1162 * Performance: we strongly hope that the JIT
1163 * compiler will notice that the "/" and "%"
1164 * can be combined into a single operation.
1165 * TODO: test whether we should replace the
1166 * carry computation with:
1167 * carry = w - (ulong)q[i] * b;
1169 ulong w = (ulong)a[i] + (carry << 32);
1170 q[i] = (uint)(w / b);
1177 * Euclidian division of a[] (unsigned) by b (single word). This
1178 * method assumes that b is not 0. a[] is modified in place. The
1179 * remainder (in the 0..b-1 range) is returned.
1181 static uint MutateDivRem(uint[] a, uint b)
1183 int n = Length(0, a);
1188 for (int i = n - 1; i >= 0; i --) {
1190 * Performance: we strongly hope that the JIT
1191 * compiler will notice that the "/" and "%"
1192 * can be combined into a single operation.
1193 * TODO: test whether we should replace the
1194 * carry computation with:
1195 * carry = w - (ulong)q[i] * b;
1197 ulong w = (ulong)a[i] + (carry << 32);
1198 a[i] = (uint)(w / b);
1205 * Euclidian division of a[] by b[] (unsigned). This method
1206 * assumes that b[] is neither 0 or 1, and a[] is not smaller
1207 * than b[] (this implies that the quotient won't be zero).
1209 static void DivRem(uint[] a, uint[] b, out uint[] q, out uint[] r)
1211 int nb = Length(0, b);
1214 * Special case when the divisor fits on one word.
1218 DivRem(a, b[0], out q, out r[0]);
1225 * We first normalize divisor and dividend such that the
1226 * most significant bit of the most significant word of
1227 * the divisor is set. We can then compute the quotient
1228 * word by word. In details:
1231 * w = 2^32 (one word)
1232 * a = (w*a0 + a1) * w^N + a2
1241 * In other words, a0 and a1 are the two upper words of a[],
1242 * b0 is the upper word of b[] and has length 32 bits exactly,
1243 * and the quotient of a by b fits in one word.
1245 * Under these conditions, define q and r such that:
1249 * We can then compute a value u this way:
1250 * if a0 = b0, then let u = w-1
1251 * otherwise, let u be such that
1252 * (w*a0 + a1) = u*b0 + v, where 0 <= v < b0
1253 * It can then be shown that all these inequations hold:
1257 * In plain words, this allows us to compute an almost-exact
1258 * estimate of the upper word of the quotient, with only
1259 * one 64-bit division.
1263 * Normalize dividend and divisor. The normalized dividend
1264 * will become the temporary array for the remainder, and
1265 * we will modify it, so we make sure we have a copy.
1267 int norm = LeadingZeros(b[nb - 1]);
1268 r = ShiftLeft(a, norm);
1270 r = new uint[a.Length];
1271 Array.Copy(a, 0, r, 0, a.Length);
1273 uint[] b2 = ShiftLeft(b, norm);
1274 int nr = Length(0, r);
1276 if (Length(0, b2) != nb) {
1277 throw new Exception("normalize error 1");
1279 if (b2[nb - 1] < 0x80000000) {
1280 throw new Exception("normalize error 2");
1283 uint[] ta = ShiftRight(r, norm);
1284 if (Compare(a, ta) != 0) {
1285 throw new Exception("normalize error 3");
1287 uint[] tb = ShiftRight(b2, norm);
1288 if (Compare(b, tb) != 0) {
1289 throw new Exception("normalize error 4");
1296 * Length of the quotient will be (at most) k words. This
1297 * is the number of iterations in the loop below.
1299 int k = (nr - nb) + 1;
1302 throw new Exception("wrong iteration count: " + k);
1308 * The upper word of a[] (the one we modify, i.e. currently
1309 * stored in r[]) is in a0; it is carried over from the
1310 * previous loop iteration. Initially, it is zero.
1313 uint b0 = b[nb - 1];
1321 ulong ah = ((ulong)a0 << 32) | (ulong)a1;
1322 u = (uint)(ah / b0);
1326 * Candidate word for the quotient:
1327 * -- if u = 0 then qw is necessarily 0, and the
1328 * rest of this iteration is trivial;
1329 * -- if u = 1 then we try qw = 1;
1330 * -- otherwise, we try qw = u-1.
1337 } else if (u == 1) {
1344 * "Trying" a candidate word means subtracting from
1345 * r[] the product qw*b (b[] being shifted by k words).
1346 * The result may be negative, in which case we
1347 * overestimated qw; it may be greater than b or
1348 * equal to b, in which case we underestimated qw;
1349 * or it may be just fine.
1353 for (int i = 0; i < nb; i ++) {
1355 ulong z = (ulong)wb * (ulong)qw + carry;
1358 uint wc = wa - (uint)z;
1369 * Once we have adjusted everything, the upper word
1370 * of r[] will be nullified; wo do it now. Note that
1371 * for the first loop iteration, that upper word
1372 * may be absent (so already zero, but also "virtual").
1379 * At that point, "carry" should be equal to a0
1380 * if we estimated right.
1382 if (carry < (ulong)a0) {
1388 if (carry + 1 != (ulong)a0) {
1389 throw new Exception("div error 1");
1391 if (!MutateSub(r, nb, b, k)) {
1392 throw new Exception("div error 2");
1395 MutateSub(r, nb, b, k);
1397 } else if (carry > (ulong)a0) {
1403 if (carry - 1 != (ulong)a0) {
1404 throw new Exception("div error 3");
1406 if (!MutateAdd(r, nb, b, k)) {
1407 throw new Exception("div error 4");
1410 MutateAdd(r, nb, b, k);
1412 } else if (tooBig) {
1414 * Underestimate, but no expected carry.
1418 if (MutateSub(r, nb, b, k)) {
1419 throw new Exception("div error 5");
1422 MutateSub(r, nb, b, k);
1431 * At that point, r[] contains the remainder but needs
1432 * to be shifted back, to account for the normalization
1433 * performed before. q[] is correct (but possibly
1436 r = ShiftRight(r, norm);
1439 // =================================================================
1442 * Conversion to sbyte; an OverflowException is thrown if the
1443 * value does not fit. Use ToInt to get a truncating conversion.
1445 public static explicit operator sbyte(ZInt val)
1448 if (x < SByte.MinValue || x > SByte.MaxValue) {
1449 throw new OverflowException();
1455 * Conversion to byte; an OverflowException is thrown if the
1456 * value does not fit. Use ToInt to get a truncating conversion.
1458 public static explicit operator byte(ZInt val)
1461 if (x > Byte.MaxValue) {
1462 throw new OverflowException();
1468 * Conversion to short; an OverflowException is thrown if the
1469 * value does not fit. Use ToInt to get a truncating conversion.
1471 public static explicit operator short(ZInt val)
1474 if (x < Int16.MinValue || x > Int16.MaxValue) {
1475 throw new OverflowException();
1481 * Conversion to ushort; an OverflowException is thrown if the
1482 * value does not fit. Use ToInt to get a truncating conversion.
1484 public static explicit operator ushort(ZInt val)
1487 if (x > UInt16.MaxValue) {
1488 throw new OverflowException();
1494 * Conversion to int; an OverflowException is thrown if the
1495 * value does not fit. Use ToInt to get a truncating conversion.
1497 public static explicit operator int(ZInt val)
1499 if (val.varray != null) {
1500 throw new OverflowException();
1506 * Conversion to uint; an OverflowException is thrown if the
1507 * value does not fit. Use ToUInt to get a truncating conversion.
1509 public static explicit operator uint(ZInt val)
1513 throw new OverflowException();
1515 uint[] m = val.varray;
1518 } else if (m.Length > 1) {
1519 throw new OverflowException();
1526 * Conversion to long; an OverflowException is thrown if the
1527 * value does not fit. Use ToLong to get a truncating conversion.
1529 public static explicit operator long(ZInt val)
1532 uint[] m = val.varray;
1535 } else if (m.Length == 1) {
1536 return (long)m[0] | ((long)s << 32);
1537 } else if (m.Length == 2) {
1540 if (((w1 ^ (uint)s) >> 31) != 0) {
1541 throw new OverflowException();
1543 return (long)w0 | ((long)w1 << 32);
1545 throw new OverflowException();
1550 * Conversion to ulong; an OverflowException is thrown if the
1551 * value does not fit. Use ToULong to get a truncating conversion.
1553 public static explicit operator ulong(ZInt val)
1557 throw new OverflowException();
1559 uint[] m = val.varray;
1562 } else if (m.Length == 1) {
1564 } else if (m.Length == 2) {
1565 return (ulong)m[0] | ((ulong)m[1] << 32);
1567 throw new OverflowException();
1572 * By definition, conversion from sbyte conserves the value.
1574 public static implicit operator ZInt(sbyte val)
1576 return new ZInt((int)val);
1580 * By definition, conversion from byte conserves the value.
1582 public static implicit operator ZInt(byte val)
1584 return new ZInt((uint)val);
1588 * By definition, conversion from short conserves the value.
1590 public static implicit operator ZInt(short val)
1592 return new ZInt((int)val);
1596 * By definition, conversion from ushort conserves the value.
1598 public static implicit operator ZInt(ushort val)
1600 return new ZInt((uint)val);
1604 * By definition, conversion from int conserves the value.
1606 public static implicit operator ZInt(int val)
1608 return new ZInt(val);
1612 * By definition, conversion from uint conserves the value.
1614 public static implicit operator ZInt(uint val)
1616 return new ZInt(val);
1620 * By definition, conversion from long conserves the value.
1622 public static implicit operator ZInt(long val)
1624 return new ZInt(val);
1628 * By definition, conversion from ulong conserves the value.
1630 public static implicit operator ZInt(ulong val)
1632 return new ZInt(val);
1636 * Unary '+' operator is a no-operation.
1638 public static ZInt operator +(ZInt a)
1646 public static ZInt operator -(ZInt a)
1649 uint[] m = a.varray;
1651 if (s == Int32.MinValue) {
1652 return new ZInt(0, new uint[1] { 0x80000000 });
1654 return new ZInt(-s, null);
1659 * Two's complement: invert all bits, then add 1. The
1660 * result array will usually have the same size, but may
1661 * be one word longer or one word shorter.
1663 int n = Length(s, m);
1664 uint[] bm = new uint[n];
1665 for (int i = 0; i < n; i ++) {
1668 if (MutateIncr(bm)) {
1670 * Extra carry. This may happen only if the source
1671 * array contained only zeros, which may happen only
1672 * for a source value -2^(32*k) for some integer k
1675 bm = new uint[n + 1];
1677 return new ZInt(0, bm);
1680 * The resulting array might be too big by one word,
1681 * so we must not assume that it is trimmed.
1683 return Make((int)~(uint)s, bm);
1688 * Addition operator.
1690 public static ZInt operator +(ZInt a, ZInt b)
1694 uint[] ma = a.varray;
1695 uint[] mb = b.varray;
1704 return new ZInt((long)sa + (long)sb);
1706 ma = new uint[1] { (uint)sa };
1708 } else if (mb == null) {
1712 mb = new uint[1] { (uint)sb };
1717 int n = Math.Max(na, nb) + 1;
1718 uint[] mc = new uint[n];
1720 for (int i = 0; i < n; i ++) {
1721 uint wa = i < na ? ma[i] : (uint)sa;
1722 uint wb = i < nb ? mb[i] : (uint)sb;
1723 ulong z = (ulong)wa + (ulong)wb + carry;
1727 return Make((-(int)carry) ^ sa ^ sb, mc);
1731 * Subtraction operator.
1733 public static ZInt operator -(ZInt a, ZInt b)
1737 uint[] ma = a.varray;
1738 uint[] mb = b.varray;
1747 return new ZInt((long)sa - (long)sb);
1749 ma = new uint[1] { (uint)sa };
1751 } else if (mb == null) {
1755 mb = new uint[1] { (uint)sb };
1760 int n = Math.Max(na, nb) + 1;
1761 uint[] mc = new uint[n];
1763 for (int i = 0; i < n; i ++) {
1764 uint wa = i < na ? ma[i] : (uint)sa;
1765 uint wb = i < nb ? mb[i] : (uint)sb;
1766 long z = (long)wa - (long)wb + carry;
1770 return Make((int)carry ^ sa ^ sb, mc);
1774 * Increment operator.
1776 public static ZInt operator ++(ZInt a)
1779 uint[] ma = a.varray;
1781 return new ZInt((long)s + 1);
1783 uint[] mb = Dup(ma);
1784 if (MutateIncr(mb)) {
1786 mb = new uint[n + 1];
1788 return new ZInt(0, mb);
1795 * Decrement operator.
1797 public static ZInt operator --(ZInt a)
1800 uint[] ma = a.varray;
1802 return new ZInt((long)s - 1);
1806 * MutateDecr() will report a carry only if the varray
1807 * contained only zeros; since this value was not small,
1808 * then it must have been negative.
1810 uint[] mb = Dup(ma);
1811 if (MutateDecr(mb)) {
1813 uint[] mc = new uint[n + 1];
1814 Array.Copy(mb, 0, mc, 0, n);
1816 return new ZInt(-1, mc);
1823 * Multiplication operator.
1825 public static ZInt operator *(ZInt a, ZInt b)
1829 uint[] ma = a.varray;
1830 uint[] mb = b.varray;
1834 * -- one of the operands is zero
1835 * -- both operands are small
1845 return new ZInt((long)sa * (long)sb);
1847 } else if (mb == null) {
1854 * Get both values in sign+magnitude representation.
1856 ToAbs(a, out sa, out ma);
1857 ToAbs(b, out sb, out mb);
1860 * Compute the product. Set the sign.
1862 ZInt r = Make(0, Mul(ma, mb));
1863 if ((sa ^ sb) < 0) {
1870 * Integer division: the quotient is returned, and the remainder
1871 * is written in 'r'.
1873 * This operation follows the C# rules for integer division:
1874 * -- rounding is towards 0
1875 * -- quotient is positive if dividend and divisor have the same
1876 * sign, negative otherwise
1877 * -- remainder has the sign of the dividend
1879 * Attempt at dividing by zero triggers a DivideByZeroException.
1881 public static ZInt DivRem(ZInt a, ZInt b, out ZInt r)
1884 DivRem(a, b, out q, out r);
1888 static void DivRem(ZInt a, ZInt b,
1889 out ZInt q, out ZInt r)
1893 uint[] ma = a.varray;
1894 uint[] mb = b.varray;
1897 * Division by zero triggers an exception.
1899 if (mb == null && sb == 0) {
1900 throw new DivideByZeroException();
1904 * If the dividend is zero, then both quotient and
1905 * remainder are zero.
1907 if (ma == null && sa == 0) {
1914 * If both dividend and divisor are small, then we
1915 * use the native 64-bit integer types. If only the
1916 * divisor is small, then we have a special fast case
1917 * for division by 1 or -1.
1919 if (ma == null && mb == null) {
1920 q = new ZInt((long)sa / (long)sb);
1921 r = new ZInt((long)sa % (long)sb);
1929 } else if (sb == -1) {
1937 * We know that the dividend is not 0, and the divisor
1938 * is not -1, 0 or 1. We now want the sign+magnitude
1939 * representations of both operands.
1941 ToAbs(a, out sa, out ma);
1942 ToAbs(b, out sb, out mb);
1945 * If the divisor is greater (in magnitude) than the
1946 * dividend, then the quotient is zero and the remainder
1947 * is equal to the dividend. If the divisor and dividend
1948 * are equal in magnitude, then the remainder is zero and
1949 * the quotient is 1 if divisor and dividend have the same
1950 * sign, -1 otherwise.
1952 int cc = Compare(ma, mb);
1957 } else if (cc == 0) {
1958 q = (sa == sb) ? One : MinusOne;
1964 * At that point, we know that the divisor is not -1, 0
1965 * or 1, and that the quotient will not be 0. We perform
1966 * the unsigned division (with the magnitudes), then
1967 * we adjust the signs.
1971 DivRem(ma, mb, out mq, out mr);
1974 * Quotient is positive if divisor and dividend have the
1975 * same sign, negative otherwise. Remainder always has
1976 * the sign of the dividend, but it may be zero.
1987 if (q * b + r != a) {
1988 throw new Exception("division error");
1994 * Division operator: see DivRem() for details.
1996 public static ZInt operator /(ZInt a, ZInt b)
1999 DivRem(a, b, out q, out r);
2004 * Remainder operator: see DivRem() for details.
2006 public static ZInt operator %(ZInt a, ZInt b)
2009 DivRem(a, b, out q, out r);
2014 * Reduce this value modulo the provided m. This differs from
2015 * '%' in that the returned value is always in the 0 to abs(m)-1
2018 public ZInt Mod(ZInt m)
2031 * Left-shift operator.
2033 public static ZInt operator <<(ZInt a, int n)
2036 if (n == Int32.MinValue) {
2042 uint[] ma = a.varray;
2045 return new ZInt((long)sa << n);
2051 uint[] mr = new uint[GetLengthForLeftShift(a, n)];
2054 * Special case when the shift is a multiple of 32.
2061 return new ZInt(sa >> 31, mr);
2063 Array.Copy(ma, 0, mr, kw, ma.Length);
2064 return new ZInt(sa, mr);
2069 * At that point, we know that the source integer does
2070 * not fit in a signed "int", or is shifted by more than
2071 * 32 bits, or both. Either way, the result will not fit
2074 * We process all input words one by one.
2086 for (int i = 0; i < j; i ++) {
2088 mr[i + kw] = rem | (wa << kb);
2094 if ((j + kw) == mr.Length - 1) {
2095 if (rem == ((uint)sa >> ikb)) {
2096 throw new Exception(
2097 "Wrong left-shift: untrimmed");
2099 } else if ((j + kw) == mr.Length) {
2100 if (rem != ((uint)sa >> ikb)) {
2101 throw new Exception(
2102 "Wrong left-shift: dropped bits");
2105 throw new Exception(
2106 "Wrong left-shift: oversized output length");
2109 if ((j + kw) < mr.Length) {
2110 mr[j + kw] = rem | ((uint)sa << kb);
2112 return new ZInt(sa, mr);
2116 * Right-shift operator.
2118 public static ZInt operator >>(ZInt a, int n)
2121 if (n == Int32.MinValue) {
2122 throw new OverflowException();
2127 uint[] ma = a.varray;
2130 * If right-shifting a "small" value, then we can
2131 * do the computation with the native ">>" operator
2132 * on "int" values, unless the shift count is 32
2133 * or more, in which case we get either 0 or -1,
2134 * depending on the source value sign.
2137 return new ZInt(sa >> n, null);
2139 return new ZInt(sa >> 31, null);
2144 * At that point, we know that the source value uses
2145 * a non-null varray. We compute the bit length of the
2146 * result. If the result would fit in an "int" (bit length
2147 * of 31 or less) then we handle it as a special case.
2151 int bl = a.BitLength - n;
2153 return new ZInt(sa, null);
2154 } else if (bl <= 31) {
2156 return new ZInt((int)ma[kw], null);
2160 uint w1 = (kw + 1) < p ? ma[kw + 1] : (uint)sa;
2161 return new ZInt((int)((w0 >> kb)
2162 | (w1 << (32 - kb))), null);
2167 * Result will require an array. Let's allocate it.
2169 uint[] mr = new uint[(bl + 31) >> 5];
2172 * Special case when the shift is a multiple of 32.
2176 if (mr.Length != (ma.Length - kw)) {
2177 throw new Exception(
2178 "Wrong right-shift: output length");
2181 Array.Copy(ma, kw, mr, 0, ma.Length - kw);
2182 return new ZInt(sa, mr);
2186 * We process all input words one by one.
2189 uint rem = ma[kw] >> kb;
2191 for (int i = kw + 1; i < j; i ++) {
2193 mr[i - kw - 1] = rem | (wa << ikb);
2197 if ((j - kw - 1) == mr.Length - 1) {
2198 if (rem == ((uint)sa >> kb)) {
2199 throw new Exception(
2200 "Wrong right-shift: untrimmed");
2202 } else if ((j - kw - 1) == mr.Length) {
2203 if (rem != ((uint)sa >> kb)) {
2204 throw new Exception(
2205 "Wrong right-shift: dropped bits");
2208 throw new Exception(
2209 "Wrong right-shift: oversized output length");
2212 if ((j - kw - 1) < mr.Length) {
2213 mr[j - kw - 1] = rem | ((uint)sa << ikb);
2215 return new ZInt(sa, mr);
2219 * NOTES ON BITWISE BOOLEAN OPERATIONS
2221 * When both operands are "small" then the result is necessarily
2222 * small: in "small" encoding, all bits beyond bit 31 of the
2223 * two's complement encoding are equal to bit 31, so the result
2224 * of computing the operation on bits 31 will also be valid for
2225 * all subsequent bits. Therefore, when the two operands are
2226 * small, we can just do the operation on the "int" values and the
2227 * result is guaranteed "small" as well.
2231 * Compute the bitwise AND between a "big" and a "small" values.
2233 static ZInt AndSmall(int s, uint[] m, int x)
2238 return new ZInt(s, r);
2240 return new ZInt(x & (int)m[0], null);
2245 * Bitwise AND operator.
2247 public static ZInt operator &(ZInt a, ZInt b)
2251 uint[] ma = a.varray;
2252 uint[] mb = b.varray;
2255 return new ZInt(sa & sb, null);
2257 return AndSmall(sb, mb, sa);
2259 } else if (mb == null) {
2260 return AndSmall(sa, ma, sb);
2264 * Both values are big. Since upper zero bits force the
2265 * result to contain zeros, the result size is that of the
2266 * positive operand (the smallest of the two if both are
2267 * positive). If both operands are negative, then the
2268 * result magnitude may be as large as the largest of the
2269 * two source magnitudes.
2271 * Result is negative if and only if both operands are
2279 nr = Math.Min(na, nb);
2287 nr = Math.Max(na, nb);
2290 uint[] mr = new uint[nr];
2291 for (int i = 0; i < nr; i ++) {
2292 uint wa = i < na ? ma[i] : (uint)sa;
2293 uint wb = i < nb ? mb[i] : (uint)sb;
2296 return Make(sa & sb, mr);
2300 * Compute the bitwise OR between a "big" value (sign and
2301 * magnitude, already normalized/trimmed), and a small value
2302 * (which can be positive or negative).
2304 static ZInt OrSmall(int s, uint[] m, int x)
2307 return new ZInt(x | (int)m[0], null);
2311 return new ZInt(s, r);
2316 * Bitwise OR operator.
2318 public static ZInt operator |(ZInt a, ZInt b)
2322 uint[] ma = a.varray;
2323 uint[] mb = b.varray;
2326 return new ZInt(sa | sb, null);
2328 return OrSmall(sb, mb, sa);
2330 } else if (mb == null) {
2331 return OrSmall(sa, ma, sb);
2335 * Both values are big. Since upper one bits force the
2336 * result to contain ones, the result size is that of
2337 * the negative operand (the greater, i.e. "smallest" of
2338 * the two if both are negative). If both operands are
2339 * positive, then the result magnitude may be as large
2340 * as the largest of the two source magnitudes.
2342 * Result is positive if and only if both operands are
2350 nr = Math.Max(na, nb);
2358 nr = Math.Min(na, nb);
2361 uint[] mr = new uint[nr];
2362 for (int i = 0; i < nr; i ++) {
2363 uint wa = i < na ? ma[i] : (uint)sa;
2364 uint wb = i < nb ? mb[i] : (uint)sb;
2367 return Make(sa | sb, mr);
2371 * Bitwise XOR operator.
2373 public static ZInt operator ^(ZInt a, ZInt b)
2377 uint[] ma = a.varray;
2378 uint[] mb = b.varray;
2379 if (ma == null && mb == null) {
2380 return new ZInt(sa ^ sb, null);
2392 uint[] mx = new uint[nx];
2393 mx[0] = ma[0] ^ (uint)sb;
2396 for (int i = 1; i < nx; i ++) {
2400 Array.Copy(ma, 1, mx, 1, nx - 1);
2403 return Make(sa ^ (sb >> 31), mx);
2407 * Both operands use varrays.
2408 * Result can be as big as the bigger of the two operands
2409 * (it _will_ be that big, necessarily, if the two operands
2410 * have distinct sizes).
2414 int nr = Math.Max(na, nb);
2415 uint[] mr = new uint[nr];
2416 for (int i = 0; i < nr; i ++) {
2417 uint wa = i < na ? ma[i] : (uint)sa;
2418 uint wb = i < nb ? mb[i] : (uint)sb;
2421 return Make((sa ^ sb) >> 31, mr);
2425 * Bitwise inversion operator.
2427 public static ZInt operator ~(ZInt a)
2430 uint[] m = a.varray;
2432 return new ZInt(~s, null);
2435 uint[] r = new uint[n];
2436 for (int i = 0; i < n; i ++) {
2439 return new ZInt(~s, r);
2444 * Basic comparison; returned value is -1, 0 or 1 depending on
2445 * whether this instance is to be considered lower then, equal
2446 * to, or greater than the provided object.
2448 * All ZInt instances are considered greater than 'null', and
2449 * lower than any non-null object that is not a ZInt.
2451 public int CompareTo(object obj)
2456 if (!(obj is ZInt)) {
2459 return CompareTo((ZInt)obj);
2463 * Basic comparison; returned value is -1, 0 or 1 depending on
2464 * whether this instance is to be considered lower then, equal
2465 * to, or greater than the provided value.
2467 public int CompareTo(ZInt v)
2470 uint[] mv = v.varray;
2471 int sign1 = small >> 31;
2472 int sign2 = sv >> 31;
2473 if (sign1 != sign2) {
2475 * One of the sign* values is -1, the other is 0.
2477 return sign1 - sign2;
2481 * Both values have the same sign. Since the varrays are
2482 * trimmed, we can use their presence and length to
2483 * quickly resolve most cases.
2486 if (varray == null) {
2488 return Compare(small, sv);
2496 int n1 = varray.Length;
2500 } else if (n1 == n2) {
2501 return Compare(varray, mv);
2508 if (varray == null) {
2510 return Compare(small, sv);
2518 return Compare(varray, mv);
2525 * Equality comparison: a ZInt instance is equal only to another
2526 * ZInt instance that encodes the same integer value.
2528 public override bool Equals(object obj)
2533 if (!(obj is ZInt)) {
2536 return CompareTo((ZInt)obj) == 0;
2540 * Equality comparison: a ZInt instance is equal only to another
2541 * ZInt instance that encodes the same integer value.
2543 public bool Equals(ZInt v)
2545 return CompareTo(v) == 0;
2549 * The hash code for a ZInt is equal to its lower 32 bits.
2551 public override int GetHashCode()
2553 if (varray == null) {
2556 return (int)varray[0];
2561 * Equality operator.
2563 public static bool operator ==(ZInt a, ZInt b)
2565 return a.CompareTo(b) == 0;
2569 * Inequality operator.
2571 public static bool operator !=(ZInt a, ZInt b)
2573 return a.CompareTo(b) != 0;
2577 * Lower-than operator.
2579 public static bool operator <(ZInt a, ZInt b)
2581 return a.CompareTo(b) < 0;
2585 * Lower-or-equal operator.
2587 public static bool operator <=(ZInt a, ZInt b)
2589 return a.CompareTo(b) <= 0;
2593 * Greater-than operator.
2595 public static bool operator >(ZInt a, ZInt b)
2597 return a.CompareTo(b) > 0;
2601 * Greater-or-equal operator.
2603 public static bool operator >=(ZInt a, ZInt b)
2605 return a.CompareTo(b) >= 0;
2609 * Power function: this raises x to the power e. The exponent e
2610 * MUST NOT be negative. If x and e are both zero, then 1 is
2613 public static ZInt Pow(ZInt x, int e)
2616 throw new ArgumentOutOfRangeException();
2621 if (e == 1 || x.IsZero || x.IsOne) {
2629 if (x.IsPowerOfTwo) {
2630 int t = x.BitLength - 1;
2631 long u = (long)t * (long)e;
2632 if (u > (long)Int32.MaxValue) {
2633 throw new OverflowException();
2650 return neg ? -x : x;
2654 * Modular exponentation: this function raises v to the power e
2655 * modulo m. The returned value is reduced modulo m: it will be
2656 * in the 0 to abs(m)-1 range.
2658 * The modulus m must be positive. If m is 1, then the result is
2659 * 0 (regardless of the values of v and e).
2661 * The exponent e must be nonnegative (this function does not
2662 * compute modular inverses). If e is zero, then the result is 1
2663 * (except if m is 1).
2665 public static ZInt ModPow(ZInt v, ZInt e, ZInt m)
2669 throw new ArgumentOutOfRangeException();
2674 } else if (sm == 0) {
2675 throw new DivideByZeroException();
2677 if (m.varray == null && m.small == 1) {
2687 // TODO: use Montgomery's multiplication when the exponent
2690 for (int n = e.BitLength - 2; n >= 0; n --) {
2700 * Get the absolute value of a ZInt.
2702 public static ZInt Abs(ZInt x)
2704 return (x.Sign < 0) ? -x : x;
2707 private static void AppendHex(StringBuilder sb, uint w, bool trim)
2711 for (; i >= 0 && (w >> i) == 0; i -= 4);
2717 for (; i >= 0; i -= 4) {
2718 sb.Append("0123456789ABCDEF"[(int)((w >> i) & 0x0F)]);
2723 * Convert this value to hexadecimal. If this instance is zero,
2724 * then "0" is returned. Otherwise, the number of digits is
2725 * minimal (no leading '0'). A leading '-' is used for negative
2726 * values. Hexadecimal digits are uppercase.
2728 public string ToHexString()
2730 if (varray == null && small == 0) {
2733 StringBuilder sb = new StringBuilder();
2739 if (x.varray == null) {
2740 AppendHex(sb, (uint)x.small, true);
2742 int n = x.varray.Length;
2743 AppendHex(sb, x.varray[n - 1], true);
2744 for (int j = n - 2; j >= 0; j --) {
2745 AppendHex(sb, x.varray[j], false);
2748 return sb.ToString();
2752 * Convert this value to decimal. A leading '-' sign is used for
2753 * negative value. The number of digits is minimal (no leading '0',
2754 * except for zero, which is returned as "0").
2756 public override string ToString()
2758 return ToString(10);
2761 private static int[] NDIGITS32 = {
2762 0, 0, 31, 20, 15, 13, 12, 11, 10, 10, 9, 9, 8, 8, 8, 8, 7, 7,
2763 7, 7, 7, 7, 7, 7, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6
2766 private static uint[] RADIXPP32 = {
2767 0, 0, 2147483648, 3486784401, 1073741824, 1220703125,
2768 2176782336, 1977326743, 1073741824, 3486784401, 1000000000,
2769 2357947691, 429981696, 815730721, 1475789056, 2562890625,
2770 268435456, 410338673, 612220032, 893871739, 1280000000,
2771 1801088541, 2494357888, 3404825447, 191102976, 244140625,
2772 308915776, 387420489, 481890304, 594823321, 729000000,
2773 887503681, 1073741824, 1291467969, 1544804416, 1838265625,
2777 private static char[] DCHAR =
2778 "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ".ToCharArray();
2780 static void AppendDigits(StringBuilder sb, uint v, int radix, int num)
2782 while (num -- > 0) {
2783 sb.Append(DCHAR[v % (uint)radix]);
2784 v = v / (uint)radix;
2788 static void AppendDigits(StringBuilder sb, uint v, int radix)
2791 sb.Append(DCHAR[v % (uint)radix]);
2792 v = v / (uint)radix;
2797 * Convert this value to a string, in the provided radix. The
2798 * radix must be in the 2 to 36 range (inclusive); uppercase
2799 * letters 'A' to 'Z' are used for digits of value 10 to 35.
2800 * If the value is zero, then "0" is returned; otherwise, leading
2801 * '0' digits are removed. A leading '-' sign is added for negative
2804 public string ToString(int radix)
2806 if (radix < 2 || radix > 36) {
2807 throw new ArgumentOutOfRangeException();
2811 * Special optimized case for base 16.
2814 return ToHexString();
2825 StringBuilder sb = new StringBuilder();
2826 if (x.varray == null) {
2827 AppendDigits(sb, (uint)x.small, radix);
2829 uint[] m = new uint[x.varray.Length];
2830 Array.Copy(x.varray, 0, m, 0, m.Length);
2831 uint prad = RADIXPP32[radix];
2832 int dnum = NDIGITS32[radix];
2834 uint v = MutateDivRem(m, prad);
2835 bool qz = Length(0, m) == 0;
2837 AppendDigits(sb, v, radix);
2840 AppendDigits(sb, v, radix, dnum);
2850 static string Reverse(StringBuilder sb)
2853 char[] tc = new char[n];
2854 sb.CopyTo(0, tc, 0, n);
2855 for (int i = 0, j = n - 1; i < j; i ++, j --) {
2860 return new string(tc);
2863 static uint DigitValue(char c, int radix)
2866 if (c >= '0' && c <= '9') {
2868 } else if (c >= 'a' && c <= 'z') {
2870 } else if (c >= 'A' && c <= 'Z') {
2875 if (d < 0 || d >= radix) {
2876 throw new ArgumentException();
2881 static ZInt ParseUnsigned(string s, int radix)
2883 if (s.Length == 0) {
2884 throw new ArgumentException();
2889 uint prad = RADIXPP32[radix];
2890 int dnum = NDIGITS32[radix];
2891 foreach (char c in s) {
2892 uint d = DigitValue(c, radix);
2893 acc = acc * (uint)radix + d;
2894 if (++ accNum == dnum) {
2895 x = x * (ZInt)prad + (ZInt)acc;
2902 while (accNum -- > 0) {
2905 x = x * (ZInt)p + (ZInt)acc;
2912 * -- A leading '-' is allowed, to denote a negative value.
2913 * -- If there is a "0b" or "0B" header (after any '-' sign),
2914 * then the value is interpreted in base 2.
2915 * -- If there is a "0x" or "0X" header (after any '-' sign),
2916 * then the value is interpreted in base 16 (hexadecimal).
2917 * Both uppercase and lowercase letters are accepted.
2918 * -- If there is no header, then decimal interpretation is used.
2920 * Unexpected characters (including spaces) trigger exceptions.
2921 * There must be at least one digit.
2923 public static ZInt Parse(string s)
2927 if (s.StartsWith("-")) {
2932 if (s.StartsWith("0b") || s.StartsWith("0B")) {
2935 } else if (s.StartsWith("0x") || s.StartsWith("0X")) {
2941 ZInt x = ParseUnsigned(s, radix);
2942 return neg ? -x : x;
2946 * Parse a string in the specified radix. The radix must be in
2947 * the 2 to 36 range (inclusive). Uppercase and lowercase letters
2948 * are accepted for digits in the 10 to 35 range.
2950 * A leading '-' sign is allowed, to denote a negative value.
2951 * Otherwise, only digits (acceptable with regards to the radix)
2952 * may appear. There must be at least one digit.
2954 public static ZInt Parse(string s, int radix)
2956 if (radix < 2 || radix > 36) {
2957 throw new ArgumentOutOfRangeException();
2961 if (s.StartsWith("-")) {
2965 ZInt x = ParseUnsigned(s, radix);
2966 return neg ? -x : x;
2969 static uint DecU32BE(byte[] buf, int off, int len)
2977 return ((uint)buf[off] << 8)
2978 | (uint)buf[off + 1];
2980 return ((uint)buf[off] << 16)
2981 | ((uint)buf[off + 1] << 8)
2982 | (uint)buf[off + 2];
2984 return ((uint)buf[off] << 24)
2985 | ((uint)buf[off + 1] << 16)
2986 | ((uint)buf[off + 2] << 8)
2987 | (uint)buf[off + 3];
2992 * Decode an integer, assuming unsigned big-endian encoding.
2993 * An empty array is decoded as 0.
2995 public static ZInt DecodeUnsignedBE(byte[] buf)
2997 return DecodeUnsignedBE(buf, 0, buf.Length);
3001 * Decode an integer, assuming unsigned big-endian encoding.
3002 * An empty array is decoded as 0.
3004 public static ZInt DecodeUnsignedBE(byte[] buf, int off, int len)
3006 while (len > 0 && buf[off] == 0) {
3012 } else if (len <= 4) {
3013 return new ZInt(DecU32BE(buf, off, len));
3015 uint[] m = new uint[(len + 3) >> 2];
3017 for (int j = len; j > 0; j -= 4) {
3021 w = DecU32BE(buf, off, j);
3023 w = DecU32BE(buf, off + k, 4);
3027 return new ZInt(0, m);
3030 static RNGCryptoServiceProvider RNG = new RNGCryptoServiceProvider();
3033 * Create a random integer of the provided size. Returned value
3034 * is in the 0 (inclusive) to 2^size (exclusive) range. A
3035 * cryptographically strong RNG is used to ensure uniform selection.
3037 public static ZInt MakeRand(int size)
3040 throw new ArgumentOutOfRangeException();
3042 byte[] buf = new byte[(size + 7) >> 3];
3046 buf[0] &= (byte)(0xFF >> (8 - kb));
3048 return DecodeUnsignedBE(buf);
3052 * Create a random integer in the 0 (inclusive) to max (exclusive)
3053 * range. 'max' must be positive. A cryptographically strong RNG
3054 * is used to ensure uniform selection.
3056 public static ZInt MakeRand(ZInt max)
3058 if (max.Sign <= 0) {
3059 throw new ArgumentOutOfRangeException();
3061 int bl = max.BitLength;
3063 ZInt x = MakeRand(bl);
3071 * Create a random integer in the min (inclusive) to max (exclusive)
3072 * range. 'max' must be greater than min. A cryptographically
3073 * strong RNG is used to ensure uniform selection.
3075 public static ZInt MakeRand(ZInt min, ZInt max)
3078 throw new ArgumentOutOfRangeException();
3080 return min + MakeRand(max - min);
3084 * Check whether this integer is prime. A probabilistic algorithm
3085 * is used, that theoretically ensures that a non-prime won't be
3086 * declared prime with probability greater than 2^(-100). Note that
3087 * this holds regardless of how the integer was generated (this
3088 * method does not assume that uniform random selection was used).
3090 * (Realistically, the probability of a computer hardware
3091 * malfunction is way greater than 2^(-100), so this property
3092 * returns the primality status with as much certainty as can be
3093 * achieved with a computer.)
3095 public bool IsPrime {
3097 return IsProbablePrime(50);
3101 static uint[] PRIMES_BF = new uint[] {
3102 0xA08A28AC, 0x28208A20, 0x02088288, 0x800228A2,
3103 0x20A00A08, 0x80282088, 0x800800A2, 0x08028228,
3104 0x0A20A082, 0x22880020, 0x28020800, 0x88208082,
3105 0x02022020, 0x08828028, 0x8008A202, 0x20880880
3108 private bool IsProbablePrime(int rounds)
3114 } else if (cc < 0) {
3117 if (x.varray == null) {
3118 if (x.small < (PRIMES_BF.Length << 5)) {
3119 return (PRIMES_BF[x.small >> 5]
3120 & ((uint)1 << (x.small & 31))) != 0;
3123 if (!x.TestBit(0)) {
3131 for (a = 0; !m.TestBit(a); a ++);
3133 while (rounds -- > 0) {
3134 ZInt b = MakeRand(Two, x);
3135 ZInt z = ModPow(b, m, x);
3136 for (int j = 0; j < a; j ++) {
3157 * Encode this integer as bytes (signed big-endian convention).
3158 * Encoding is of minimal length that still contains a sign bit
3159 * (compatible with ASN.1 DER encoding).
3161 public byte[] ToBytesBE()
3163 byte[] r = new byte[(BitLength + 8) >> 3];
3164 ToBytesBE(r, 0, r.Length);
3169 * Encode this integer as bytes (signed little-endian convention).
3170 * Encoding is of minimal length that still contains a sign bit.
3172 public byte[] ToBytesLE()
3174 byte[] r = new byte[(BitLength + 8) >> 3];
3175 ToBytesLE(r, 0, r.Length);
3180 * Encode this integer as bytes (signed big-endian convention).
3181 * Output length is provided; exactly that many bytes will be
3182 * written. The value is sign-extended or truncated if needed.
3184 public void ToBytesBE(byte[] buf, int off, int len)
3186 ToBytes(true, buf, off, len);
3190 * Encode this integer as bytes (signed little-endian convention).
3191 * Output length is provided; exactly that many bytes will be
3192 * written. The value is sign-extended or truncated if needed.
3194 public void ToBytesLE(byte[] buf, int off, int len)
3196 ToBytes(false, buf, off, len);
3200 * Encode this integer as bytes (unsigned big-endian convention).
3201 * Encoding is of minimal length, possibly without a sign bit. If
3202 * this value is zero, then an empty array is returned. If this
3203 * value is negative, then an ArgumentOutOfRangeException is thrown.
3205 public byte[] ToBytesUnsignedBE()
3208 throw new ArgumentOutOfRangeException();
3210 byte[] r = new byte[(BitLength + 7) >> 3];
3211 ToBytesBE(r, 0, r.Length);
3216 * Encode this integer as bytes (unsigned little-endian convention).
3217 * Encoding is of minimal length, possibly without a sign bit. If
3218 * this value is zero, then an empty array is returned. If this
3219 * value is negative, then an ArgumentOutOfRangeException is thrown.
3221 public byte[] ToBytesUnsignedLE()
3224 throw new ArgumentOutOfRangeException();
3226 byte[] r = new byte[(BitLength + 7) >> 3];
3227 ToBytesLE(r, 0, r.Length);
3231 void ToBytes(bool be, byte[] buf, int off, int len)
3233 uint iw = (uint)small >> 31;
3234 for (int i = 0; i < len; i ++) {
3237 if (varray == null) {
3238 w = (j == 0) ? (uint)small : iw;
3240 w = (j < varray.Length) ? varray[j] : iw;
3242 byte v = (byte)(w >> ((i & 3) << 3));
3244 buf[off + len - 1 - i] = v;