mp_is_square.c

   1 #include "tommath_private.h"
   2 #ifdef MP_IS_SQUARE_C
   3 /* LibTomMath, multiple-precision integer library -- Tom St Denis */
   4 /* SPDX-License-Identifier: Unlicense */
   5
   6 /* Check if remainders are possible squares - fast exclude non-squares */
   7 static const char rem_128[128] = {
   8    0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
   9    0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
  10    1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
  11    1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
  12    0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
  13    1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
  14    1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
  15    1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1
  16 };
  17
  18 static const char rem_105[105] = {
  19    0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1,
  20    0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1,
  21    0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1,
  22    1, 0, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1,
  23    0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1,
  24    1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 0, 1, 1, 1, 1,
  25    1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1
  26 };
  27
  28 /* Store non-zero to ret if arg is square, and zero if not */
  29 mp_err mp_is_square(const mp_int *arg, bool *ret)
  30 {
  31    mp_err   err;
  32    mp_digit c;
  33    mp_int   t;
  34    uint32_t r;
  35
  36    /* Default to Non-square :) */
  37    *ret = false;
  38
  39    if (mp_isneg(arg)) {
  40       return MP_VAL;
  41    }
  42
  43    if (mp_iszero(arg)) {
  44       *ret = true;
  45       return MP_OKAY;
  46    }
  47
  48    /* First check mod 128 (suppose that MP_DIGIT_BIT is at least 7) */
  49    if (rem_128[127u & arg->dp[0]] == (char)1) {
  50       return MP_OKAY;
  51    }
  52
  53    /* Next check mod 105 (3*5*7) */
  54    if ((err = mp_mod_d(arg, 105uL, &c)) != MP_OKAY) {
  55       return err;
  56    }
  57    if (rem_105[c] == (char)1) {
  58       return MP_OKAY;
  59    }
  60
  61
  62    if ((err = mp_init_u32(&t, 11u*13u*17u*19u*23u*29u*31u)) != MP_OKAY) {
  63       return err;
  64    }
  65    if ((err = mp_mod(arg, &t, &t)) != MP_OKAY) {
  66       goto LBL_ERR;
  67    }
  68    r = mp_get_u32(&t);
  69    /* Check for other prime modules, note it's not an ERROR but we must
  70     * free "t" so the easiest way is to goto LBL_ERR.  We know that err
  71     * is already equal to MP_OKAY from the mp_mod call
  72     */
  73    if (((1uL<<(r%11uL)) & 0x5C4uL) != 0uL)         goto LBL_ERR;
  74    if (((1uL<<(r%13uL)) & 0x9E4uL) != 0uL)         goto LBL_ERR;
  75    if (((1uL<<(r%17uL)) & 0x5CE8uL) != 0uL)        goto LBL_ERR;
  76    if (((1uL<<(r%19uL)) & 0x4F50CuL) != 0uL)       goto LBL_ERR;
  77    if (((1uL<<(r%23uL)) & 0x7ACCA0uL) != 0uL)      goto LBL_ERR;
  78    if (((1uL<<(r%29uL)) & 0xC2EDD0CuL) != 0uL)     goto LBL_ERR;
  79    if (((1uL<<(r%31uL)) & 0x6DE2B848uL) != 0uL)    goto LBL_ERR;
  80
  81    /* Final check - is sqr(sqrt(arg)) == arg ? */
  82    if ((err = mp_sqrt(arg, &t)) != MP_OKAY) {
  83       goto LBL_ERR;
  84    }
  85    if ((err = mp_sqr(&t, &t)) != MP_OKAY) {
  86       goto LBL_ERR;
  87    }
  88
  89    *ret = (mp_cmp_mag(&t, arg) == MP_EQ);
  90 LBL_ERR:
  91    mp_clear(&t);
  92    return err;
  93 }
  94 #endif