1 /*
2
3 Copyright 2012, 2014, 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 <limits.h>
21 #include <stdlib.h>
22 #include <stdio.h>
23
24 #include "testutils.h"
25
26 #define MAXBITS 400
27 #define COUNT 9000
28
29 /* Called when s is supposed to be floor(sqrt(u)), and r = u - s^2 */
30 static int
31 sqrtrem_valid_p (const mpz_t u, const mpz_t s, const mpz_t r)
32 {
33 mpz_t t;
34
35 mpz_init (t);
36 mpz_mul (t, s, s);
37 mpz_sub (t, u, t);
38 if (mpz_sgn (t) < 0 || mpz_cmp (t, r) != 0)
39 {
40 mpz_clear (t);
41 return 0;
42 }
43 mpz_add_ui (t, s, 1);
44 mpz_mul (t, t, t);
45 if (mpz_cmp (t, u) <= 0)
46 {
47 mpz_clear (t);
48 return 0;
49 }
50
51 mpz_clear (t);
52 return 1;
53 }
54
55 void
56 mpz_mpn_sqrtrem (mpz_t s, mpz_t r, const mpz_t u)
57 {
58 mp_limb_t *sp, *rp;
59 mp_size_t un, sn, ret;
60
61 un = mpz_size (u);
62
63 mpz_xor (s, s, u);
64 sn = (un + 1) / 2;
65 sp = mpz_limbs_write (s, sn + 1);
66 sp [sn] = 11;
67
68 if (un & 1)
69 rp = NULL; /* Exploits the fact that r already is correct. */
70 else {
71 mpz_add (r, u, s);
72 rp = mpz_limbs_write (r, un + 1);
73 rp [un] = 19;
74 }
75
76 ret = mpn_sqrtrem (sp, rp, mpz_limbs_read (u), un);
77
78 if (sp [sn] != 11)
79 {
80 fprintf (stderr, "mpn_sqrtrem buffer overrun on sp.\n");
81 abort ();
82 }
83 if (un & 1) {
84 if ((ret != 0) != (mpz_size (r) != 0)) {
85 fprintf (stderr, "mpn_sqrtrem wrong return value with NULL.\n");
86 abort ();
87 }
88 } else {
89 mpz_limbs_finish (r, ret);
90 if ((size_t) ret != mpz_size (r)) {
91 fprintf (stderr, "mpn_sqrtrem wrong return value.\n");
92 abort ();
93 }
94 if (rp [un] != 19)
95 {
96 fprintf (stderr, "mpn_sqrtrem buffer overrun on rp.\n");
97 abort ();
98 }
99 }
100
101 mpz_limbs_finish (s, (un + 1) / 2);
102 }
103
104 void
105 testmain (int argc, char **argv)
106 {
107 unsigned i;
108 mpz_t u, s, r;
109
110 mpz_init (s);
111 mpz_init (r);
112
113 mpz_init_set_si (u, -1);
114 if (mpz_perfect_square_p (u))
115 {
116 fprintf (stderr, "mpz_perfect_square_p failed on -1.\n");
117 abort ();
118 }
119
120 if (!mpz_perfect_square_p (s))
121 {
122 fprintf (stderr, "mpz_perfect_square_p failed on 0.\n");
123 abort ();
124 }
125
126 for (i = 0; i < COUNT; i++)
127 {
128 mini_rrandomb (u, MAXBITS - (i & 0xFF));
129 mpz_sqrtrem (s, r, u);
130
131 if (!sqrtrem_valid_p (u, s, r))
132 {
133 fprintf (stderr, "mpz_sqrtrem failed:\n");
134 dump ("u", u);
135 dump ("sqrt", s);
136 dump ("rem", r);
137 abort ();
138 }
139
140 mpz_mpn_sqrtrem (s, r, u);
141
142 if (!sqrtrem_valid_p (u, s, r))
143 {
144 fprintf (stderr, "mpn_sqrtrem failed:\n");
145 dump ("u", u);
146 dump ("sqrt", s);
147 dump ("rem", r);
148 abort ();
149 }
150
151 if (mpz_sgn (r) == 0) {
152 mpz_neg (u, u);
153 mpz_sub_ui (u, u, 1);
154 }
155
156 if ((mpz_sgn (u) <= 0 || (i & 1)) ?
157 mpz_perfect_square_p (u) :
158 mpn_perfect_square_p (mpz_limbs_read (u), mpz_size (u)))
159 {
160 fprintf (stderr, "mp%s_perfect_square_p failed on non square:\n",
161 (mpz_sgn (u) <= 0 || (i & 1)) ? "z" : "n");
162 dump ("u", u);
163 abort ();
164 }
165
166 mpz_mul (u, s, s);
167 if (!((mpz_sgn (u) <= 0 || (i & 1)) ?
168 mpz_perfect_square_p (u) :
169 mpn_perfect_square_p (mpz_limbs_read (u), mpz_size (u))))
170 {
171 fprintf (stderr, "mp%s_perfect_square_p failed on square:\n",
172 (mpz_sgn (u) <= 0 || (i & 1)) ? "z" : "n");
173 dump ("u", u);
174 abort ();
175 }
176
177 }
178 mpz_clear (u);
179 mpz_clear (s);
180 mpz_clear (r);
181 }