#include <tommath.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, [email protected], http://math.libtomcrypt.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(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[127 & DIGIT(arg,0)] == 1) {
     return MP_OKAY;
  }

  /* Next check mod 105 (3*5*7) */
  if ((res = mp_mod_d(arg,105,&c)) != MP_OKAY) {
     return res;
  }
  if (rem_105[c] == 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 ( (1L<<(r%11)) & 0x5C4L )             goto ERR;
  if ( (1L<<(r%13)) & 0x9E4L )             goto ERR;
  if ( (1L<<(r%17)) & 0x5CE8L )            goto ERR;
  if ( (1L<<(r%19)) & 0x4F50CL )           goto ERR;
  if ( (1L<<(r%23)) & 0x7ACCA0L )          goto ERR;
  if ( (1L<<(r%29)) & 0xC2EDD0CL )         goto ERR;
  if ( (1L<<(r%31)) & 0x6DE2B848L )        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