1 /* mpz_perfect_square_p(arg) -- Return non-zero if ARG is a pefect square,
4 Copyright (C) 1991 Free Software Foundation, Inc.
6 This file is part of the GNU MP Library.
8 The GNU MP Library is free software; you can redistribute it and/or modify
9 it under the terms of the GNU General Public License as published by
10 the Free Software Foundation; either version 2, or (at your option)
13 The GNU MP Library is distributed in the hope that it will be useful,
14 but WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 GNU General Public License for more details.
18 You should have received a copy of the GNU General Public License
19 along with the GNU MP Library; see the file COPYING. If not, write to
20 the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
26 #if BITS_PER_MP_LIMB == 32
27 static unsigned int primes[] = {3, 5, 7, 11, 13, 17, 19, 23, 29};
28 static unsigned long int residue_map[] =
29 {0x3, 0x13, 0x17, 0x23b, 0x161b, 0x1a317, 0x30af3, 0x5335f, 0x13d122f3};
31 #define PP 0xC0CFD797L /* 3 x 5 x 7 x 11 x 13 x ... x 29 */
34 /* sq_res_0x100[x mod 0x100] == 1 iff x mod 0x100 is a quadratic residue
36 static char sq_res_0x100[0x100] =
38 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,
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 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,
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 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,
43 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,
44 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,
45 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,
50 mpz_perfect_square_p (const MP_INT *a)
52 mpz_perfect_square_p (a)
58 mp_size asize = a->size;
59 mp_srcptr aptr = a->d;
63 /* No negative numbers are perfect squares. */
67 /* The first test excludes 55/64 (85.9%) of the perfect square candidates
69 if (sq_res_0x100[aptr[0] % 0x100] == 0)
72 #if BITS_PER_MP_LIMB == 32
73 /* The second test excludes 30652543/30808063 (99.5%) of the remaining
74 perfect square candidates in O(n) time. */
76 /* Firstly, compute REM = A mod PP. */
91 udiv_qrnnd (dummy, n1, n1, n0, PP);
95 /* We have A mod PP in REM. Now decide if REM is a quadratic residue
96 modulo the factors in PP. */
97 for (i = 0; i < (sizeof primes) / sizeof (int); i++)
103 if ((residue_map[i] & (1L << rem)) == 0)
108 /* For the third and last test, we finally compute the square root,
109 to make sure we've really got a perfect square. */
110 root_ptr = (mp_ptr) alloca ((asize + 1) / 2 * BYTES_PER_MP_LIMB);
112 /* Iff mpn_sqrt returns zero, the square is perfect. */
114 int retval = !mpn_sqrt (root_ptr, NULL, aptr, asize);