1 4 package gnu.math; 5 6 12 13 class MPN 14 { 15 20 public static int add_1 (int[] dest, int[] x, int size, int y) 21 { 22 long carry = (long) y & 0xffffffffL; 23 for (int i = 0; i < size; i++) 24 { 25 carry += ((long) x[i] & 0xffffffffL); 26 dest[i] = (int) carry; 27 carry >>= 32; 28 } 29 return (int) carry; 30 } 31 32 38 public static int add_n (int dest[], int[] x, int[] y, int len) 39 { 40 long carry = 0; 41 for (int i = 0; i < len; i++) 42 { 43 carry += ((long) x[i] & 0xffffffffL) 44 + ((long) y[i] & 0xffffffffL); 45 dest[i] = (int) carry; 46 carry >>>= 32; 47 } 48 return (int) carry; 49 } 50 51 56 57 public static int sub_n (int[] dest, int[] X, int[] Y, int size) 58 { 59 int cy = 0; 60 for (int i = 0; i < size; i++) 61 { 62 int y = Y[i]; 63 int x = X[i]; 64 y += cy; 65 cy = (y^0x80000000) < (cy^0x80000000) ? 1 : 0; 68 y = x - y; 69 cy += (y^0x80000000) > (x ^ 0x80000000) ? 1 : 0; 70 dest[i] = y; 71 } 72 return cy; 73 } 74 75 83 84 public static int mul_1 (int[] dest, int[] x, int len, int y) 85 { 86 long yword = (long) y & 0xffffffffL; 87 long carry = 0; 88 for (int j = 0; j < len; j++) 89 { 90 carry += ((long) x[j] & 0xffffffffL) * yword; 91 dest[j] = (int) carry; 92 carry >>>= 32; 93 } 94 return (int) carry; 95 } 96 97 106 107 public static void mul (int[] dest, 108 int[] x, int xlen, 109 int[] y, int ylen) 110 { 111 dest[xlen] = MPN.mul_1 (dest, x, xlen, y[0]); 112 113 for (int i = 1; i < ylen; i++) 114 { 115 long yword = (long) y[i] & 0xffffffffL; 116 long carry = 0; 117 for (int j = 0; j < xlen; j++) 118 { 119 carry += ((long) x[j] & 0xffffffffL) * yword 120 + ((long) dest[i+j] & 0xffffffffL); 121 dest[i+j] = (int) carry; 122 carry >>>= 32; 123 } 124 dest[i+xlen] = (int) carry; 125 } 126 } 127 128 133 public static long udiv_qrnnd (long N, int D) 134 { 135 long q, r; 136 long a1 = N >>> 32; 137 long a0 = N & 0xffffffffL; 138 if (D >= 0) 139 { 140 if (a1 < ((D - a1 - (a0 >>> 31)) & 0xffffffffL)) 141 { 142 143 q = N / D; 144 r = N % D; 145 } 146 else 147 { 148 149 long c = N - ((long) D << 31); 150 151 q = c / D; 152 r = c % D; 153 154 q += 1 << 31; 155 } 156 } 157 else 158 { 159 long b1 = D >>> 1; 160 long c = N >>> 1; 163 if (a1 < b1 || (a1 >> 1) < b1) 164 { 165 if (a1 < b1) 166 { 167 q = c / b1; 168 r = c % b1; 169 } 170 else 171 { 172 c = ~(c - (b1 << 32)); 173 q = c / b1; 174 r = c % b1; 175 q = (~q) & 0xffffffffL; 176 r = (b1 - 1) - r; 177 } 178 r = 2 * r + (a0 & 1); 179 if ((D & 1) != 0) 180 { 181 if (r >= q) { 182 r = r - q; 183 } else if (q - r <= ((long) D & 0xffffffffL)) { 184 r = r - q + D; 185 q -= 1; 186 } else { 187 r = r - q + D + D; 188 q -= 2; 189 } 190 } 191 } 192 else 193 { 194 if (a0 >= ((long)(-D) & 0xffffffffL)) 195 { 196 q = -1; 197 r = a0 + D; 198 } 199 else 200 { 201 q = -2; 202 r = a0 + D + D; 203 } 204 } 205 } 206 207 return (r << 32) | (q & 0xFFFFFFFFl); 208 } 209 210 215 216 public static int divmod_1 (int[] quotient, int[] dividend, 217 int len, int divisor) 218 { 219 int i = len - 1; 220 long r = dividend[i]; 221 if ((r & 0xffffffffL) >= ((long)divisor & 0xffffffffL)) 222 r = 0; 223 else 224 { 225 quotient[i--] = 0; 226 r <<= 32; 227 } 228 229 for (; i >= 0; i--) 230 { 231 int n0 = dividend[i]; 232 r = (r & ~0xffffffffL) | (n0 & 0xffffffffL); 233 r = udiv_qrnnd (r, divisor); 234 quotient[i] = (int) r; 235 } 236 return (int)(r >> 32); 237 } 238 239 244 public static int submul_1 (int[] dest, int offset, int[] x, int len, int y) 245 { 246 long yl = (long) y & 0xffffffffL; 247 int carry = 0; 248 int j = 0; 249 do 250 { 251 long prod = ((long) x[j] & 0xffffffffL) * yl; 252 int prod_low = (int) prod; 253 int prod_high = (int) (prod >> 32); 254 prod_low += carry; 255 carry = ((prod_low ^ 0x80000000) < (carry ^ 0x80000000) ? 1 : 0) 258 + prod_high; 259 int x_j = dest[offset+j]; 260 prod_low = x_j - prod_low; 261 if ((prod_low ^ 0x80000000) > (x_j ^ 0x80000000)) 262 carry++; 263 dest[offset+j] = prod_low; 264 } 265 while (++j < len); 266 return carry; 267 } 268 269 275 276 public static void divide (int[] zds, int nx, int[] y, int ny) 277 { 278 281 288 291 int j = nx; 292 do 293 { int qhat; if (zds[j]==y[ny-1]) 298 qhat = -1; else 300 { 301 long w = (((long)(zds[j])) << 32) + ((long)zds[j-1] & 0xffffffffL); 302 qhat = (int) udiv_qrnnd (w, y[ny-1]); 303 } 304 if (qhat != 0) 305 { 306 int borrow = submul_1 (zds, j - ny, y, ny, qhat); 307 int save = zds[j]; 308 long num = ((long)save&0xffffffffL) - ((long)borrow&0xffffffffL); 309 while (num != 0) 310 { 311 qhat--; 312 long carry = 0; 313 for (int i = 0; i < ny; i++) 314 { 315 carry += ((long) zds[j-ny+i] & 0xffffffffL) 316 + ((long) y[i] & 0xffffffffL); 317 zds[j-ny+i] = (int) carry; 318 carry >>>= 32; 319 } 320 zds[j] += carry; 321 num = carry - 1; 322 } 323 } 324 zds[j] = qhat; 325 } while (--j >= ny); 326 } 327 328 334 public static int chars_per_word (int radix) 335 { 336 if (radix < 10) 337 { 338 if (radix < 8) 339 { 340 if (radix <= 2) 341 return 32; 342 else if (radix == 3) 343 return 20; 344 else if (radix == 4) 345 return 16; 346 else 347 return 18 - radix; 348 } 349 else 350 return 10; 351 } 352 else if (radix < 12) 353 return 9; 354 else if (radix <= 16) 355 return 8; 356 else if (radix <= 23) 357 return 7; 358 else if (radix <= 40) 359 return 6; 360 else if (radix <= 256) 362 return 4; 363 else 364 return 1; 365 } 366 367 368 public static int count_leading_zeros (int i) 369 { 370 if (i == 0) 371 return 32; 372 int count = 0; 373 for (int k = 16; k > 0; k = k >> 1) { 374 int j = i >>> k; 375 if (j == 0) 376 count += k; 377 else 378 i = j; 379 } 380 return count; 381 } 382 383 public static int set_str (int dest[], byte[] str, int str_len, int base) 384 { 385 int size = 0; 386 if ((base & (base - 1)) == 0) 387 { 388 391 int next_bitpos = 0; 392 int bits_per_indigit = 0; 393 for (int i = base; (i >>= 1) != 0; ) bits_per_indigit++; 394 int res_digit = 0; 395 396 for (int i = str_len; --i >= 0; ) 397 { 398 int inp_digit = str[i]; 399 res_digit |= inp_digit << next_bitpos; 400 next_bitpos += bits_per_indigit; 401 if (next_bitpos >= 32) 402 { 403 dest[size++] = res_digit; 404 next_bitpos -= 32; 405 res_digit = inp_digit >> (bits_per_indigit - next_bitpos); 406 } 407 } 408 409 if (res_digit != 0) 410 dest[size++] = res_digit; 411 } 412 else 413 { 414 int indigits_per_limb = MPN.chars_per_word (base); 416 int str_pos = 0; 417 418 while (str_pos < str_len) 419 { 420 int chunk = str_len - str_pos; 421 if (chunk > indigits_per_limb) 422 chunk = indigits_per_limb; 423 int res_digit = str[str_pos++]; 424 int big_base = base; 425 426 while (--chunk > 0) 427 { 428 res_digit = res_digit * base + str[str_pos++]; 429 big_base *= base; 430 } 431 432 int cy_limb; 433 if (size == 0) 434 cy_limb = res_digit; 435 else 436 { 437 cy_limb = MPN.mul_1 (dest, dest, size, big_base); 438 cy_limb += MPN.add_1 (dest, dest, size, res_digit); 439 } 440 if (cy_limb != 0) 441 dest[size++] = cy_limb; 442 } 443 } 444 return size; 445 } 446 447 451 public static int cmp (int[] x, int[] y, int size) 452 { 453 while (--size >= 0) 454 { 455 int x_word = x[size]; 456 int y_word = y[size]; 457 if (x_word != y_word) 458 { 459 return (x_word ^ 0x80000000) > (y_word ^0x80000000) ? 1 : -1; 463 } 464 } 465 return 0; 466 } 467 468 472 public static int cmp (int[] x, int xlen, int[] y, int ylen) 473 { 474 return xlen > ylen ? 1 : xlen < ylen ? -1 : cmp (x, y, xlen); 475 } 476 477 484 485 public static int rshift (int[] dest, int[] x, int x_start, 486 int len, int count) 487 { 488 int count_2 = 32 - count; 489 int low_word = x[x_start]; 490 int retval = low_word << count_2; 491 int i = 1; 492 for (; i < len; i++) 493 { 494 int high_word = x[x_start+i]; 495 dest[i-1] = (low_word >>> count) | (high_word << count_2); 496 low_word = high_word; 497 } 498 dest[i-1] = low_word >>> count; 499 return retval; 500 } 501 502 509 public static void rshift0 (int[] dest, int[] x, int x_start, 510 int len, int count) 511 { 512 if (count > 0) 513 rshift(dest, x, x_start, len, count); 514 else 515 for (int i = 0; i < len; i++) 516 dest[i] = x[i + x_start]; 517 } 518 519 525 public static long rshift_long (int[] x, int len, int count) 526 { 527 int wordno = count >> 5; 528 count &= 31; 529 int sign = x[len-1] < 0 ? -1 : 0; 530 int w0 = wordno >= len ? sign : x[wordno]; 531 wordno++; 532 int w1 = wordno >= len ? sign : x[wordno]; 533 if (count != 0) 534 { 535 wordno++; 536 int w2 = wordno >= len ? sign : x[wordno]; 537 w0 = (w0 >>> count) | (w1 << (32-count)); 538 w1 = (w1 >>> count) | (w2 << (32-count)); 539 } 540 return ((long)w1 << 32) | ((long)w0 & 0xffffffffL); 541 } 542 543 549 550 public static int lshift (int[] dest, int d_offset, 551 int[] x, int len, int count) 552 { 553 int count_2 = 32 - count; 554 int i = len - 1; 555 int high_word = x[i]; 556 int retval = high_word >>> count_2; 557 d_offset++; 558 while (--i >= 0) 559 { 560 int low_word = x[i]; 561 dest[d_offset+i] = (high_word << count) | (low_word >>> count_2); 562 high_word = low_word; 563 } 564 dest[d_offset+i] = high_word << count; 565 return retval; 566 } 567 568 569 570 static int findLowestBit (int word) 571 { 572 int i = 0; 573 while ((word & 0xF) == 0) 574 { 575 word >>= 4; 576 i += 4; 577 } 578 if ((word & 3) == 0) 579 { 580 word >>= 2; 581 i += 2; 582 } 583 if ((word & 1) == 0) 584 i += 1; 585 return i; 586 } 587 588 589 590 static int findLowestBit (int[] words) 591 { 592 for (int i = 0; ; i++) 593 { 594 if (words[i] != 0) 595 return 32 * i + findLowestBit (words[i]); 596 } 597 } 598 599 603 604 public static int gcd (int[] x, int[] y, int len) 605 { 606 int i, word; 607 for (i = 0; ; i++) 609 { 610 word = x[i] | y[i]; 611 if (word != 0) 612 { 613 break; 615 } 616 } 617 int initShiftWords = i; 618 int initShiftBits = findLowestBit (word); 619 621 len -= initShiftWords; 623 MPN.rshift0 (x, x, initShiftWords, len, initShiftBits); 624 MPN.rshift0 (y, y, initShiftWords, len, initShiftBits); 625 626 int[] odd_arg; 627 int[] other_arg; 628 if ((x[0] & 1) != 0) 629 { 630 odd_arg = x; 631 other_arg = y; 632 } 633 else 634 { 635 odd_arg = y; 636 other_arg = x; 637 } 638 639 for (;;) 640 { 641 for (i = 0; other_arg[i] == 0; ) i++; 645 if (i > 0) 646 { 647 int j; 648 for (j = 0; j < len-i; j++) 649 other_arg[j] = other_arg[j+i]; 650 for ( ; j < len; j++) 651 other_arg[j] = 0; 652 } 653 i = findLowestBit(other_arg[0]); 654 if (i > 0) 655 MPN.rshift (other_arg, other_arg, 0, len, i); 656 657 659 i = MPN.cmp(odd_arg, other_arg, len); 662 if (i == 0) 663 break; 664 if (i > 0) 665 { MPN.sub_n (odd_arg, odd_arg, other_arg, len); 667 int[] tmp = odd_arg; odd_arg = other_arg; other_arg = tmp; 669 } 670 else 671 { MPN.sub_n (other_arg, other_arg, odd_arg, len); 673 } 674 while (odd_arg[len-1] == 0 && other_arg[len-1] == 0) 675 len--; 676 } 677 if (initShiftWords + initShiftBits > 0) 678 { 679 if (initShiftBits > 0) 680 { 681 int sh_out = MPN.lshift (x, initShiftWords, x, len, initShiftBits); 682 if (sh_out != 0) 683 x[(len++)+initShiftWords] = sh_out; 684 } 685 else 686 { 687 for (i = len; --i >= 0;) 688 x[i+initShiftWords] = x[i]; 689 } 690 for (i = initShiftWords; --i >= 0; ) 691 x[i] = 0; 692 len += initShiftWords; 693 } 694 return len; 695 } 696 697 public static int intLength (int i) 698 { 699 return 32 - count_leading_zeros (i < 0 ? ~i : i); 700 } 701 702 704 public static int intLength (int[] words, int len) 705 { 706 len--; 707 return intLength (words[len]) + 32 * len; 708 } 709 710 732 } 733 | Popular Tags |