libs/libtommath/bn_mp_is_square.c

   1 #include <tommath.h>
   2 #ifdef BN_MP_IS_SQUARE_C
   3 /* LibTomMath, multiple-precision integer library -- Tom St Denis
   4  *
   5  * LibTomMath is a library that provides multiple-precision
   6  * integer arithmetic as well as number theoretic functionality.
   7  *
   8  * The library was designed directly after the MPI library by
   9  * Michael Fromberger but has been written from scratch with
  10  * additional optimizations in place.
  11  *
  12  * The library is free for all purposes without any express
  13  * guarantee it works.
  14  *
  15  * Tom St Denis, tomstdenis@gmail.com, http://math.libtomcrypt.com
  16  */
  17
  18 /* Check if remainders are possible squares - fast exclude non-squares */
  19 static const char rem_128[128] = {
  20  0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
  21  0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
  22  1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
  23  1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
  24  0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
  25  1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
  26  1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
  27  1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1
  28 };
  29
  30 static const char rem_105[105] = {
  31  0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1,
  32  0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1,
  33  0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1,
  34  1, 0, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1,
  35  0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1,
  36  1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 0, 1, 1, 1, 1,
  37  1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1
  38 };
  39
  40 /* Store non-zero to ret if arg is square, and zero if not */
  41 int mp_is_square(mp_int *arg,int *ret)
  42 {
  43   int           res;
  44   mp_digit      c;
  45   mp_int        t;
  46   unsigned long r;
  47
  48   /* Default to Non-square :) */
  49   *ret = MP_NO;
  50
  51   if (arg->sign == MP_NEG) {
  52     return MP_VAL;
  53   }
  54
  55   /* digits used?  (TSD) */
  56   if (arg->used == 0) {
  57      return MP_OKAY;
  58   }
  59
  60   /* First check mod 128 (suppose that DIGIT_BIT is at least 7) */
  61   if (rem_128[127 & DIGIT(arg,0)] == 1) {
  62      return MP_OKAY;
  63   }
  64
  65   /* Next check mod 105 (3*5*7) */
  66   if ((res = mp_mod_d(arg,105,&c)) != MP_OKAY) {
  67      return res;
  68   }
  69   if (rem_105[c] == 1) {
  70      return MP_OKAY;
  71   }
  72
  73
  74   if ((res = mp_init_set_int(&t,11L*13L*17L*19L*23L*29L*31L)) != MP_OKAY) {
  75      return res;
  76   }
  77   if ((res = mp_mod(arg,&t,&t)) != MP_OKAY) {
  78      goto ERR;
  79   }
  80   r = mp_get_int(&t);
  81   /* Check for other prime modules, note it's not an ERROR but we must
  82    * free "t" so the easiest way is to goto ERR.  We know that res
  83    * is already equal to MP_OKAY from the mp_mod call
  84    */
  85   if ( (1L<<(r%11)) & 0x5C4L )             goto ERR;
  86   if ( (1L<<(r%13)) & 0x9E4L )             goto ERR;
  87   if ( (1L<<(r%17)) & 0x5CE8L )            goto ERR;
  88   if ( (1L<<(r%19)) & 0x4F50CL )           goto ERR;
  89   if ( (1L<<(r%23)) & 0x7ACCA0L )          goto ERR;
  90   if ( (1L<<(r%29)) & 0xC2EDD0CL )         goto ERR;
  91   if ( (1L<<(r%31)) & 0x6DE2B848L )        goto ERR;
  92
  93   /* Final check - is sqr(sqrt(arg)) == arg ? */
  94   if ((res = mp_sqrt(arg,&t)) != MP_OKAY) {
  95      goto ERR;
  96   }
  97   if ((res = mp_sqr(&t,&t)) != MP_OKAY) {
  98      goto ERR;
  99   }
 100
 101   *ret = (mp_cmp_mag(&t,arg) == MP_EQ) ? MP_YES : MP_NO;
 102 ERR:mp_clear(&t);
 103   return res;
 104 }
 105 #endif
 106
 107 /* $Source: /cvs/libtom/libtommath/bn_mp_is_square.c,v $ */
 108 /* $Revision: 1.3 $ */
 109 /* $Date: 2006/03/31 14:18:44 $ */