libtommath/bn_mp_is_square.c

#include <tommath_private.h>
#ifdef BN_MP_IS_SQUARE_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 was designed directly after the MPI library by
 * Michael Fromberger but has been written from scratch with
 * additional optimizations in place.
 *
 * The library is free for all purposes without any express
 * guarantee it works.
 *
 * Tom St Denis, tstdenis82@gmail.com, http://libtom.org
 */

/* Check if remainders are possible squares - fast exclude non-squares */
static const char rem_128[128] = {
   0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
   0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
   1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
   1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
   0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
   1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
   1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
   1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1
};

static const char rem_105[105] = {
   0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1,
   0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1,
   0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1,
   1, 0, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1,
   0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1,
   1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 0, 1, 1, 1, 1,
   1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1
};

/* Store non-zero to ret if arg is square, and zero if not */
int mp_is_square(const mp_int *arg, int *ret)
{
   int           res;
   mp_digit      c;
   mp_int        t;
   unsigned long r;

   /* Default to Non-square :) */
   *ret = MP_NO;

   if (arg->sign == MP_NEG) {
      return MP_VAL;
   }

   /* digits used?  (TSD) */
   if (arg->used == 0) {
      return MP_OKAY;
   }

   /* First check mod 128 (suppose that DIGIT_BIT is at least 7) */
   if (rem_128[127u & DIGIT(arg, 0)] == (char)1) {
      return MP_OKAY;
   }

   /* Next check mod 105 (3*5*7) */
   if ((res = mp_mod_d(arg, 105uL, &c)) != MP_OKAY) {
      return res;
   }
   if (rem_105[c] == (char)1) {
      return MP_OKAY;
   }


   if ((res = mp_init_set_int(&t, 11L*13L*17L*19L*23L*29L*31L)) != MP_OKAY) {
      return res;
   }
   if ((res = mp_mod(arg, &t, &t)) != MP_OKAY) {
      goto ERR;
   }
   r = mp_get_int(&t);
   /* Check for other prime modules, note it's not an ERROR but we must
    * free "t" so the easiest way is to goto ERR.  We know that res
    * is already equal to MP_OKAY from the mp_mod call
    */
   if (((1uL<<(r%11uL)) & 0x5C4uL) != 0uL)         goto ERR;
   if (((1uL<<(r%13uL)) & 0x9E4uL) != 0uL)         goto ERR;
   if (((1uL<<(r%17uL)) & 0x5CE8uL) != 0uL)        goto ERR;
   if (((1uL<<(r%19uL)) & 0x4F50CuL) != 0uL)       goto ERR;
   if (((1uL<<(r%23uL)) & 0x7ACCA0uL) != 0uL)      goto ERR;
   if (((1uL<<(r%29uL)) & 0xC2EDD0CuL) != 0uL)     goto ERR;
   if (((1uL<<(r%31uL)) & 0x6DE2B848uL) != 0uL)    goto ERR;

   /* Final check - is sqr(sqrt(arg)) == arg ? */
   if ((res = mp_sqrt(arg, &t)) != MP_OKAY) {
      goto ERR;
   }
   if ((res = mp_sqr(&t, &t)) != MP_OKAY) {
      goto ERR;
   }

   *ret = (mp_cmp_mag(&t, arg) == MP_EQ) ? MP_YES : MP_NO;
ERR:
   mp_clear(&t);
   return res;
}
#endif

/* ref:         $Format:%D$ */
/* git commit:  $Format:%H$ */
/* commit time: $Format:%ai$ */