github.com/aergoio/aergo@v1.3.1/libtool/src/gmp-6.1.2/mini-gmp/tests/t-sqrt.c (about)

     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 (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  }