rts/gmp/mpn/generic/perfsqr.c

   1 /* mpn_perfect_square_p(u,usize) -- Return non-zero if U is a perfect square,
   2    zero otherwise.
   3
   4 Copyright (C) 1991, 1993, 1994, 1996, 1997, 2000 Free Software Foundation,
   5 Inc.
   6
   7 This file is part of the GNU MP Library.
   8
   9 The GNU MP Library is free software; you can redistribute it and/or modify
  10 it under the terms of the GNU Lesser General Public License as published by
  11 the Free Software Foundation; either version 2.1 of the License, or (at your
  12 option) any later version.
  13
  14 The GNU MP Library is distributed in the hope that it will be useful, but
  15 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  16 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  17 License for more details.
  18
  19 You should have received a copy of the GNU Lesser General Public License
  20 along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
  21 the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
  22 MA 02111-1307, USA. */
  23
  24 #include <stdio.h> /* for NULL */
  25 #include "gmp.h"
  26 #include "gmp-impl.h"
  27 #include "longlong.h"
  28
  29
  30 /* sq_res_0x100[x mod 0x100] == 1 iff x mod 0x100 is a quadratic residue
  31    modulo 0x100.  */
  32 static unsigned char const sq_res_0x100[0x100] =
  33 {
  34   1,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
  35   0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
  36   1,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
  37   0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
  38   0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
  39   0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
  40   0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
  41   0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
  42 };
  43
  44 int
  45 #if __STDC__
  46 mpn_perfect_square_p (mp_srcptr up, mp_size_t usize)
  47 #else
  48 mpn_perfect_square_p (up, usize)
  49      mp_srcptr up;
  50      mp_size_t usize;
  51 #endif
  52 {
  53   mp_limb_t rem;
  54   mp_ptr root_ptr;
  55   int res;
  56   TMP_DECL (marker);
  57
  58   /* The first test excludes 55/64 (85.9%) of the perfect square candidates
  59      in O(1) time.  */
  60   if ((sq_res_0x100[(unsigned int) up[0] % 0x100] & 1) == 0)
  61     return 0;
  62
  63 #if defined (PP)
  64   /* The second test excludes 30652543/30808063 (99.5%) of the remaining
  65      perfect square candidates in O(n) time.  */
  66
  67   /* Firstly, compute REM = A mod PP.  */
  68   if (UDIV_TIME > (2 * UMUL_TIME + 6))
  69     rem = mpn_preinv_mod_1 (up, usize, (mp_limb_t) PP, (mp_limb_t) PP_INVERTED);
  70   else
  71     rem = mpn_mod_1 (up, usize, (mp_limb_t) PP);
  72
  73   /* Now decide if REM is a quadratic residue modulo the factors in PP.  */
  74
  75   /* If A is just a few limbs, computing the square root does not take long
  76      time, so things might run faster if we limit this loop according to the
  77      size of A.  */
  78
  79 #if BITS_PER_MP_LIMB == 64
  80   if (((CNST_LIMB(0x12DD703303AED3) >> rem % 53) & 1) == 0)
  81     return 0;
  82   if (((CNST_LIMB(0x4351B2753DF) >> rem % 47) & 1) == 0)
  83     return 0;
  84   if (((CNST_LIMB(0x35883A3EE53) >> rem % 43) & 1) == 0)
  85     return 0;
  86   if (((CNST_LIMB(0x1B382B50737) >> rem % 41) & 1) == 0)
  87     return 0;
  88   if (((CNST_LIMB(0x165E211E9B) >> rem % 37) & 1) == 0)
  89     return 0;
  90   if (((CNST_LIMB(0x121D47B7) >> rem % 31) & 1) == 0)
  91     return 0;
  92 #endif
  93   if (((0x13D122F3L >> rem % 29) & 1) == 0)
  94     return 0;
  95   if (((0x5335FL >> rem % 23) & 1) == 0)
  96     return 0;
  97   if (((0x30AF3L >> rem % 19) & 1) == 0)
  98     return 0;
  99   if (((0x1A317L >> rem % 17) & 1) == 0)
 100     return 0;
 101   if (((0x161BL >> rem % 13) & 1) == 0)
 102     return 0;
 103   if (((0x23BL >> rem % 11) & 1) == 0)
 104     return 0;
 105   if (((0x017L >> rem % 7) & 1) == 0)
 106     return 0;
 107   if (((0x13L >> rem % 5) & 1) == 0)
 108     return 0;
 109   if (((0x3L >> rem % 3) & 1) == 0)
 110     return 0;
 111 #endif
 112
 113   TMP_MARK (marker);
 114
 115   /* For the third and last test, we finally compute the square root,
 116      to make sure we've really got a perfect square.  */
 117   root_ptr = (mp_ptr) TMP_ALLOC ((usize + 1) / 2 * BYTES_PER_MP_LIMB);
 118
 119   /* Iff mpn_sqrtrem returns zero, the square is perfect.  */
 120   res = ! mpn_sqrtrem (root_ptr, NULL, up, usize);
 121   TMP_FREE (marker);
 122   return res;
 123 }