mp_sqrtmod_prime can get caught in an infinite loop if "prime" isn't prime

[friedrichsenm]

If *prime isn't prime but is composite that is 1 mod 4, you can hit infinite loops in two places. (My use case only involves numbers 1 mod 4, so this may happen other times too)

The first is in lines 62-69

mp_set(&Z, 2uL);
   /* Z = 2 */
   for (;;) {
      if ((err = mp_kronecker(&Z, prime, &legendre)) != MP_OKAY)     goto LBL_END;
      if (legendre == -1) break;
      if ((err = mp_add_d(&Z, 1uL, &Z)) != MP_OKAY)               goto LBL_END;
      /* Z = Z + 1 */
   }

I was able to fix this in a local copy adding the following check after line 66

if (legendre == 0) {err = MP_VAL; goto LBL_END;}

The second infinite loop was in lines 84-111

for (;;) {
      if ((err = mp_copy(&T, &t1)) != MP_OKAY)                    goto LBL_END;
      i = 0;
      for (;;) {
         if (mp_cmp_d(&t1, 1uL) == MP_EQ) break;
         if ((err = mp_exptmod(&t1, &two, prime, &t1)) != MP_OKAY) goto LBL_END;
         i++;
      }
      if (i == 0u) {
         if ((err = mp_copy(&R, ret)) != MP_OKAY)                  goto LBL_END;
         err = MP_OKAY;
         goto LBL_END;
      }
      if ((err = mp_sub_d(&M, i, &t1)) != MP_OKAY)                goto LBL_END;
      if ((err = mp_sub_d(&t1, 1uL, &t1)) != MP_OKAY)             goto LBL_END;
      if ((err = mp_exptmod(&two, &t1, prime, &t1)) != MP_OKAY)   goto LBL_END;
      /* t1 = 2 ^ (M - i - 1) */
      if ((err = mp_exptmod(&C, &t1, prime, &t1)) != MP_OKAY)     goto LBL_END;
      /* t1 = C ^ (2 ^ (M - i - 1)) mod prime */
      if ((err = mp_sqrmod(&t1, prime, &C)) != MP_OKAY)           goto LBL_END;
      /* C = (t1 * t1) mod prime */
      if ((err = mp_mulmod(&R, &t1, prime, &R)) != MP_OKAY)       goto LBL_END;
      /* R = (R * t1) mod prime */
      if ((err = mp_mulmod(&T, &C, prime, &T)) != MP_OKAY)        goto LBL_END;
      /* T = (T * C) mod prime */
      mp_set(&M, i);
      /* M = i */
   }

I was able to fix this by adding the following after line 88

if (mp_cmp_d(&M, i) == MP_EQ) {err = MP_VAL; goto LBL_END;}

[czurnieden]

My use case only involves numbers 1 mod 4, so this may happen other times too

Unlikely, p \cong 1 mod 4 with p a perfect square of a prime A001248 is a special case for the Jacobi symbol. A001248. LibTomMath's primetest checks for perfect squares, so checking for primality would catch that.

I'm a bit torn here.
On one side: the Tonelli-Shanks algorithm used here assumes the modulus to be prime and the user is expected to check for that condition.
On the other side: LibTomMath should fail safely in all cases of error (if mathematically possible and economically feasible, of course!) especially in the case of wrong input.

The solution for composites p \cong 1 mod 4 for the first loop works but it can still run a long time if the factor is large.
The solution for composites p \cong 1 mod 4) for the first loop works but it can still run a long time if M is large (can be up to 2^INT_MAX large) but I have not yet checked if that is possible at all and if such a very large M is itself a red flag.

Still: both patches are cheap and allow the algorithm to halt in case of bad input.
(Can be made even cheaper because S and therefore M cannot get bigger than 2^INT_MAX)

What do my dear colleagues say? Is it a proper bug or not?

Sorry for the delay but the current situation is an unknown kind of stress; we were used to sales promising the moon to the customers, pointy-haired bosses insisting on implementing aforementioned sattelite and we won't even twitch anymore when we hear the supersonic boom of a passing deadline but this is a new one. A lot of us are also self-employed/have a little company and with some of our clients in financial troubles, way too often fatal nowadays, while others take advantage of the situation (to put it mildly!) we might not be as prompt in our reactions as the users of LibTomMath became accustomed to.

[friedrichsenm]

No problem with the delay!

FYI, the use case is for an algorithm that takes in a positive number and outputs its representation as the sum of four squares. High level, it makes a couple of random guesses for the first two square and then hopes that what's left is is a prime number 1 mod 4 so you can find its representation as the sum of two squares. You could do a primality test here before calling the function to make sure it will work, but adding that check might slow down the algorithm more than just trying to find the other two squares.