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.1

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