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
30 * Mutable container for a modular integer.
34 * - Each instance is initialised over a given modulus, or by
35 * duplicating an existing instance.
37 * - All operands for a given operation must share the same modulus,
38 * i.e. the instances must have been created with Dup() calls from
41 * - Modulus must be odd and greater than 1.
47 * Dedicated class for a modulus. It maintains some useful
48 * values that depend only on the modulus.
53 * val = modulus value, 31 bits per word, little-endian
54 * bitLen = modulus bit length (exact)
55 * n0i is such that n0i * val[0] == -1 mod 2^31
56 * R = 2^(31*val.Length) (modular)
57 * R2 = 2^(31*(val.Length*2)) (modular)
58 * Rx[i] = 2^(31*(val.Length+(2^i))) (modular)
68 internal ZMod(byte[] bmod, int off, int len)
70 bitLen = BigInt.BitLength(bmod, off, len);
71 if (bitLen <= 1 || (bmod[off + len - 1] & 1) == 0) {
72 throw new CryptoException("invalid modulus");
74 int n = (bitLen + 30) / 31;
76 DecodeBE(bmod, off, len, val);
83 n0i = (1 + ~y) & 0x7FFFFFFF;
86 * Compute the Rx[] values.
89 while ((n >> zk) != 0) {
94 for (int i = 0; i < 31; i ++) {
99 Array.Copy(R, 0, Rx[0], 0, n);
100 for (int i = 0; i < 31; i ++) {
103 for (int k = 1; k < zk; k ++) {
105 MontyMul(Rx[k - 1], Rx[k - 1], Rx[k], val, n0i);
109 * Compute R2 by multiplying the relevant Rx[]
113 uint[] tt = new uint[n];
114 for (int k = 0; k < zk; k ++) {
115 if (((n >> k) & 1) == 0) {
120 Array.Copy(Rx[k], 0, R2, 0, n);
122 MontyMul(Rx[k], R2, tt, val, n0i);
123 Array.Copy(tt, 0, R2, 0, n);
128 * Compute N-2; used as an exponent for modular
129 * inversion. Since modulus N is odd and greater
130 * than 1, N-2 is necessarily positive.
132 vm2 = new byte[(bitLen + 7) >> 3];
134 for (int i = 0; i < vm2.Length; i ++) {
135 int b = bmod[off + len - 1 - i];
137 vm2[vm2.Length - 1 - i] = (byte)b;
148 * Get the modulus bit length.
150 public int ModBitLength {
157 * Test whether this value is zero.
161 return IsZeroCT != 0;
165 public uint IsZeroCT {
168 for (int i = 1; i < val.Length; i ++) {
171 return ~(uint)((z | -z) >> 31);
176 * Test whether this value is one.
184 public uint IsOneCT {
186 int z = (int)val[0] ^ 1;
187 for (int i = 1; i < val.Length; i ++) {
190 return ~(uint)((z | -z) >> 31);
200 * Create a new instance by decoding the provided modulus
201 * (unsigned big-endian). Value is zero.
203 public ModInt(byte[] modulus)
204 : this(modulus, 0, modulus.Length)
209 * Create a new instance by decoding the provided modulus
210 * (unsigned big-endian). Value is zero.
212 public ModInt(byte[] modulus, int off, int len)
214 Init(new ZMod(modulus, off, len));
220 int n = mod.val.Length;
227 * Duplicate this instance. The new instance uses the same
228 * modulus, and its value is initialized to the same value as
233 ModInt m = new ModInt(mod);
234 Array.Copy(val, 0, m.val, 0, val.Length);
239 * Set the value in this instance to a copy of the value in
240 * the provided instance. The 'm' instance may use a different
241 * modulus, in which case the value may incur modular reduction.
243 public void Set(ModInt m)
246 * If the instances use the same modulus array, then
247 * the value can simply be copied as is.
250 Array.Copy(m.val, 0, val, 0, val.Length);
257 * Set the value in this instance to a copy of the value in
258 * the provided instance, but do so only conditionally. If
259 * 'ctl' is -1, then the copy is done. If 'ctl' is 0, then the
260 * copy is NOT done. This instance and the source instance
261 * MUST use the same modulus.
263 public void CondCopy(ModInt m, uint ctl)
266 throw new CryptoException("Not same modulus");
268 for (int i = 0; i < val.Length; i ++) {
270 val[i] = w ^ (ctl & (m.val[i] ^ w));
275 * Set the value in this instance to a copy of either 'a'
276 * (if ctl is -1) or 'b' (if ctl is 0).
278 public void CopyMux(uint ctl, ModInt a, ModInt b)
280 if (mod != a.mod || mod != b.mod) {
281 throw new CryptoException("Not same modulus");
283 for (int i = 0; i < val.Length; i ++) {
286 val[i] = bw ^ (ctl & (aw ^ bw));
291 * Set the value in this instance to the provided integer value
292 * (which MUST be nonnegative).
294 public void Set(int x)
297 for (int i = 1; i < val.Length; i ++) {
303 * Set this value to either 0 (if ctl is 0) or 1 (if ctl is -1),
304 * in Montgomery representation.
306 public void SetMonty(uint ctl)
308 for (int i = 0; i < val.Length; i ++) {
309 val[i] = ctl & mod.R[i];
314 * Set this value by decoding the provided integer (big-endian
315 * unsigned encoding). If the source value does not fit (it is
316 * not lower than the modulus), then this method return 0, and
317 * this instance is set to 0. Otherwise, it returns -1.
319 public uint Decode(byte[] buf)
321 return Decode(buf, 0, buf.Length);
325 * Set this value by decoding the provided integer (big-endian
326 * unsigned encoding). If the source value does not fit (it is
327 * not lower than the modulus), then this method return 0, and
328 * this instance is set to 0. Otherwise, it returns -1.
330 public uint Decode(byte[] buf, int off, int len)
333 * Decode into val[]; if the truncation loses some
334 * non-zero bytes, then this returns 0.
336 uint x = DecodeBE(buf, off, len, val);
339 * Compare with modulus. We want to get a 0 if the
340 * subtraction would not yield a carry, i.e. the
341 * source value is not lower than the modulus.
343 x &= Sub(val, mod.val, 0);
346 * At that point, x is -1 if the value fits, 0
349 for (int i = 0; i < val.Length; i ++) {
356 * Set this value by decoding the provided integer (big-endian
357 * unsigned encoding). The source value is reduced if necessary.
359 public void DecodeReduce(byte[] buf)
361 DecodeReduce(buf, 0, buf.Length);
365 * Set this value by decoding the provided integer (big-endian
366 * unsigned encoding). The source value is reduced if necessary.
368 public void DecodeReduce(byte[] buf, int off, int len)
370 uint[] x = new uint[((len << 3) + 30) / 31];
371 DecodeBE(buf, off, len, x);
376 * Set the value from the provided source array, with modular
379 void Reduce(uint[] b)
382 * If the modulus uses only one word then we must use
383 * a special code path.
385 if (mod.val.Length == 1) {
391 * Fast copy of words that do not incur modular
394 int aLen = mod.val.Length;
396 int cLen = Math.Min(aLen - 1, bLen);
397 Array.Copy(b, bLen - cLen, val, 0, cLen);
398 for (int i = cLen; i < aLen; i ++) {
403 * Inject extra words. We use the pre-computed
404 * Rx[] values to do shifts.
406 for (int j = bLen - cLen; j > 0;) {
408 * We can add power-of-2 words, but less
409 * than the modulus size. Note that the modulus
410 * uses at least two words, so this process works.
414 int nk = 1 << (k + 1);
415 if (nk >= aLen || nk > j) {
420 MontyMul(val, mod.Rx[k], tmp1, mod.val, mod.n0i);
422 Array.Copy(b, j, tmp2, 0, num);
423 for (int i = num; i < tmp2.Length; i ++) {
426 ModAdd(tmp1, tmp2, val, mod.val, 0);
431 * Modular reduction in case the modulus fits on a single
434 void ReduceSmallMod(uint[] b)
438 int nlen = mod.bitLen;
441 for (int i = b.Length - 1; i >= 0; i --) {
443 * Multiply x by R (Montgomery multiplication by R^2).
445 ulong z = (ulong)x * (ulong)r2;
446 uint u = ((uint)z * n0i) & 0x7FFFFFFF;
447 z += (ulong)u * (ulong)n;
451 * Ensure x fits on 31 bits (it may be up to twice
452 * the modulus at that point). If x >= 2^31 then,
453 * necessarily, x is greater than the modulus, and
454 * the subtraction is sound; moreover, in that case,
455 * subtracting the modulus brings back x to less
456 * than the modulus, hence fitting on 31 bits.
458 x -= (uint)((int)x >> 31) & n;
461 * Add the next word, then reduce. The addition
462 * does not overflow since both operands fit on
465 * Since the modulus could be much smaller than
466 * 31 bits, we need a full remainder operation here.
471 * Constant-time modular reduction.
472 * We first perform two subtraction of the
473 * shifted modulus to ensure that the high bit
474 * is cleared. This allows the loop to work
477 x -= (uint)((int)x >> 31) & (n << (31 - nlen));
478 x -= (uint)((int)x >> 31) & (n << (31 - nlen));
479 for (int j = 31 - nlen; j >= 0; j --) {
481 x += (uint)((int)x >> 31) & (n << j);
488 * Encode into bytes. Big-endian unsigned encoding is used, the
489 * returned array having the minimal length to encode the modulus.
491 public byte[] Encode()
493 return Encode(false);
497 * Encode into bytes. Big-endian encoding is used; if 'signed' is
498 * true, then signed encoding is used: returned value will have a
499 * leading bit set to 0. Returned array length is the minimal size
500 * for encoding the modulus (with a sign bit if using signed
503 public byte[] Encode(bool signed)
509 byte[] buf = new byte[(x + 7) >> 3];
510 Encode(buf, 0, buf.Length);
515 * Encode into bytes. The provided array is fully set; big-endian
516 * encoding is used, and extra leading bytes of value 0 are added
517 * if necessary. If the destination array is too small, then the
518 * value is silently truncated.
520 public void Encode(byte[] buf)
522 Encode(buf, 0, buf.Length);
526 * Encode into bytes. The provided array chunk is fully set;
527 * big-endian encoding is used, and extra leading bytes of value
528 * 0 are added if necessary. If the destination array is too
529 * small, then the value is silently truncated.
531 public void Encode(byte[] buf, int off, int len)
533 EncodeBE(val, buf, off, len);
537 * Get the least significant bit of the value (0 or 1).
541 return val[0] & (uint)1;
545 * Add a small integer to this instance. The small integer
546 * 'x' MUST be lower than 2^31 and MUST be lower than the modulus.
548 public void Add(uint x)
551 for (int i = 1; i < tmp1.Length; i ++) {
554 ModAdd(val, tmp1, val, mod.val, 0);
558 * Add another value to this instance. The operand 'b' may
559 * be the same as this instance.
561 public void Add(ModInt b)
564 throw new CryptoException("Not same modulus");
566 ModAdd(val, b.val, val, mod.val, 0);
570 * Subtract a small integer from this instance. The small integer
571 * 'x' MUST be lower than 2^31 and MUST be lower than the modulus.
573 public void Sub(uint x)
576 for (int i = 1; i < tmp1.Length; i ++) {
579 ModSub(val, tmp1, val, mod.val);
583 * Subtract another value from this instance. The operand 'b'
584 * may be the same as this instance.
586 public void Sub(ModInt b)
589 throw new CryptoException("Not same modulus");
591 ModSub(val, b.val, val, mod.val);
599 ModSub(null, val, val, mod.val);
603 * Convert this instance to Montgomery representation.
605 public void ToMonty()
607 MontyMul(val, mod.R2, tmp1, mod.val, mod.n0i);
608 Array.Copy(tmp1, 0, val, 0, val.Length);
612 * Convert this instance back from Montgomery representation to
613 * normal representation.
615 public void FromMonty()
618 for (int i = 1; i < tmp1.Length; i ++) {
621 MontyMul(val, tmp1, tmp2, mod.val, mod.n0i);
622 Array.Copy(tmp2, 0, val, 0, val.Length);
626 * Compute a Montgomery multiplication with the provided
627 * value. The other operand may be this instance.
629 public void MontyMul(ModInt b)
632 throw new CryptoException("Not same modulus");
634 MontyMul(val, b.val, tmp1, mod.val, mod.n0i);
635 Array.Copy(tmp1, 0, val, 0, val.Length);
639 * Montgomery-square this instance.
641 public void MontySquare()
643 MontyMul(val, val, tmp1, mod.val, mod.n0i);
644 Array.Copy(tmp1, 0, val, 0, val.Length);
648 * Perform modular exponentiation. Exponent is in big-endian
651 public void Pow(byte[] exp)
653 Pow(exp, 0, exp.Length);
657 * Perform modular exponentiation. Exponent is in big-endian
660 public void Pow(byte[] exp, int off, int len)
662 MontyMul(val, mod.R2, tmp1, mod.val, mod.n0i);
664 for (int i = 1; i < val.Length; i ++) {
667 for (int i = 0; i < len; i ++) {
668 int x = exp[off + len - 1 - i];
669 for (int j = 0; j < 8; j ++) {
670 MontyMul(val, tmp1, tmp2, mod.val, mod.n0i);
671 uint ctl = (uint)-((x >> j) & 1);
672 for (int k = 0; k < val.Length; k ++) {
673 val[k] = (tmp2[k] & ctl)
676 MontyMul(tmp1, tmp1, tmp2, mod.val, mod.n0i);
677 Array.Copy(tmp2, 0, tmp1, 0, tmp2.Length);
683 * Compute modular inverse of this value. If this instance is
684 * zero, then it remains equal to zero. If the modulus is not
685 * prime, then this function computes wrong values.
693 * Compute the square root for this value. This method assumes
694 * that the modulus is prime, greater than or equal to 7, and
695 * equal to 3 modulo 4; if it is not, then the returned value
696 * and the contents of this instance are indeterminate.
698 * Returned value is -1 if the value was indeed a square, 0
699 * otherwise. In the latter case, array contents are the square
700 * root of the opposite of the original value.
702 public uint SqrtBlum()
705 * We suppose that p = 3 mod 4; we raise to the power
706 * (p-3)/4, then do an extra multiplication to go to
709 * Since we know the modulus bit length, we can do
710 * the exponentiation from the high bits downwards.
713 Array.Copy(val, 0, tmp1, 0, val.Length);
714 int k = (mod.bitLen - 2) / 31;
715 int j = mod.bitLen - 2 - k * 31;
716 uint ew = mod.val[k];
717 for (int i = mod.bitLen - 2; i >= 2; i --) {
718 uint ctl = ~(uint)-((int)(ew >> j) & 1);
719 MontyMul(tmp1, tmp1, tmp2, mod.val, mod.n0i);
720 MontyMul(val, tmp2, tmp1, mod.val, mod.n0i);
721 for (int m = 0; m < tmp1.Length; m ++) {
723 tmp1[m] = w ^ (ctl & (w ^ tmp2[m]));
732 * The extra multiplication. Square root is written in
733 * tmp2 (in Montgomery representation).
735 MontyMul(val, tmp1, tmp2, mod.val, mod.n0i);
738 * Square it back in tmp1, to see if it indeed yields
741 MontyMul(tmp2, tmp2, tmp1, mod.val, mod.n0i);
743 for (int i = 0; i < val.Length; i ++) {
744 z |= (int)(val[i] ^ tmp1[i]);
746 uint good = ~(uint)((z | -z) >> 31);
749 * Convert back the result to normal representation.
751 Array.Copy(tmp2, 0, val, 0, val.Length);
757 * Conditionally swap this instance with the one provided in
758 * parameter. The swap is performed if ctl is -1, not performed
761 public void CondSwap(ModInt b, uint ctl)
764 throw new CryptoException("Not same modulus");
766 for (int i = 0; i < val.Length; i ++) {
769 uint m = ctl & (x ^ y);
776 * Compare for equality this value with another. Comparison still
777 * works if the two values use distinct moduli.
779 public bool Eq(ModInt b)
785 * Compare for equality this value with another. Comparison still
786 * works if the two values use distinct moduli. Returned value is
787 * -1 on equality, 0 otherwise.
789 public uint EqCT(ModInt b)
792 if (b.val.Length > val.Length) {
793 for (int i = 0; i < val.Length; i ++) {
794 z |= val[i] ^ b.val[i];
796 for (int i = val.Length; i < b.val.Length; i ++) {
800 for (int i = 0; i < b.val.Length; i ++) {
801 z |= val[i] ^ b.val[i];
803 for (int i = b.val.Length; i < val.Length; i ++) {
808 return ~(uint)((x | -x) >> 31);
811 /* ============================================================== */
814 * Decode value (unsigned big-endian) into dst[] (little-endian,
815 * 31 bits per word). All words of dst[] are initialised.
816 * Returned value is -1 on success, 0 if some non-zero source
817 * bits had to be ignored.
819 static uint DecodeBE(byte[] buf, int off, int len, uint[] dst)
825 while (i < dst.Length) {
833 acc |= (b << accLen);
836 dst[i ++] = acc & 0x7FFFFFFF;
838 acc = b >> (8 - accLen);
845 return ~(uint)((x | -x) >> 31);
849 * Encode an integer (array of words, little-endian, 31 bits per
850 * word) into big-endian encoding, with the provided length (in
853 static byte[] EncodeBE(uint[] x, int len)
855 byte[] val = new byte[len];
856 EncodeBE(x, val, 0, len);
861 * Encode an integer (array of words, little-endian, 31 bits per
862 * word) into big-endian encoding, with the provided length (in
865 static void EncodeBE(uint[] x, byte[] val)
867 EncodeBE(x, val, 0, val.Length);
871 * Encode an integer (array of words, little-endian, 31 bits per
872 * word) into big-endian encoding, with the provided length (in
875 static void EncodeBE(uint[] x, byte[] val, int off, int len)
880 for (int i = len - 1; i >= 0; i --) {
883 uint z = (j < x.Length) ? x[j ++] : 0;
884 b = acc | (z << accLen);
885 acc = z >> (8 - accLen);
892 val[off + i] = (byte)b;
897 * Subtract b from a; the carry is returned (-1 if carry, 0
898 * otherwise). The operation is done only if ctl is -1; if
899 * ctl is 0, then a[] is unmodified, but the carry is still
900 * computed and returned.
902 * The two operand arrays MUST have the same size.
904 static uint Sub(uint[] a, uint[] b, uint ctl)
909 for (int i = 0; i < n; i ++) {
912 uint cw = (uint)cc + aw - bw;
914 a[i] = (aw & ~ctl) | (cw & ctl);
920 * Left-shift value by one bit; the extra bit (carry) is
923 static uint LShift(uint[] a)
927 for (int i = 0; i < n; i ++) {
929 a[i] = (cc | (aw << 1)) & 0x7FFFFFFF;
932 return (uint)-(int)cc;
936 * Modular left-shift value by one bit. Value and modulus MUST
937 * have the same length. Value MUST be lower than modulus.
939 static void ModMul2(uint[] a, uint[] mod)
944 * First pass: compute 2*a-mod, but don't keep the
945 * result, only the final carry (0 or -1).
949 for (int i = 0; i < n; i ++) {
951 uint aws = ((aw << 1) | cc1) & 0x7FFFFFFF;
953 uint z = aws - mod[i] + (uint)cc2;
959 * If cc2 is 0, then the subtraction yields no carry and
960 * must be done. Otherwise, cc2 is -1, and the subtraction
963 uint ctl = ~(uint)cc2;
966 for (int i = 0; i < n; i ++) {
968 uint aws = ((aw << 1) | cc1) & 0x7FFFFFFF;
970 uint z = aws - (mod[i] & ctl) + (uint)cc2;
972 a[i] = z & 0x7FFFFFFF;
979 * If 'hi' is zero, then this computes a+b-n if a+b >= n,
982 * If 'hi' is non-zero, then this computes 2^k+a+b-n, where
985 * Result is written in d[]. The same array may be used in
988 static void ModAdd(uint[] a, uint[] b, uint[] d, uint[] n, uint hi)
991 * Set ctl to -1 if hi is non-zero, 0 otherwise.
992 * 'ctl' computes whether the subtraction with n[]
993 * is needed in the second pass.
996 uint ctl = (uint)((x | -x) >> 31);
998 for (int pass = 0; pass < 2; pass ++) {
1000 * cc1 is the carry for a+b (0 or 1)
1002 * cc2 is the carry for the modulus
1003 * subtraction (0 or -1)
1007 for (int i = 0; i < n.Length; i ++) {
1011 uint sw = aw + bw + cc1;
1014 uint dw = sw - nw + (uint)cc2;
1015 cc2 = (int)dw >> 31;
1018 d[i] = (dw & ctl) | (sw & ~ctl);
1023 * Compute aggregate subtraction carry. This should
1024 * not be 1 if the operands are correct (it would
1025 * mean that a+b-n overflows, so both a and b are
1026 * larger than n). If it is 0, then a+b-n >= 0,
1027 * so the subtraction must be done; otherwise, the
1028 * aggregate carry will be -1, and the subtraction
1029 * should not be done (unless forced by a non-zero
1033 ctl |= ~(uint)(cc2 >> 31);
1038 * Modular subtraction.
1040 * This computes a-b, then adds n if the result is negative.
1041 * If a is null, then it is assumed to be all-zeros.
1043 * Result is written in d[]. The same array may be used in
1046 static void ModSub(uint[] a, uint[] b, uint[] d, uint[] n)
1049 for (int pass = 0; pass < 2; pass ++) {
1051 * cc1 = carry for a-b (0 or -1)
1052 * cc2 = carry for modulus addition (0 or 1)
1056 for (int i = 0; i < n.Length; i ++) {
1057 uint aw = (a == null) ? 0 : a[i];
1060 uint sw = aw - bw + (uint)cc1;
1061 cc1 = (int)sw >> 31;
1063 uint dw = sw + nw + cc2;
1067 d[i] = (dw & ctl) | (sw & ~ctl);
1072 * Modulus addition must be done if and only if
1073 * a-b had a final carry.
1080 * Compute Montgomery multiplication of a[] by b[], result in
1081 * d[], with modulus n[]. All arrays must have the same length.
1082 * d[] must be distinct from a[], b[] and n[]. Modulus must be
1083 * odd. n0i must be such that n[0]*n0i = -1 mod 2^31. Values a[]
1084 * and b[] must be lower than n[] (so, in particular, a[] and
1085 * b[] cannot be the same array as n[]).
1087 static void MontyMul(uint[] a, uint[] b, uint[] d, uint[] n, uint n0i)
1090 for (int i = 0; i < len; i ++) {
1094 for (int i = 0; i < len; i ++) {
1096 uint u = ((d[0] + ai * b[0]) * n0i) & 0x7FFFFFFF;
1098 for (int j = 0; j < len; j ++) {
1099 ulong z = (ulong)d[j]
1100 + (ulong)ai * (ulong)b[j]
1101 + (ulong)u * (ulong)n[j] + cc;
1104 d[j - 1] = (uint)z & 0x7FFFFFFF;
1107 if (((uint)z & 0x7FFFFFFF) != 0) {
1108 throw new Exception("BAD!");
1113 d[len - 1] = (uint)dh & 0x7FFFFFFF;
1117 uint ctl = (uint)((x | -x) >> 31);
1118 Sub(d, n, ctl | ~Sub(d, n, 0));