1/*
2Copyright 2017 Free Software Foundation, Inc.
3
4This file is part of the GNU MP Library test suite.
5
6The GNU MP Library test suite is free software; you can redistribute it
7and/or modify it under the terms of the GNU General Public License as
8published by the Free Software Foundation; either version 3 of the License,
9or (at your option) any later version.
10
11The GNU MP Library test suite is distributed in the hope that it will be
12useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
13MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General
14Public License for more details.
15
16You should have received a copy of the GNU General Public License along with
17the GNU MP Library test suite.  If not, see https://www.gnu.org/licenses/.  */
18
19/* Usage:
20
21   ./sqrtrem_1_2 x
22
23     Checks mpn_sqrtrem() exhaustively, starting from 0, incrementing
24     the operand by a single unit, until all values handled by
25     mpn_sqrtrem{1,2} are tested. SLOW.
26
27   ./sqrtrem_1_2 s 1
28
29     Checks some special cases for mpn_sqrtrem(). I.e. values of the form
30     2^k*i and 2^k*(i+1)-1, with k=2^n and 0<i<2^k, until all such values,
31     handled by mpn_sqrtrem{1,2}, are tested.
32     Currently supports only the test of values that fits in one limb.
33     Less slow than the exhaustive test.
34
35   ./sqrtrem_1_2 c
36
37     Checks all corner cases for mpn_sqrtrem(). I.e. values of the form
38     i*i and (i+1)*(i+1)-1, for each value of i, until all such values,
39     handled by mpn_sqrtrem{1,2}, are tested.
40     Slightly faster than the special cases test.
41
42   For larger values, use
43   ./try mpn_sqrtrem
44
45 */
46
47#include <stdlib.h>
48#include <stdio.h>
49#include "gmp-impl.h"
50#include "longlong.h"
51#include "tests.h"
52#define STOP(x) return (x)
53/* #define STOP(x) x */
54#define SPINNER(v)					\
55  do {							\
56    MPN_SIZEINBASE_2EXP (spinner_count, q, v, 1);	\
57    --spinner_count;					\
58    spinner();						\
59  } while (0)
60
61int something_wrong (mp_limb_t er, mp_limb_t ec, mp_limb_t es)
62{
63  fprintf (stderr, "root = %lu , rem = {%lu , %lu}\n", (long unsigned) es,(long unsigned) ec,(long unsigned) er);
64  return -1;
65}
66
67int
68check_all_values (int justone, int quick)
69{
70  mp_limb_t es, mer, er, s[1], r[2], q[2];
71  mp_size_t x;
72  unsigned bits;
73
74  es=1;
75  if (quick) {
76    printf ("Quick, skipping some... (%u)\n", GMP_NUMB_BITS - 2);
77    es <<= GMP_NUMB_BITS / 2 - 1;
78  }
79  er=0;
80  mer= es << 1;
81  *q = es * es;
82  printf ("All values tested, up to bits:\n");
83  do {
84    x = mpn_sqrtrem (s, r, q, 1);
85    if (UNLIKELY (x != (er != 0)) || UNLIKELY (*s != es)
86	|| UNLIKELY ((x == 1) && (er != *r)))
87      STOP (something_wrong (er, 0, es));
88
89    if (UNLIKELY (er == mer)) {
90      ++es;
91      if (UNLIKELY ((es & 0xff) == 0))
92	SPINNER(1);
93      mer +=2; /* mer = es * 2 */
94      er = 0;
95    } else
96      ++er;
97    ++*q;
98  } while (*q != 0);
99  q[1] = 1;
100  SPINNER(2);
101  printf ("\nValues of a single limb, tested.\n");
102  if (justone) return 0;
103  printf ("All values tested, up to bits:\n");
104  do {
105    x = mpn_sqrtrem (s, r, q, 2);
106    if (UNLIKELY (x != (er != 0)) || UNLIKELY (*s != es)
107	|| UNLIKELY ((x == 1) && (er != *r)))
108      STOP (something_wrong (er, 0, es));
109
110    if (UNLIKELY (er == mer)) {
111      ++es;
112      if (UNLIKELY ((es & 0x7f) == 0))
113	SPINNER(2);
114      mer +=2; /* mer = es * 2 */
115      if (UNLIKELY (mer == 0))
116	break;
117      er = 0;
118    } else
119      ++er;
120    q[1] += (++*q == 0);
121  } while (1);
122  SPINNER(2);
123  printf ("\nValues with at most a limb for reminder, tested.\n");
124  printf ("Testing more values not supported, jet.\n");
125  return 0;
126}
127
128mp_limb_t
129upd (mp_limb_t *s, mp_limb_t k)
130{
131  mp_limb_t _s = *s;
132
133  while (k > _s * 2)
134    {
135      k -= _s * 2 + 1;
136      ++_s;
137    }
138  *s = _s;
139  return k;
140}
141
142mp_limb_t
143upd1 (mp_limb_t *s, mp_limb_t k)
144{
145  mp_limb_t _s = *s;
146
147  if (LIKELY (k < _s * 2)) return k + 1;
148  *s = _s + 1;
149  return k - _s * 2;
150}
151
152int
153check_some_values (int justone, int quick)
154{
155  mp_limb_t es, her, er, k, s[1], r[2], q[2];
156  mp_size_t x;
157  unsigned bits;
158
159  es = 1 << 1;
160  if (quick) {
161    es <<= GMP_NUMB_BITS / 4 - 1;
162    printf ("Quick, skipping some... (%u)\n", GMP_NUMB_BITS / 2);
163  }
164  er = 0;
165  *q = es * es;
166  printf ("High-half values tested, up to bits:\n");
167  do {
168    k  = *q - 1;
169    do {
170      x = mpn_sqrtrem (s, r, q, 1);
171      if (UNLIKELY (x != (er != 0)) || UNLIKELY (*s != es)
172	  || UNLIKELY ((x == 1) && (er != *r)))
173	STOP (something_wrong (er, 0, es));
174
175      if (UNLIKELY ((es & 0xffff) == 0))
176	SPINNER(1);
177      if ((*q & k) == 0) {
178	*q |= k;
179	er = upd (&es, k + er);
180      } else {
181	++*q;
182	er = upd1 (&es, er);
183      }
184    } while (es & k);
185  } while (*q != 0);
186  q[1] = 1;
187  SPINNER(2);
188  printf ("\nValues of a single limb, tested.\n");
189  if (justone) return 0;
190  if (quick) {
191    es <<= GMP_NUMB_BITS / 2 - 1;
192    q[1] <<= GMP_NUMB_BITS - 2;
193    printf ("Quick, skipping some... (%u)\n", GMP_NUMB_BITS - 2);
194  }
195  printf ("High-half values tested, up to bits:\n");
196  do {
197    x = mpn_sqrtrem (s, r, q, 2);
198    if (UNLIKELY (x != (er != 0)) || UNLIKELY (*s != es)
199	|| UNLIKELY ((x == 1) && (er != *r)))
200      STOP (something_wrong (er, 0, es));
201
202    if (*q == 0) {
203      *q = GMP_NUMB_MAX;
204      if (UNLIKELY ((es & 0xffff) == 0)) {
205	if (UNLIKELY (es == GMP_NUMB_HIGHBIT))
206	  break;
207	SPINNER(2);
208      }
209      /* er = er + GMP_NUMB_MAX - 1 - es*2 // postponed */
210      ++es;
211      /* er = er + GMP_NUMB_MAX - 1 - 2*(es-1) =
212            = er +(GMP_NUMB_MAX + 1)- 2* es = er - 2*es */
213      er = upd (&es, er - 2 * es);
214    } else {
215      *q = 0;
216      ++q[1];
217      er = upd1 (&es, er);
218    }
219  } while (1);
220  SPINNER(2);
221  printf ("\nValues with at most a limb for reminder, tested.\n");
222  er = GMP_NUMB_MAX; her = 0;
223
224  printf ("High-half values tested, up to bits:\n");
225  do {
226    x = mpn_sqrtrem (s, r, q, 2);
227    if (UNLIKELY (x != (her?2:(er != 0))) || UNLIKELY (*s != es)
228	|| UNLIKELY ((x != 0) && ((er != *r) || ((x == 2) && (r[1] != 1)))))
229      STOP (something_wrong (er, her, es));
230
231    if (*q == 0) {
232      *q = GMP_NUMB_MAX;
233      if (UNLIKELY ((es & 0xffff) == 0)) {
234	SPINNER(2);
235      }
236      if (her) {
237	++es;
238	her = 0;
239	er = er - 2 * es;
240      } else {
241	her = --er != GMP_NUMB_MAX;
242	if (her & (er > es * 2)) {
243	  er -= es * 2 + 1;
244	  her = 0;
245	  ++es;
246	}
247      }
248    } else {
249      *q = 0;
250      if (++q[1] == 0) break;
251      if ((her == 0) | (er < es * 2)) {
252	her += ++er == 0;
253      }	else {
254	  er -= es * 2;
255	  her = 0;
256	  ++es;
257      }
258    }
259  } while (1);
260  printf ("| %u\nValues of at most two limbs, tested.\n", GMP_NUMB_BITS*2);
261  return 0;
262}
263
264int
265check_corner_cases (int justone, int quick)
266{
267  mp_limb_t es, er, s[1], r[2], q[2];
268  mp_size_t x;
269  unsigned bits;
270
271  es = 1;
272  if (quick) {
273    es <<= GMP_NUMB_BITS / 2 - 1;
274    printf ("Quick, skipping some... (%u)\n", GMP_NUMB_BITS - 2);
275  }
276  er = 0;
277  *q = es*es;
278  printf ("Corner cases tested, up to bits:\n");
279  do {
280    x = mpn_sqrtrem (s, r, q, 1);
281    if (UNLIKELY (x != (er != 0)) || UNLIKELY (*s != es)
282	|| UNLIKELY ((x == 1) && (er != *r)))
283      STOP (something_wrong (er, 0, es));
284
285    if (er != 0) {
286      ++es;
287      if (UNLIKELY ((es & 0xffff) == 0))
288	SPINNER(1);
289      er = 0;
290      ++*q;
291    } else {
292      er = es * 2;
293      *q += er;
294    }
295  } while (*q != 0);
296  q[1] = 1;
297  SPINNER(2);
298  printf ("\nValues of a single limb, tested.\n");
299  if (justone) return 0;
300  if (quick) {
301    es <<= GMP_NUMB_BITS / 2 - 1;
302    q[1] <<= GMP_NUMB_BITS - 2;
303    printf ("Quick, skipping some... (%u)\n", GMP_NUMB_BITS - 2);
304    --es;
305    --q[1];
306    q[0] -= es*2+1;
307  }
308  printf ("Corner cases tested, up to bits:\n");
309  do {
310    x = mpn_sqrtrem (s, r, q, 2);
311    if (UNLIKELY (x != (er != 0)) || UNLIKELY (*s != es)
312	|| UNLIKELY ((x == 1) && (er != *r)))
313      STOP (something_wrong (er, 0, es));
314
315    if (er != 0) {
316      ++es;
317      if (UNLIKELY ((es & 0xff) == 0))
318	SPINNER(2);
319      er = 0;
320      q[1] += (++*q == 0);
321      if (UNLIKELY (es == GMP_NUMB_HIGHBIT))
322	break;
323    } else {
324      er = es * 2;
325      add_ssaaaa (q[1], *q, q[1], *q, 0, er);
326    }
327  } while (1);
328  SPINNER(2);
329  printf ("\nValues with at most a limb for reminder, tested.\nCorner cases tested, up to bits:\n");
330  x = mpn_sqrtrem (s, r, q, 2);
331  if ((*s != es) || (x != 0))
332    STOP (something_wrong (0, 0, es));
333  q[1] += 1;
334  x = mpn_sqrtrem (s, r, q, 2);
335  if ((*s != es) || (x != 2) || (*r != 0) || (r[1] != 1))
336    STOP (something_wrong (0, 1, es));
337  ++es;
338  q[1] += (++*q == 0);
339  do {
340    x = mpn_sqrtrem (s, r, q, 2);
341    if (UNLIKELY (x != (er != 0) * 2) || UNLIKELY (*s != es)
342	|| UNLIKELY ((x == 2) && ((er != *r) || (r[1] != 1))))
343      STOP (something_wrong (er, er != 0, es));
344
345    if (er != 0) {
346      ++es;
347      if (UNLIKELY (es == 0))
348	break;
349      if (UNLIKELY ((es & 0xff) == 0))
350	SPINNER(2);
351      er = 0;
352      q[1] += (++*q == 0);
353    } else {
354      er = es * 2;
355      add_ssaaaa (q[1], *q, q[1], *q, 1, er);
356    }
357  } while (1);
358  printf ("| %u\nValues of at most two limbs, tested.\n", GMP_NUMB_BITS*2);
359  return 0;
360}
361
362int
363main (int argc, char **argv)
364{
365  int mode = 0;
366  int justone = 0;
367  int quick = 0;
368
369  for (;argc > 1;--argc,++argv)
370    switch (*argv[1]) {
371    default:
372      fprintf (stderr, "usage: sqrtrem_1_2 [x|c|s] [1|2] [q]\n");
373      exit (1);
374    case 'x':
375      mode = 0;
376      break;
377    case 'c':
378      mode = 1;
379      break;
380    case 's':
381      mode = 2;
382      break;
383    case 'q':
384      quick = 1;
385      break;
386    case '1':
387      justone = 1;
388      break;
389    case '2':
390      justone = 0;
391    }
392
393  switch (mode) {
394  default:
395    return check_all_values (justone, quick);
396  case 1:
397    return check_corner_cases (justone, quick);
398  case 2:
399    return check_some_values (justone, quick);
400  }
401}
402