home *** CD-ROM | disk | FTP | other *** search
/ InfoMagic Source Code 1993 July / THE_SOURCE_CODE_CD_ROM.iso / gnu / gmp-1.3.2 / mpz_perfsqr.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-05-18  |  3.4 KB  |  119 lines
  1. /* mpz_perfect_square_p(arg) -- Return non-zero if ARG is a pefect square,
  2.    zero otherwise.
  3.  
  4. Copyright (C) 1991 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 General Public License as published by
  10. the Free Software Foundation; either version 2, or (at your option)
  11. any later version.
  12.  
  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.
  17.  
  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.  */
  21.  
  22. #include "gmp.h"
  23. #include "gmp-impl.h"
  24. #include "longlong.h"
  25.  
  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};
  30.  
  31. #define PP 0xC0CFD797L        /* 3 x 5 x 7 x 11 x 13 x ... x 29 */
  32. #endif
  33.  
  34. /* sq_res_0x100[x mod 0x100] == 1 iff x mod 0x100 is a quadratic residue
  35.    modulo 0x100.  */
  36. static char sq_res_0x100[0x100] =
  37. {
  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,
  46. };
  47.  
  48. int
  49. #ifdef __STDC__
  50. mpz_perfect_square_p (const MP_INT *a)
  51. #else
  52. mpz_perfect_square_p (a)
  53.      const MP_INT *a;
  54. #endif
  55. {
  56.   mp_limb n1, n0;
  57.   mp_size i;
  58.   mp_size asize = a->size;
  59.   mp_srcptr aptr = a->d;
  60.   mp_limb rem;
  61.   mp_ptr root_ptr;
  62.  
  63.   /* No negative numbers are perfect squares.  */
  64.   if (asize < 0)
  65.     return 0;
  66.  
  67.   /* The first test excludes 55/64 (85.9%) of the perfect square candidates
  68.      in O(1) time.  */
  69.   if (sq_res_0x100[aptr[0] % 0x100] == 0)
  70.     return 0;
  71.  
  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.  */
  75.  
  76.   /* Firstly, compute REM = A mod PP.  */
  77.   n1 = aptr[asize - 1];
  78.   if (n1 >= PP)
  79.     {
  80.       n1 = 0;
  81.       i = asize - 1;
  82.     }
  83.   else
  84.     i = asize - 2;
  85.  
  86.   for (; i >= 0; i--)
  87.     {
  88.       mp_limb dummy;
  89.  
  90.       n0 = aptr[i];
  91.       udiv_qrnnd (dummy, n1, n1, n0, PP);
  92.     }
  93.   rem = n1;
  94.  
  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++)
  98.     {
  99.       unsigned int p;
  100.  
  101.       p = primes[i];
  102.       rem %= p;
  103.       if ((residue_map[i] & (1L << rem)) == 0)
  104.     return 0;
  105.     }
  106. #endif
  107.  
  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);
  111.  
  112.   /* Iff mpn_sqrt returns zero, the square is perfect.  */
  113.   {
  114.     int retval = !mpn_sqrt (root_ptr, NULL, aptr, asize);
  115.     alloca (0);
  116.     return retval;
  117.   }
  118. }
  119.