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