new function mp_sqrtmod_prime

[karel-m]

Hi,

it would be nice to have mp_sqrtmod_prime function - it is necessary for handling compressed ECC keys as mentioned in this issue libtom/libtomcrypt#34

here is my suggestion (I'll send it as pull request if you find it worth considering):

diff --git a/src/ltm/bn_mp_sqrtmod_prime.c b/src/ltm/bn_mp_sqrtmod_prime.c
--- /dev/null
+++ b/src/ltm/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;
+
+  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 */
+  }
+
+  res = MP_VAL;
+cleanup:
+  mp_clear_multi(&t1, &C, &Q, &S, &Z, &M, &T, &R, &two, NULL);
+  return res;
+}
+
+#endif
diff --git a/src/ltm/tommath.h b/src/ltm/tommath.h
index 58a26c8..14da465 100644
--- a/src/ltm/tommath.h
+++ b/src/ltm/tommath.h
@@ -394,6 +394,9 @@ int mp_n_root(mp_int *a, mp_digit b, mp_int *c);
 /* 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);

diff --git a/src/ltm/tommath_class.h b/src/ltm/tommath_class.h
index 68b88b9..a82fa44 100644
--- a/src/ltm/tommath_class.h
+++ b/src/ltm/tommath_class.h
@@ -104,6 +104,7 @@
 #define BN_MP_SQR_C
 #define BN_MP_SQRMOD_C
 #define BN_MP_SQRT_C
+#define BN_MP_SQRTMOD_PRIME_C
 #define BN_MP_SUB_C
 #define BN_MP_SUB_D_C
 #define BN_MP_SUBMOD_C
@@ -823,6 +824,18 @@
    #define BN_MP_CLEAR_C
 #endif

+#if defined(BN_MP_SQRTMOD_PRIME_C)
+   #define BN_MP_JACOBI_C
+   #define BN_MP_ZERO_C
+   #define BN_MP_SET_INT_C
+   #define BN_MP_COPY_C
+   #define BN_MP_SUB_C
+   #define BN_MP_SUB_D_C
+   #define BN_MP_DIV_2_C
+   #define BN_MP_ADD_D_C
+   #define BN_MP_EXPTMOD_C
+#endif
+
 #if defined(BN_MP_SUB_C)
    #define BN_S_MP_ADD_C
    #define BN_MP_CMP_MAG_C

[tomstdenis]

Cool, care to write a bit in the BN user manual and/or the tommath text?

[sjaeckel]

+1 for user manual & tommath text
also, please create one commit with the changes in the c and tommath.h files and afterwards just do make new_file and commit these changes separately
then please open a PR

[tomstdenis]

To be fair I wasn't volunteering to do that work :-) (at least I won't be getting involved anytime soon beyond just chipping in my two cents). Should be easy to add to the BN manual and for the tommath book I would just add it towards the end near the discussion about Jacobi.

[sjaeckel]

Oh, I was pointing to @karel-m :)

[tomstdenis]

Phew, dodged work. +1 dexterity.

[sjaeckel]

I'm always trying to follow the approach 'the person who had the fun part to do the implementation, then please also takes the bad part and writes documentation' ;)

[tomstdenis]

Ya I played that role for 5+ years. I don't mind documentation/etc but I'm not in a place to contribute commits to OSS just yet and my 2 year old at home won't let me do it in my spare time :-)

That being said at the very least make sure to update the BN manual since as an SDK that's more important.

[sjaeckel]

I just don't merge PR's that introduce major features but miss documentation.

[tomstdenis]

Good lad. And yes, it would be good to see this ported to TFM for completeness.

[karel-m]

Ad creating pull request see #34

By doc update you mean a patch for bn.tex? If yes, I have to admit that "I don't speak TeX" but for sure I can write couple of paragraphs that somebody can put into the documentation.

[karel-m]

Here is some basic doc:

Modular square root

Function mp_sqrtmod_prime(n, p, x) solves the modular equatioon x^2 = n (mod p) where p is a prime number greater than 2 (odd prime). The result is returned in the third argument x, the function returns MP_OKAY on success, MP_VAL or another error on failure.

The implementation is split for two different cases:

1. if p mod 4 == 3 we apply Handbook of Applied Cryptography algorithm 3.36 and compute x directly as x = n^(p+1)/4 (mod p)
2. otherwise we use Tonelli-Shanks algorithm https://en.wikipedia.org/wiki/Tonelli–Shanks_algorithm

The function does not check the primality of parameter p thus it is up to the caller to assure that this parameter is a prime number. When p is a composite the function behaviour is undefined, it may even return a false-positive MP_OKAY.

[buggywhip]

On Apr 18, 2015, at 1:39 PM, karel-m notifications@github.com wrote:

Here is some basic doc:

Modular square root

Function mp_sqrtmod_prime(n, p, x) solves the modular equatioon x^2 = n (mod p) where p is a prime number greater than 2 (odd prime). The result is returned in the third argument x, the function returns MP_OKAY on success, MP_VAL or another error on failure.

The implementation is split for two different cases:

if p mod 4 == 3 we apply Handbook of Applied Cryptography algorithm 3.36 and compute x directly as x = n^(p+1)/4 (mod p)

otherwise we use Tonelli-Shanks algorithm https://en.wikipedia.org/wiki/Tonelli–Shanks_algorithm

The function does not check the primality of parameter p thus it is up to the caller to assure that this parameter is a prime number. When p is a composite the function behaviour is undefined, it may even return a false-positive MP_OKAY.

This last bothers me. Given that many/most prime numbers are only probabilistically prime, what are the downside problems caused by the number not actually being a prime? True, done right, the odds are small this would be a problem, but wouldn't it make sense to understand the risks of a bad answer? ...or somehow check the answer? ???

[karel-m]

The thing is that neither Handbook of Applied Cryptography algorithm 3.36 nor Tonelli-Shanks algorithm work for composite p and testing primality in each mp_sqrtmod_prime call is too expensive.

In fact I have not tested what exactly happens when p is not a prime. IMO it is not even guaranteed that the function will return an error. Therefore I propose to clearly state it in documentation.

[buggywhip]

On Apr 18, 2015, at 2:15 PM, karel-m notifications@github.com wrote:

testing primality in each mp_sqrtmod_prime call is too expensive

I need to get closer to this problem, but what about just squaring and comparing the answer?

[sjaeckel]

Do we have to do anything about this?
As the corresponding PR is already merged and a new release will be pushed-out soon we should act now.

[sjaeckel]

I just realized that this is related to #31 ... @czurnieden what do you think as you created #31 ? does the jacobi issue somehow influence proper functioning of this function?

[czurnieden]

I thought it was already resolved? It's two lines: one to check for 0/1 and return the proper result "1" and the other one in the docs to remark the lack of support for negative numerators with the only problem: legacy. I would do it myself but my version has gone quite off over the last year when I wasn't able to keep it up-to-date regularly. I think I'll follow Randall Munroe's advice he gave in http://xkcd.com/1597/ ;-)

But a bit more serious:

After careful consumption of a large can of coffee and a much smaller copy of Cohen's "Algebraic Number Theory" the answer is: No, it won't make much trouble.

If I understand it correct, the composites that cause the largest trouble are perfect powers k^n (with k,n in N+ and n>1)and which are quite expensive to test for because you need to do it brute force (there's a slightly faster algorithm listed in Cohen for p^n and p prime) and check for every prime up to log_2(k). N-th integer root is O(log n) as is exponentiation and we could take the table in is_small_prime for the primes to make the prime generation O(1) and ignore that part of the complexity calculation, but it won't make it cheap enough to "just slap it on". Doing sqrt(x)%p "by hand" would still be cheaper.

If you want to do it you should place the code (an implementation is in bn_mp_n_root.c in my fork) before the primality test(s) where it makes the most sense because all non-deterministic tests in use (Rabin-Miller and Lucas-Selfridge) can get in trouble with perfect powers.

I would just put a note in the docs to highlight the fact that the result cannot be trusted 100% if the primality of the input cannot be guaranteed.

@sjaeckel Yes, the documentation of my stuff is a work in progress. Although a very slow progress but it goes forward.

[tomstdenis]

Is this function "ready" and just needs enhancements for edge cases or is it considered defective?

[sjaeckel]

Should be ready AFAIK.
I implemented the handling of the edge case of negative values in mp_jacobi in e8d2609

[tomstdenis]

With that then let's close this issue.