Annotation of OpenXM_contrib/gmp/mpn/generic/perfsqr.c, Revision 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: }
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