xref: /netbsd-src/external/lgpl3/gmp/dist/tests/mpz/t-perfsqr.c (revision 72c7faa4dbb41dbb0238d6b4a109da0d4b236dd4)
1 /* Test mpz_perfect_square_p.
2 
3 Copyright 2000-2002 Free Software Foundation, Inc.
4 
5 This file is part of the GNU MP Library test suite.
6 
7 The GNU MP Library test suite is free software; you can redistribute it
8 and/or modify it under the terms of the GNU General Public License as
9 published by the Free Software Foundation; either version 3 of the License,
10 or (at your option) any later version.
11 
12 The GNU MP Library test suite is distributed in the hope that it will be
13 useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General
15 Public License for more details.
16 
17 You should have received a copy of the GNU General Public License along with
18 the 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 
42 void
check_modulo(void)43 check_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. */
88 void
check_sqrt(int reps)89 check_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 
138 int
main(int argc,char ** argv)139 main (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