Annotation of OpenXM_contrib/gmp/mpn/generic/perfsqr.c, Revision 1.1.1.3
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.3 ! ohara 4: Copyright 1991, 1993, 1994, 1996, 1997, 2000, 2001 Free Software Foundation,
1.1.1.2 maekawa 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,
1.1.1.3 ! ohara 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
1.1 maekawa 42: };
43:
44: int
45: mpn_perfect_square_p (mp_srcptr up, mp_size_t usize)
46: {
47: mp_limb_t rem;
48: mp_ptr root_ptr;
49: int res;
50: TMP_DECL (marker);
51:
1.1.1.3 ! ohara 52: ASSERT (usize >= 1);
! 53:
! 54: /* The first test excludes 212/256 (82.8%) of the perfect square candidates
1.1 maekawa 55: in O(1) time. */
1.1.1.3 ! ohara 56: if ((sq_res_0x100[(unsigned int) up[0] % 0x100]) == 0)
1.1 maekawa 57: return 0;
58:
1.1.1.3 ! ohara 59: #if 0
! 60: /* Check that we have even multiplicity of 2, and then check that the rest is
! 61: a possible perfect square. Leave disabled until we can determine this
! 62: really is an improvement. It it is, it could completely replace the
! 63: simple probe above, since this should through out more non-squares, but at
! 64: the expense of somewhat more cycles. */
! 65: {
! 66: mp_limb_t lo;
! 67: int cnt;
! 68: lo = up[0];
! 69: while (lo == 0)
! 70: up++, lo = up[0], usize--;
! 71: count_trailing_zeros (cnt, lo);
! 72: if ((cnt & 1) != 0)
! 73: return 0; /* return of not even multiplicity of 2 */
! 74: lo >>= cnt; /* shift down to align lowest non-zero bit */
! 75: lo >>= 1; /* shift away lowest non-zero bit */
! 76: if ((lo & 3) != 0)
! 77: return 0;
! 78: }
! 79: #endif
! 80:
1.1 maekawa 81: #if defined (PP)
1.1.1.3 ! ohara 82: /* The second test excludes most of the remaining perfect square candidates
! 83: in O(n) time. Each residue tested keeps (p+1)/2 out of p, and with no
! 84: common factors the tests are independent so the total passed is
! 85:
! 86: p=M (p+1)/2
! 87: product -------
! 88: p=3 p
! 89:
! 90: which, for M=29 (for 32-bit hosts) is 155520/30808063, or 99.5% excluded,
! 91: and for M=53 (for 64-bit hosts) is 67722117120/742518990138757, or 99.99%
! 92: excluded.
! 93:
! 94: If U is just a few limbs, computing the square root does not take long
! 95: time, so things might run faster if we suppress this for A less than some
! 96: appropriate threshold. */
! 97:
! 98:
! 99: /* Firstly, compute REM = U mod PP. */
! 100: #if defined (PP_INVERTED)
! 101: rem = MPN_MOD_OR_PREINV_MOD_1 (up, usize, (mp_limb_t) PP,
! 102: (mp_limb_t) PP_INVERTED);
! 103: #else
! 104: rem = mpn_mod_1 (up, usize, (mp_limb_t) PP);
! 105: #endif
1.1 maekawa 106:
1.1.1.3 ! ohara 107: /* Now decide if REM is a quadratic residue modulo the factors in PP. We
! 108: test prime pairs multiplied together, as long as their products don't
! 109: overflow a limb. The tests that exclude most false squares are put
! 110: first. */
1.1 maekawa 111:
1.1.1.3 ! ohara 112: if (0
1.1 maekawa 113: #if BITS_PER_MP_LIMB == 64
1.1.1.3 ! ohara 114: || ((CNST_LIMB(0x0C22C90530902D3) >> rem % (3 * 19)) & 1) == 0
! 115: || ((CNST_LIMB(0x0230148611CA33) >> rem % (5 * 11)) & 1) == 0
! 116: || ((CNST_LIMB(0x12DD703303AED3) >> rem % 53) & 1) == 0
! 117: || ((CNST_LIMB(0x04351B2753DF) >> rem % 47) & 1) == 0
! 118: || ((CNST_LIMB(0x35883A3EE53) >> rem % 43) & 1) == 0
! 119: || ((CNST_LIMB(0x1B382B50737) >> rem % 41) & 1) == 0
! 120: || ((CNST_LIMB(0x165E211E9B) >> rem % 37) & 1) == 0
! 121: || ((0x121D47B7L >> rem % 31) & 1) == 0
! 122: || ((0x13D122F3L >> rem % 29) & 1) == 0
! 123: || ((0x05335FL >> rem % 23) & 1) == 0
! 124: /* 19 handled above */
! 125: || ((0x1A317L >> rem % 17) & 1) == 0
! 126: || ((0x161BL >> rem % 13) & 1) == 0
! 127: /* 11 handled above */
! 128: || ((0x17L >> rem % 7) & 1) == 0
! 129: /* 5 handled above */
! 130: /* 3 handled above */
1.1 maekawa 131: #endif
1.1.1.3 ! ohara 132: #if BITS_PER_MP_LIMB == 32
! 133: || ((0x058293L >> rem % (3 * 7)) & 1) == 0
! 134: || ((0x13D122F3L >> rem % 29) & 1) == 0
! 135: || ((0x05335FL >> rem % 23) & 1) == 0
! 136: || ((0x30AF3L >> rem % 19) & 1) == 0
! 137: || ((0x1A317L >> rem % 17) & 1) == 0
! 138: || ((0x161BL >> rem % 13) & 1) == 0
! 139: || ((0x23BL >> rem % 11) & 1) == 0
! 140: /* 7 handled above */
! 141: || ((0x13L >> rem % 5) & 1) == 0
! 142: /* 3 handled above */
! 143: #endif
! 144: #if BITS_PER_MP_LIMB == 16
! 145: || ((0x0653L >> rem % (3 * 5)) & 1) == 0
! 146: || ((0x161BL >> rem % 13) & 1) == 0
! 147: || ((0x23BL >> rem % 11) & 1) == 0
! 148: || ((0x17L >> rem % 7) & 1) == 0
! 149: /* 5 handled above */
! 150: /* 3 handled above */
! 151: #endif
! 152: #if BITS_PER_MP_LIMB == 8
! 153: || ((0x17L >> rem % 7) & 1) == 0
! 154: || ((0x13L >> rem % 5) & 1) == 0
! 155: || ((0x3L >> rem % 3) & 1) == 0
1.1 maekawa 156: #endif
1.1.1.3 ! ohara 157: #if BITS_PER_MP_LIMB == 4
! 158: || ((0x13L >> rem % 5) & 1) == 0
! 159: || ((0x3L >> rem % 3) & 1) == 0
! 160: #endif
! 161: #if BITS_PER_MP_LIMB == 2
! 162: || ((0x3L >> rem % 3) & 1) == 0
! 163: #endif
! 164: )
! 165: return 0;
! 166: #endif /* PP */
1.1 maekawa 167:
168: TMP_MARK (marker);
169:
170: /* For the third and last test, we finally compute the square root,
171: to make sure we've really got a perfect square. */
172: root_ptr = (mp_ptr) TMP_ALLOC ((usize + 1) / 2 * BYTES_PER_MP_LIMB);
173:
174: /* Iff mpn_sqrtrem returns zero, the square is perfect. */
175: res = ! mpn_sqrtrem (root_ptr, NULL, up, usize);
176: TMP_FREE (marker);
177: return res;
178: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to perfsqr.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpn / generic