1/* 2 3Copyright 2012, 2014, Free Software Foundation, Inc. 4 5This file is part of the GNU MP Library test suite. 6 7The GNU MP Library test suite is free software; you can redistribute it 8and/or modify it under the terms of the GNU General Public License as 9published by the Free Software Foundation; either version 3 of the License, 10or (at your option) any later version. 11 12The GNU MP Library test suite is distributed in the hope that it will be 13useful, but WITHOUT ANY WARRANTY; without even the implied warranty of 14MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General 15Public License for more details. 16 17You should have received a copy of the GNU General Public License along with 18the 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 */ 30static int 31sqrtrem_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 55void 56mpz_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 104void 105testmain (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