ghc/runtime/gmp/mpz_perfsqr.c

   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 }