tests/mpz/t-perfsqr.c

   1 /* Test mpz_perfect_square_p.
   2
   3 Copyright 2000, 2001, 2002 Free Software Foundation, Inc.
   4
   5 This file is part of the GNU MP Library.
   6
   7 The GNU MP Library is free software; you can redistribute it and/or modify
   8 it under the terms of the GNU Lesser General Public License as published by
   9 the Free Software Foundation; either version 3 of the License, or (at your
  10 option) any later version.
  11
  12 The GNU MP Library is distributed in the hope that it will be useful, but
  13 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  14 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  15 License for more details.
  16
  17 You should have received a copy of the GNU Lesser General Public License
  18 along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  19
  20 #include <stdio.h>
  21 #include <stdlib.h>
  22
  23 #include "gmp.h"
  24 #include "gmp-impl.h"
  25 #include "tests.h"
  26
  27 #include "mpn/perfsqr.h"
  28
  29
  30 /* check_modulo() exercises mpz_perfect_square_p on squares which cover each
  31    possible quadratic residue to each divisor used within
  32    mpn_perfect_square_p, ensuring those residues aren't incorrectly claimed
  33    to be non-residues.
  34
  35    Each divisor is taken separately.  It's arranged that n is congruent to 0
  36    modulo the other divisors, 0 of course being a quadratic residue to any
  37    modulus.
  38
  39    The values "(j*others)^2" cover all quadratic residues mod divisor[i],
  40    but in no particular order.  j is run from 1<=j<=divisor[i] so that zero
  41    is excluded.  A literal n==0 doesn't reach the residue tests.  */
  42
  43 void
  44 check_modulo (void)
  45 {
  46   static const unsigned long  divisor[] = PERFSQR_DIVISORS;
  47   unsigned long  i, j;
  48
  49   mpz_t  alldiv, others, n;
  50
  51   mpz_init (alldiv);
  52   mpz_init (others);
  53   mpz_init (n);
  54
  55   /* product of all divisors */
  56   mpz_set_ui (alldiv, 1L);
  57   for (i = 0; i < numberof (divisor); i++)
  58     mpz_mul_ui (alldiv, alldiv, divisor[i]);
  59
  60   for (i = 0; i < numberof (divisor); i++)
  61     {
  62       /* product of all divisors except i */
  63       mpz_set_ui (others, 1L);
  64       for (j = 0; j < numberof (divisor); j++)
  65         if (i != j)
  66           mpz_mul_ui (others, others, divisor[j]);
  67
  68       for (j = 1; j <= divisor[i]; j++)
  69         {
  70           /* square */
  71           mpz_mul_ui (n, others, j);
  72           mpz_mul (n, n, n);
  73           if (! mpz_perfect_square_p (n))
  74             {
  75               printf ("mpz_perfect_square_p got 0, want 1\n");
  76               mpz_trace ("  n", n);
  77               abort ();
  78             }
  79         }
  80     }
  81
  82   mpz_clear (alldiv);
  83   mpz_clear (others);
  84   mpz_clear (n);
  85 }
  86
  87
  88 /* Exercise mpz_perfect_square_p compared to what mpz_sqrt says. */
  89 void
  90 check_sqrt (int reps)
  91 {
  92   mpz_t x2, x2t, x;
  93   mp_size_t x2n;
  94   int res;
  95   int i;
  96   /* int cnt = 0; */
  97   gmp_randstate_ptr rands = RANDS;
  98   mpz_t bs;
  99
 100   mpz_init (bs);
 101
 102   mpz_init (x2);
 103   mpz_init (x);
 104   mpz_init (x2t);
 105
 106   for (i = 0; i < reps; i++)
 107     {
 108       mpz_urandomb (bs, rands, 9);
 109       x2n = mpz_get_ui (bs);
 110       mpz_rrandomb (x2, rands, x2n);
 111       /* mpz_out_str (stdout, -16, x2); puts (""); */
 112
 113       res = mpz_perfect_square_p (x2);
 114       mpz_sqrt (x, x2);
 115       mpz_mul (x2t, x, x);
 116
 117       if (res != (mpz_cmp (x2, x2t) == 0))
 118         {
 119           printf    ("mpz_perfect_square_p and mpz_sqrt differ\n");
 120           mpz_trace ("   x  ", x);
 121           mpz_trace ("   x2 ", x2);
 122           mpz_trace ("   x2t", x2t);
 123           printf    ("   mpz_perfect_square_p %d\n", res);
 124           printf    ("   mpz_sqrt             %d\n", mpz_cmp (x2, x2t) == 0);
 125           abort ();
 126         }
 127
 128       /* cnt += res != 0; */
 129     }
 130   /* printf ("%d/%d perfect squares\n", cnt, reps); */
 131
 132   mpz_clear (bs);
 133   mpz_clear (x2);
 134   mpz_clear (x);
 135   mpz_clear (x2t);
 136 }
 137
 138
 139 int
 140 main (int argc, char **argv)
 141 {
 142   int reps = 200000;
 143
 144   tests_start ();
 145   mp_trace_base = -16;
 146
 147   if (argc == 2)
 148      reps = atoi (argv[1]);
 149
 150   check_modulo ();
 151   check_sqrt (reps);
 152
 153   tests_end ();
 154   exit (0);
 155 }