Home | History | Annotate | Line # | Download | only in tests 
      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 ((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 
    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 }
    182