Commit 1c85a28372bdd010137a742e9af19f861ab4f06a
2015-04-17T22:46:11
adding bn_mp_sqrtmod_prime.c
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142
diff --git a/bn_mp_sqrtmod_prime.c b/bn_mp_sqrtmod_prime.c
new file mode 100755
index 0000000..1399fe0
--- /dev/null
+++ b/bn_mp_sqrtmod_prime.c
@@ -0,0 +1,122 @@
+#include <tommath.h>
+#ifdef BN_MP_SQRTMOD_PRIME_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis
+ *
+ * LibTomMath is a library that provides multiple-precision
+ * integer arithmetic as well as number theoretic functionality.
+ *
+ * The library is free for all purposes without any express
+ * guarantee it works.
+ */
+
+/* Tonelli-Shanks algorithm
+ * https://en.wikipedia.org/wiki/Tonelli%E2%80%93Shanks_algorithm
+ * https://gmplib.org/list-archives/gmp-discuss/2013-April/005300.html
+ *
+ */
+
+int mp_sqrtmod_prime(mp_int *n, mp_int *prime, mp_int *ret)
+{
+ int res, legendre;
+ mp_int t1, C, Q, S, Z, M, T, R, two;
+ unsigned long i;
+
+ /* first handle the simple cases */
+ if (mp_cmp_d(n, 0) == MP_EQ) {
+ mp_zero(ret);
+ return MP_OKAY;
+ }
+ if (mp_cmp_d(prime, 2) == MP_EQ) return MP_VAL; /* prime must be odd */
+ if ((res = mp_jacobi(n, prime, &legendre)) != MP_OKAY) return res;
+ if (legendre == -1) return MP_VAL; /* quadratic non-residue mod prime */
+
+ mp_init_multi(&t1, &C, &Q, &S, &Z, &M, &T, &R, &two, NULL);
+
+ /* SPECIAL CASE: if prime mod 4 == 3
+ * compute directly: res = n^(prime+1)/4 mod prime
+ * Handbook of Applied Cryptography algorithm 3.36
+ */
+ if ((res = mp_mod_d(prime, 4, &i)) != MP_OKAY) goto cleanup;
+ if (i == 3) {
+ if ((res = mp_add_d(prime, 1, &t1)) != MP_OKAY) goto cleanup;
+ if ((res = mp_div_2(&t1, &t1)) != MP_OKAY) goto cleanup;
+ if ((res = mp_div_2(&t1, &t1)) != MP_OKAY) goto cleanup;
+ if ((res = mp_exptmod(n, &t1, prime, ret)) != MP_OKAY) goto cleanup;
+ res = MP_OKAY;
+ goto cleanup;
+ }
+
+ /* NOW: Tonelli-Shanks algorithm */
+
+ /* factor out powers of 2 from prime-1, defining Q and S as: prime-1 = Q*2^S */
+ if ((res = mp_copy(prime, &Q)) != MP_OKAY) goto cleanup;
+ if ((res = mp_sub_d(&Q, 1, &Q)) != MP_OKAY) goto cleanup;
+ /* Q = prime - 1 */
+ mp_zero(&S);
+ /* S = 0 */
+ while (mp_iseven(&Q)) {
+ if ((res = mp_div_2(&Q, &Q)) != MP_OKAY) goto cleanup;
+ /* Q = Q / 2 */
+ if ((res = mp_add_d(&S, 1, &S)) != MP_OKAY) goto cleanup;
+ /* S = S + 1 */
+ }
+
+ /* find a Z such that the Legendre symbol (Z|prime) == -1 */
+ mp_set_int(&Z, 2);
+ /* Z = 2 */
+ while(1) {
+ if ((res = mp_jacobi(&Z, prime, &legendre)) != MP_OKAY) goto cleanup;
+ if (legendre == -1) break;
+ if ((res = mp_add_d(&Z, 1, &Z)) != MP_OKAY) goto cleanup;
+ /* Z = Z + 1 */
+ }
+
+ if ((res = mp_exptmod(&Z, &Q, prime, &C)) != MP_OKAY) goto cleanup;
+ /* C = Z ^ Q mod prime */
+ if ((res = mp_add_d(&Q, 1, &t1)) != MP_OKAY) goto cleanup;
+ if ((res = mp_div_2(&t1, &t1)) != MP_OKAY) goto cleanup;
+ /* t1 = (Q + 1) / 2 */
+ if ((res = mp_exptmod(n, &t1, prime, &R)) != MP_OKAY) goto cleanup;
+ /* R = n ^ ((Q + 1) / 2) mod prime */
+ if ((res = mp_exptmod(n, &Q, prime, &T)) != MP_OKAY) goto cleanup;
+ /* T = n ^ Q mod prime */
+ if ((res = mp_copy(&S, &M)) != MP_OKAY) goto cleanup;
+ /* M = S */
+ if ((res = mp_set_int(&two, 2)) != MP_OKAY) goto cleanup;
+
+ res = MP_VAL;
+ while (1) {
+ if ((res = mp_copy(&T, &t1)) != MP_OKAY) goto cleanup;
+ i = 0;
+ while (1) {
+ if (mp_cmp_d(&t1, 1) == MP_EQ) break;
+ if ((res = mp_exptmod(&t1, &two, prime, &t1)) != MP_OKAY) goto cleanup;
+ i++;
+ }
+ if (i == 0) {
+ mp_copy(&R, ret);
+ res = MP_OKAY;
+ goto cleanup;
+ }
+ if ((res = mp_sub_d(&M, i, &t1)) != MP_OKAY) goto cleanup;
+ if ((res = mp_sub_d(&t1, 1, &t1)) != MP_OKAY) goto cleanup;
+ if ((res = mp_exptmod(&two, &t1, prime, &t1)) != MP_OKAY) goto cleanup;
+ /* t1 = 2 ^ (M - i - 1) */
+ if ((res = mp_exptmod(&C, &t1, prime, &t1)) != MP_OKAY) goto cleanup;
+ /* t1 = C ^ (2 ^ (M - i - 1)) mod prime */
+ if ((res = mp_sqrmod(&t1, prime, &C)) != MP_OKAY) goto cleanup;
+ /* C = (t1 * t1) mod prime */
+ if ((res = mp_mulmod(&R, &t1, prime, &R)) != MP_OKAY) goto cleanup;
+ /* R = (R * t1) mod prime */
+ if ((res = mp_mulmod(&T, &C, prime, &T)) != MP_OKAY) goto cleanup;
+ /* T = (T * C) mod prime */
+ mp_set(&M, i);
+ /* M = i */
+ }
+
+cleanup:
+ mp_clear_multi(&t1, &C, &Q, &S, &Z, &M, &T, &R, &two, NULL);
+ return res;
+}
+
+#endif
diff --git a/tommath.h b/tommath.h
index 3fa5c3f..fa69688 100644
--- a/tommath.h
+++ b/tommath.h
@@ -430,6 +430,9 @@ int mp_n_root_ex (mp_int * a, mp_digit b, mp_int * c, int fast);
/* special sqrt algo */
int mp_sqrt(mp_int *arg, mp_int *ret);
+/* special sqrt (mod prime) */
+int mp_sqrtmod_prime(mp_int *arg, mp_int *prime, mp_int *ret);
+
/* is number a square? */
int mp_is_square(mp_int *arg, int *ret);