1 /* Copyright 2006, 2007, 2009, 2010 Free Software Foundation, Inc. 2 3 This program is free software; you can redistribute it and/or modify it under 4 the terms of the GNU General Public License as published by the Free Software 5 Foundation; either version 3 of the License, or (at your option) any later 6 version. 7 8 This program is distributed in the hope that it will be useful, but WITHOUT ANY 9 WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A 10 PARTICULAR PURPOSE. See the GNU General Public License for more details. 11 12 You should have received a copy of the GNU General Public License along with 13 this program. If not, see http://www.gnu.org/licenses/. */ 14 15 16 #include <stdlib.h> /* for strtol */ 17 #include <stdio.h> /* for printf */ 18 19 #include "gmp.h" 20 #include "gmp-impl.h" 21 #include "longlong.h" 22 #include "tests/tests.h" 23 24 25 static void 26 dumpy (mp_srcptr p, mp_size_t n) 27 { 28 mp_size_t i; 29 if (n > 20) 30 { 31 for (i = n - 1; i >= n - 4; i--) 32 { 33 printf ("%0*lx", (int) (2 * sizeof (mp_limb_t)), p[i]); 34 printf (" "); 35 } 36 printf ("... "); 37 for (i = 3; i >= 0; i--) 38 { 39 printf ("%0*lx", (int) (2 * sizeof (mp_limb_t)), p[i]); 40 printf (" " + (i == 0)); 41 } 42 } 43 else 44 { 45 for (i = n - 1; i >= 0; i--) 46 { 47 printf ("%0*lx", (int) (2 * sizeof (mp_limb_t)), p[i]); 48 printf (" " + (i == 0)); 49 } 50 } 51 puts (""); 52 } 53 54 static unsigned long test; 55 56 void 57 check_one (mp_ptr qp, mp_srcptr rp, mp_limb_t rh, 58 mp_srcptr np, mp_size_t nn, mp_srcptr dp, mp_size_t dn, char *fname) 59 { 60 mp_size_t qn; 61 int cmp; 62 mp_ptr tp; 63 mp_limb_t cy = 4711; /* silence warnings */ 64 TMP_DECL; 65 66 qn = nn - dn; 67 68 if (qn == 0) 69 return; 70 71 TMP_MARK; 72 73 tp = TMP_ALLOC_LIMBS (nn + 1); 74 75 if (dn >= qn) 76 mpn_mul (tp, dp, dn, qp, qn); 77 else 78 mpn_mul (tp, qp, qn, dp, dn); 79 80 if (rp != NULL) 81 { 82 cy = mpn_add_n (tp + qn, tp + qn, rp, dn); 83 cmp = cy != rh || mpn_cmp (tp, np, nn) != 0; 84 } 85 else 86 cmp = mpn_cmp (tp, np, nn - dn) != 0; 87 88 if (cmp != 0) 89 { 90 printf ("\r*******************************************************************************\n"); 91 printf ("%s inconsistent in test %lu\n", fname, test); 92 printf ("N= "); dumpy (np, nn); 93 printf ("D= "); dumpy (dp, dn); 94 printf ("Q= "); dumpy (qp, qn); 95 if (rp != NULL) 96 { 97 printf ("R= "); dumpy (rp, dn); 98 printf ("Rb= %d, Cy=%d\n", (int) cy, (int) rh); 99 } 100 printf ("T= "); dumpy (tp, nn); 101 printf ("nn = %ld, dn = %ld, qn = %ld", nn, dn, qn); 102 printf ("\n*******************************************************************************\n"); 103 abort (); 104 } 105 106 TMP_FREE; 107 } 108 109 110 /* These are *bit* sizes. */ 111 #define SIZE_LOG 16 112 #define MAX_DN (1L << SIZE_LOG) 113 #define MAX_NN (1L << (SIZE_LOG + 1)) 114 115 #define COUNT 500 116 117 mp_limb_t 118 random_word (gmp_randstate_ptr rs) 119 { 120 mpz_t x; 121 mp_limb_t r; 122 TMP_DECL; 123 TMP_MARK; 124 125 MPZ_TMP_INIT (x, 2); 126 mpz_urandomb (x, rs, 32); 127 r = mpz_get_ui (x); 128 TMP_FREE; 129 return r; 130 } 131 132 int 133 main (int argc, char **argv) 134 { 135 gmp_randstate_ptr rands; 136 unsigned long maxnbits, maxdbits, nbits, dbits; 137 mpz_t n, d, tz; 138 mp_size_t maxnn, maxdn, nn, dn, clearn, i; 139 mp_ptr np, dp, qp, rp; 140 mp_limb_t rh; 141 mp_limb_t t; 142 mp_limb_t dinv; 143 int count = COUNT; 144 mp_ptr scratch; 145 mp_limb_t ran; 146 mp_size_t alloc, itch; 147 mp_limb_t rran0, rran1, qran0, qran1; 148 TMP_DECL; 149 150 if (argc > 1) 151 { 152 char *end; 153 count = strtol (argv[1], &end, 0); 154 if (*end || count <= 0) 155 { 156 fprintf (stderr, "Invalid test count: %s.\n", argv[1]); 157 return 1; 158 } 159 } 160 161 162 maxdbits = MAX_DN; 163 maxnbits = MAX_NN; 164 165 tests_start (); 166 rands = RANDS; 167 168 mpz_init (n); 169 mpz_init (d); 170 mpz_init (tz); 171 172 maxnn = maxnbits / GMP_NUMB_BITS + 1; 173 maxdn = maxdbits / GMP_NUMB_BITS + 1; 174 175 TMP_MARK; 176 177 qp = TMP_ALLOC_LIMBS (maxnn + 2) + 1; 178 rp = TMP_ALLOC_LIMBS (maxnn + 2) + 1; 179 180 alloc = 1; 181 scratch = __GMP_ALLOCATE_FUNC_LIMBS (alloc); 182 183 for (test = 0; test < count;) 184 { 185 nbits = random_word (rands) % (maxnbits - GMP_NUMB_BITS) + 2 * GMP_NUMB_BITS; 186 if (maxdbits > nbits) 187 dbits = random_word (rands) % nbits + 1; 188 else 189 dbits = random_word (rands) % maxdbits + 1; 190 191 #if RAND_UNIFORM 192 #define RANDFUNC mpz_urandomb 193 #else 194 #define RANDFUNC mpz_rrandomb 195 #endif 196 197 do 198 { 199 RANDFUNC (n, rands, nbits); 200 do 201 { 202 RANDFUNC (d, rands, dbits); 203 } 204 while (mpz_sgn (d) == 0); 205 206 np = PTR (n); 207 dp = PTR (d); 208 nn = SIZ (n); 209 dn = SIZ (d); 210 } 211 while (nn < dn); 212 213 dp[0] |= 1; 214 215 mpz_urandomb (tz, rands, 32); 216 t = mpz_get_ui (tz); 217 218 if (t % 17 == 0) 219 dp[0] = GMP_NUMB_MAX; 220 221 switch ((int) t % 16) 222 { 223 case 0: 224 clearn = random_word (rands) % nn; 225 for (i = 0; i <= clearn; i++) 226 np[i] = 0; 227 break; 228 case 1: 229 mpn_sub_1 (np + nn - dn, dp, dn, random_word (rands)); 230 break; 231 case 2: 232 mpn_add_1 (np + nn - dn, dp, dn, random_word (rands)); 233 break; 234 } 235 236 test++; 237 238 binvert_limb (dinv, dp[0]); 239 240 rran0 = random_word (rands); 241 rran1 = random_word (rands); 242 qran0 = random_word (rands); 243 qran1 = random_word (rands); 244 245 qp[-1] = qran0; 246 qp[nn - dn + 1] = qran1; 247 rp[-1] = rran0; 248 249 ran = random_word (rands); 250 251 if ((double) (nn - dn) * dn < 1e5) 252 { 253 if (nn > dn) 254 { 255 /* Test mpn_sbpi1_bdiv_qr */ 256 MPN_ZERO (qp, nn - dn); 257 MPN_ZERO (rp, dn); 258 MPN_COPY (rp, np, nn); 259 rh = mpn_sbpi1_bdiv_qr (qp, rp, nn, dp, dn, -dinv); 260 ASSERT_ALWAYS (qp[-1] == qran0); ASSERT_ALWAYS (qp[nn - dn + 1] == qran1); 261 ASSERT_ALWAYS (rp[-1] == rran0); 262 check_one (qp, rp + nn - dn, rh, np, nn, dp, dn, "mpn_sbpi1_bdiv_qr"); 263 } 264 265 if (nn > dn) 266 { 267 /* Test mpn_sbpi1_bdiv_q */ 268 MPN_COPY (rp, np, nn); 269 MPN_ZERO (qp, nn - dn); 270 mpn_sbpi1_bdiv_q (qp, rp, nn - dn, dp, MIN(dn,nn-dn), -dinv); 271 ASSERT_ALWAYS (qp[-1] == qran0); ASSERT_ALWAYS (qp[nn - dn + 1] == qran1); 272 ASSERT_ALWAYS (rp[-1] == rran0); 273 check_one (qp, NULL, 0, np, nn, dp, dn, "mpn_sbpi1_bdiv_q"); 274 } 275 } 276 277 if (dn >= 4 && nn - dn >= 2) 278 { 279 /* Test mpn_dcpi1_bdiv_qr */ 280 MPN_COPY (rp, np, nn); 281 MPN_ZERO (qp, nn - dn); 282 rh = mpn_dcpi1_bdiv_qr (qp, rp, nn, dp, dn, -dinv); 283 ASSERT_ALWAYS (qp[-1] == qran0); ASSERT_ALWAYS (qp[nn - dn + 1] == qran1); 284 ASSERT_ALWAYS (rp[-1] == rran0); 285 check_one (qp, rp + nn - dn, rh, np, nn, dp, dn, "mpn_dcpi1_bdiv_qr"); 286 } 287 288 if (dn >= 4 && nn - dn >= 2) 289 { 290 /* Test mpn_dcpi1_bdiv_q */ 291 MPN_COPY (rp, np, nn); 292 MPN_ZERO (qp, nn - dn); 293 mpn_dcpi1_bdiv_q (qp, rp, nn - dn, dp, MIN(dn,nn-dn), -dinv); 294 ASSERT_ALWAYS (qp[-1] == qran0); ASSERT_ALWAYS (qp[nn - dn + 1] == qran1); 295 ASSERT_ALWAYS (rp[-1] == rran0); 296 check_one (qp, NULL, 0, np, nn, dp, dn, "mpn_dcpi1_bdiv_q"); 297 } 298 299 if (nn - dn < 2 || dn < 2) 300 continue; 301 302 /* Test mpn_mu_bdiv_qr */ 303 itch = mpn_mu_bdiv_qr_itch (nn, dn); 304 if (itch + 1 > alloc) 305 { 306 scratch = __GMP_REALLOCATE_FUNC_LIMBS (scratch, alloc, itch + 1); 307 alloc = itch + 1; 308 } 309 scratch[itch] = ran; 310 MPN_ZERO (qp, nn - dn); 311 MPN_ZERO (rp, dn); 312 rp[dn] = rran1; 313 rh = mpn_mu_bdiv_qr (qp, rp, np, nn, dp, dn, scratch); 314 ASSERT_ALWAYS (ran == scratch[itch]); 315 ASSERT_ALWAYS (qp[-1] == qran0); ASSERT_ALWAYS (qp[nn - dn + 1] == qran1); 316 ASSERT_ALWAYS (rp[-1] == rran0); ASSERT_ALWAYS (rp[dn] == rran1); 317 check_one (qp, rp, rh, np, nn, dp, dn, "mpn_mu_bdiv_qr"); 318 319 /* Test mpn_mu_bdiv_q */ 320 itch = mpn_mu_bdiv_q_itch (nn, dn); 321 if (itch + 1 > alloc) 322 { 323 scratch = __GMP_REALLOCATE_FUNC_LIMBS (scratch, alloc, itch + 1); 324 alloc = itch + 1; 325 } 326 scratch[itch] = ran; 327 MPN_ZERO (qp, nn - dn + 1); 328 mpn_mu_bdiv_q (qp, np, nn - dn, dp, dn, scratch); 329 ASSERT_ALWAYS (ran == scratch[itch]); 330 ASSERT_ALWAYS (qp[-1] == qran0); ASSERT_ALWAYS (qp[nn - dn + 1] == qran1); 331 check_one (qp, NULL, 0, np, nn, dp, dn, "mpn_mu_bdiv_q"); 332 } 333 334 __GMP_FREE_FUNC_LIMBS (scratch, alloc); 335 336 TMP_FREE; 337 338 mpz_clear (n); 339 mpz_clear (d); 340 mpz_clear (tz); 341 342 tests_end (); 343 return 0; 344 } 345