1/*
2
3Copyright 2012, 2014, 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 <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 */
30static int
31sqrtrem_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
55void
56mpz_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
104void
105testmain (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}
182