(root)/
gmp-6.3.0/
tests/
devel/
sqrtrem_1_2.c
       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  
      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
      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
     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
     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
     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
     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
     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  }