Return to perfsqr.c CVS log

Up to [local] / OpenXM_contrib / gmp / mpn / generic

Annotation of OpenXM_contrib/gmp/mpn/generic/perfsqr.c, Revision 1.1.1.2

1.1       maekawa     1: /* mpn_perfect_square_p(u,usize) -- Return non-zero if U is a perfect square,
                      2:    zero otherwise.
                      3:
1.1.1.2 ! maekawa     4: Copyright (C) 1991, 1993, 1994, 1996, 1997, 2000 Free Software Foundation,
        !             5: Inc.
1.1       maekawa     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
1.1.1.2 ! maekawa    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
1.1       maekawa    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
1.1.1.2 ! maekawa    16: or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
1.1       maekawa    17: License for more details.
                     18:
1.1.1.2 ! maekawa    19: You should have received a copy of the GNU Lesser General Public License
1.1       maekawa    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:
1.1.1.2 ! maekawa    24: #include <stdio.h> /* for NULL */
1.1       maekawa    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
1.1.1.2 ! maekawa    80:   if (((CNST_LIMB(0x12DD703303AED3) >> rem % 53) & 1) == 0)
1.1       maekawa    81:     return 0;
1.1.1.2 ! maekawa    82:   if (((CNST_LIMB(0x4351B2753DF) >> rem % 47) & 1) == 0)
1.1       maekawa    83:     return 0;
1.1.1.2 ! maekawa    84:   if (((CNST_LIMB(0x35883A3EE53) >> rem % 43) & 1) == 0)
1.1       maekawa    85:     return 0;
1.1.1.2 ! maekawa    86:   if (((CNST_LIMB(0x1B382B50737) >> rem % 41) & 1) == 0)
1.1       maekawa    87:     return 0;
1.1.1.2 ! maekawa    88:   if (((CNST_LIMB(0x165E211E9B) >> rem % 37) & 1) == 0)
1.1       maekawa    89:     return 0;
1.1.1.2 ! maekawa    90:   if (((CNST_LIMB(0x121D47B7) >> rem % 31) & 1) == 0)
1.1       maekawa    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: }