xref: /netbsd-src/external/lgpl3/gmp/dist/tests/devel/sqrtrem_1_2.c (revision 72c7faa4dbb41dbb0238d6b4a109da0d4b236dd4)
1 /*
2 Copyright 2017 Free Software Foundation, Inc.
3 
4 This file is part of the GNU MP Library test suite.
5 
6 The GNU MP Library test suite is free software; you can redistribute it
7 and/or modify it under the terms of the GNU General Public License as
8 published by the Free Software Foundation; either version 3 of the License,
9 or (at your option) any later version.
10 
11 The GNU MP Library test suite is distributed in the hope that it will be
12 useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General
14 Public License for more details.
15 
16 You should have received a copy of the GNU General Public License along with
17 the 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 
something_wrong(mp_limb_t er,mp_limb_t ec,mp_limb_t es)61 int 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 
67 int
check_all_values(int justone,int quick)68 check_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 
128 mp_limb_t
upd(mp_limb_t * s,mp_limb_t k)129 upd (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 
142 mp_limb_t
upd1(mp_limb_t * s,mp_limb_t k)143 upd1 (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 
152 int
check_some_values(int justone,int quick)153 check_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 
264 int
check_corner_cases(int justone,int quick)265 check_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 
362 int
main(int argc,char ** argv)363 main (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