Source for gnu.java.math.MPN

   1: /* gnu.java.math.MPN
   2:    Copyright (C) 1999, 2000, 2001, 2004  Free Software Foundation, Inc.
   3: 
   4: This file is part of GNU Classpath.
   5: 
   6: GNU Classpath is free software; you can redistribute it and/or modify
   7: it under the terms of the GNU General Public License as published by
   8: the Free Software Foundation; either version 2, or (at your option)
   9: any later version.
  10: 
  11: GNU Classpath is distributed in the hope that it will be useful, but
  12: WITHOUT ANY WARRANTY; without even the implied warranty of
  13: MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14: General Public License for more details.
  15: 
  16: You should have received a copy of the GNU General Public License
  17: along with GNU Classpath; see the file COPYING.  If not, write to the
  18: Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
  19: 02110-1301 USA.
  20: 
  21: Linking this library statically or dynamically with other modules is
  22: making a combined work based on this library.  Thus, the terms and
  23: conditions of the GNU General Public License cover the whole
  24: combination.
  25: 
  26: As a special exception, the copyright holders of this library give you
  27: permission to link this library with independent modules to produce an
  28: executable, regardless of the license terms of these independent
  29: modules, and to copy and distribute the resulting executable under
  30: terms of your choice, provided that you also meet, for each linked
  31: independent module, the terms and conditions of the license of that
  32: module.  An independent module is a module which is not derived from
  33: or based on this library.  If you modify this library, you may extend
  34: this exception to your version of the library, but you are not
  35: obligated to do so.  If you do not wish to do so, delete this
  36: exception statement from your version. */
  37: 
  38: // Included from Kawa 1.6.62 with permission of the author,
  39: // Per Bothner <per@bothner.com>.
  40: 
  41: package gnu.java.math;
  42: 
  43: /** This contains various low-level routines for unsigned bigints.
  44:  * The interfaces match the mpn interfaces in gmp,
  45:  * so it should be easy to replace them with fast native functions
  46:  * that are trivial wrappers around the mpn_ functions in gmp
  47:  * (at least on platforms that use 32-bit "limbs").
  48:  */
  49: 
  50: public class MPN
  51: {
  52:   /** Add x[0:size-1] and y, and write the size least
  53:    * significant words of the result to dest.
  54:    * Return carry, either 0 or 1.
  55:    * All values are unsigned.
  56:    * This is basically the same as gmp's mpn_add_1. */
  57:   public static int add_1 (int[] dest, int[] x, int size, int y)
  58:   {
  59:     long carry = (long) y & 0xffffffffL;
  60:     for (int i = 0;  i < size;  i++)
  61:       {
  62:         carry += ((long) x[i] & 0xffffffffL);
  63:         dest[i] = (int) carry;
  64:         carry >>= 32;
  65:       }
  66:     return (int) carry;
  67:   }
  68: 
  69:   /** Add x[0:len-1] and y[0:len-1] and write the len least
  70:    * significant words of the result to dest[0:len-1].
  71:    * All words are treated as unsigned.
  72:    * @return the carry, either 0 or 1
  73:    * This function is basically the same as gmp's mpn_add_n.
  74:    */
  75:   public static int add_n (int dest[], int[] x, int[] y, int len)
  76:   {
  77:     long carry = 0;
  78:     for (int i = 0; i < len;  i++)
  79:       {
  80:         carry += ((long) x[i] & 0xffffffffL)
  81:           + ((long) y[i] & 0xffffffffL);
  82:         dest[i] = (int) carry;
  83:         carry >>>= 32;
  84:       }
  85:     return (int) carry;
  86:   }
  87: 
  88:   /** Subtract Y[0:size-1] from X[0:size-1], and write
  89:    * the size least significant words of the result to dest[0:size-1].
  90:    * Return borrow, either 0 or 1.
  91:    * This is basically the same as gmp's mpn_sub_n function.
  92:    */
  93: 
  94:   public static int sub_n (int[] dest, int[] X, int[] Y, int size)
  95:   {
  96:     int cy = 0;
  97:     for (int i = 0;  i < size;  i++)
  98:       {
  99:         int y = Y[i];
 100:         int x = X[i];
 101:         y += cy;        /* add previous carry to subtrahend */
 102:         // Invert the high-order bit, because: (unsigned) X > (unsigned) Y
 103:         // iff: (int) (X^0x80000000) > (int) (Y^0x80000000).
 104:         cy = (y^0x80000000) < (cy^0x80000000) ? 1 : 0;
 105:         y = x - y;
 106:         cy += (y^0x80000000) > (x ^ 0x80000000) ? 1 : 0;
 107:         dest[i] = y;
 108:       }
 109:     return cy;
 110:   }
 111: 
 112:   /** Multiply x[0:len-1] by y, and write the len least
 113:    * significant words of the product to dest[0:len-1].
 114:    * Return the most significant word of the product.
 115:    * All values are treated as if they were unsigned
 116:    * (i.e. masked with 0xffffffffL).
 117:    * OK if dest==x (not sure if this is guaranteed for mpn_mul_1).
 118:    * This function is basically the same as gmp's mpn_mul_1.
 119:    */
 120: 
 121:   public static int mul_1 (int[] dest, int[] x, int len, int y)
 122:   {
 123:     long yword = (long) y & 0xffffffffL;
 124:     long carry = 0;
 125:     for (int j = 0;  j < len; j++)
 126:       {
 127:         carry += ((long) x[j] & 0xffffffffL) * yword;
 128:         dest[j] = (int) carry;
 129:         carry >>>= 32;
 130:       }
 131:     return (int) carry;
 132:   }
 133: 
 134:   /**
 135:    * Multiply x[0:xlen-1] and y[0:ylen-1], and
 136:    * write the result to dest[0:xlen+ylen-1].
 137:    * The destination has to have space for xlen+ylen words,
 138:    * even if the result might be one limb smaller.
 139:    * This function requires that xlen >= ylen.
 140:    * The destination must be distinct from either input operands.
 141:    * All operands are unsigned.
 142:    * This function is basically the same gmp's mpn_mul. */
 143: 
 144:   public static void mul (int[] dest,
 145:                           int[] x, int xlen,
 146:                           int[] y, int ylen)
 147:   {
 148:     dest[xlen] = MPN.mul_1 (dest, x, xlen, y[0]);
 149: 
 150:     for (int i = 1;  i < ylen; i++)
 151:       {
 152:         long yword = (long) y[i] & 0xffffffffL;
 153:         long carry = 0;
 154:         for (int j = 0;  j < xlen; j++)
 155:           {
 156:             carry += ((long) x[j] & 0xffffffffL) * yword
 157:               + ((long) dest[i+j] & 0xffffffffL);
 158:             dest[i+j] = (int) carry;
 159:             carry >>>= 32;
 160:           }
 161:         dest[i+xlen] = (int) carry;
 162:       }
 163:   }
 164: 
 165:   /* Divide (unsigned long) N by (unsigned int) D.
 166:    * Returns (remainder << 32)+(unsigned int)(quotient).
 167:    * Assumes (unsigned int)(N>>32) < (unsigned int)D.
 168:    * Code transcribed from gmp-2.0's mpn_udiv_w_sdiv function.
 169:    */
 170:   public static long udiv_qrnnd (long N, int D)
 171:   {
 172:     long q, r;
 173:     long a1 = N >>> 32;
 174:     long a0 = N & 0xffffffffL;
 175:     if (D >= 0)
 176:       {
 177:         if (a1 < ((D - a1 - (a0 >>> 31)) & 0xffffffffL))
 178:           {
 179:             /* dividend, divisor, and quotient are nonnegative */
 180:             q = N / D;
 181:             r = N % D;
 182:           }
 183:         else
 184:           {
 185:             /* Compute c1*2^32 + c0 = a1*2^32 + a0 - 2^31*d */
 186:             long c = N - ((long) D << 31);
 187:             /* Divide (c1*2^32 + c0) by d */
 188:             q = c / D;
 189:             r = c % D;
 190:             /* Add 2^31 to quotient */
 191:             q += 1 << 31;
 192:           }
 193:       }
 194:     else
 195:       {
 196:         long b1 = D >>> 1;      /* d/2, between 2^30 and 2^31 - 1 */
 197:         //long c1 = (a1 >> 1); /* A/2 */
 198:         //int c0 = (a1 << 31) + (a0 >> 1);
 199:         long c = N >>> 1;
 200:         if (a1 < b1 || (a1 >> 1) < b1)
 201:           {
 202:             if (a1 < b1)
 203:               {
 204:                 q = c / b1;
 205:                 r = c % b1;
 206:               }
 207:             else /* c1 < b1, so 2^31 <= (A/2)/b1 < 2^32 */
 208:               {
 209:                 c = ~(c - (b1 << 32));
 210:                 q = c / b1;  /* (A/2) / (d/2) */
 211:                 r = c % b1;
 212:                 q = (~q) & 0xffffffffL;    /* (A/2)/b1 */
 213:                 r = (b1 - 1) - r; /* r < b1 => new r >= 0 */
 214:               }
 215:             r = 2 * r + (a0 & 1);
 216:             if ((D & 1) != 0)
 217:               {
 218:                 if (r >= q) {
 219:                         r = r - q;
 220:                 } else if (q - r <= ((long) D & 0xffffffffL)) {
 221:                        r = r - q + D;
 222:                         q -= 1;
 223:                 } else {
 224:                        r = r - q + D + D;
 225:                         q -= 2;
 226:                 }
 227:               }
 228:           }
 229:         else                            /* Implies c1 = b1 */
 230:           {                             /* Hence a1 = d - 1 = 2*b1 - 1 */
 231:             if (a0 >= ((long)(-D) & 0xffffffffL))
 232:               {
 233:                 q = -1;
 234:                 r = a0 + D;
 235:               }
 236:             else
 237:               {
 238:                 q = -2;
 239:                 r = a0 + D + D;
 240:               }
 241:           }
 242:       }
 243: 
 244:     return (r << 32) | (q & 0xFFFFFFFFl);
 245:   }
 246: 
 247:     /** Divide divident[0:len-1] by (unsigned int)divisor.
 248:      * Write result into quotient[0:len-1.
 249:      * Return the one-word (unsigned) remainder.
 250:      * OK for quotient==dividend.
 251:      */
 252: 
 253:   public static int divmod_1 (int[] quotient, int[] dividend,
 254:                               int len, int divisor)
 255:   {
 256:     int i = len - 1;
 257:     long r = dividend[i];
 258:     if ((r & 0xffffffffL) >= ((long)divisor & 0xffffffffL))
 259:       r = 0;
 260:     else
 261:       {
 262:         quotient[i--] = 0;
 263:         r <<= 32;
 264:       }
 265: 
 266:     for (;  i >= 0;  i--)
 267:       {
 268:         int n0 = dividend[i];
 269:         r = (r & ~0xffffffffL) | (n0 & 0xffffffffL);
 270:         r = udiv_qrnnd (r, divisor);
 271:         quotient[i] = (int) r;
 272:       }
 273:     return (int)(r >> 32);
 274:   }
 275: 
 276:   /* Subtract x[0:len-1]*y from dest[offset:offset+len-1].
 277:    * All values are treated as if unsigned.
 278:    * @return the most significant word of
 279:    * the product, minus borrow-out from the subtraction.
 280:    */
 281:   public static int submul_1 (int[] dest, int offset, int[] x, int len, int y)
 282:   {
 283:     long yl = (long) y & 0xffffffffL;
 284:     int carry = 0;
 285:     int j = 0;
 286:     do
 287:       {
 288:         long prod = ((long) x[j] & 0xffffffffL) * yl;
 289:         int prod_low = (int) prod;
 290:         int prod_high = (int) (prod >> 32);
 291:         prod_low += carry;
 292:         // Invert the high-order bit, because: (unsigned) X > (unsigned) Y
 293:         // iff: (int) (X^0x80000000) > (int) (Y^0x80000000).
 294:         carry = ((prod_low ^ 0x80000000) < (carry ^ 0x80000000) ? 1 : 0)
 295:           + prod_high;
 296:         int x_j = dest[offset+j];
 297:         prod_low = x_j - prod_low;
 298:         if ((prod_low ^ 0x80000000) > (x_j ^ 0x80000000))
 299:           carry++;
 300:         dest[offset+j] = prod_low;
 301:       }
 302:     while (++j < len);
 303:     return carry;
 304:   }
 305: 
 306:   /** Divide zds[0:nx] by y[0:ny-1].
 307:    * The remainder ends up in zds[0:ny-1].
 308:    * The quotient ends up in zds[ny:nx].
 309:    * Assumes:  nx>ny.
 310:    * (int)y[ny-1] < 0  (i.e. most significant bit set)
 311:    */
 312: 
 313:   public static void divide (int[] zds, int nx, int[] y, int ny)
 314:   {
 315:     // This is basically Knuth's formulation of the classical algorithm,
 316:     // but translated from in scm_divbigbig in Jaffar's SCM implementation.
 317: 
 318:     // Correspondance with Knuth's notation:
 319:     // Knuth's u[0:m+n] == zds[nx:0].
 320:     // Knuth's v[1:n] == y[ny-1:0]
 321:     // Knuth's n == ny.
 322:     // Knuth's m == nx-ny.
 323:     // Our nx == Knuth's m+n.
 324: 
 325:     // Could be re-implemented using gmp's mpn_divrem:
 326:     // zds[nx] = mpn_divrem (&zds[ny], 0, zds, nx, y, ny).
 327: 
 328:     int j = nx;
 329:     do
 330:       {                          // loop over digits of quotient
 331:         // Knuth's j == our nx-j.
 332:         // Knuth's u[j:j+n] == our zds[j:j-ny].
 333:         int qhat;  // treated as unsigned
 334:         if (zds[j]==y[ny-1])
 335:           qhat = -1;  // 0xffffffff
 336:         else
 337:           {
 338:             long w = (((long)(zds[j])) << 32) + ((long)zds[j-1] & 0xffffffffL);
 339:             qhat = (int) udiv_qrnnd (w, y[ny-1]);
 340:           }
 341:         if (qhat != 0)
 342:           {
 343:             int borrow = submul_1 (zds, j - ny, y, ny, qhat);
 344:             int save = zds[j];
 345:             long num = ((long)save&0xffffffffL) - ((long)borrow&0xffffffffL);
 346:             while (num != 0)
 347:               {
 348:                 qhat--;
 349:                 long carry = 0;
 350:                 for (int i = 0;  i < ny; i++)
 351:                   {
 352:                     carry += ((long) zds[j-ny+i] & 0xffffffffL)
 353:                       + ((long) y[i] & 0xffffffffL);
 354:                     zds[j-ny+i] = (int) carry;
 355:                     carry >>>= 32;
 356:                   }
 357:                 zds[j] += carry;
 358:                 num = carry - 1;
 359:               }
 360:           }
 361:         zds[j] = qhat;
 362:       } while (--j >= ny);
 363:   }
 364: 
 365:   /** Number of digits in the conversion base that always fits in a word.
 366:    * For example, for base 10 this is 9, since 10**9 is the
 367:    * largest number that fits into a words (assuming 32-bit words).
 368:    * This is the same as gmp's __mp_bases[radix].chars_per_limb.
 369:    * @param radix the base
 370:    * @return number of digits */
 371:   public static int chars_per_word (int radix)
 372:   {
 373:     if (radix < 10)
 374:       {
 375:         if (radix < 8)
 376:           {
 377:             if (radix <= 2)
 378:               return 32;
 379:             else if (radix == 3)
 380:               return 20;
 381:             else if (radix == 4)
 382:               return 16;
 383:             else
 384:               return 18 - radix;
 385:           }
 386:         else
 387:           return 10;
 388:       }
 389:     else if (radix < 12)
 390:       return 9;
 391:     else if (radix <= 16)
 392:       return 8;
 393:     else if (radix <= 23)
 394:       return 7;
 395:     else if (radix <= 40)
 396:       return 6;
 397:     // The following are conservative, but we don't care.
 398:     else if (radix <= 256)
 399:       return 4;
 400:     else
 401:       return 1;
 402:   }
 403: 
 404:   /** Count the number of leading zero bits in an int. */
 405:   public static int count_leading_zeros (int i)
 406:   {
 407:     if (i == 0)
 408:       return 32;
 409:     int count = 0;
 410:     for (int k = 16;  k > 0;  k = k >> 1) {
 411:       int j = i >>> k;
 412:       if (j == 0)
 413:         count += k;
 414:       else
 415:         i = j;
 416:     }
 417:     return count;
 418:   }
 419: 
 420:   public static int set_str (int dest[], byte[] str, int str_len, int base)
 421:   {
 422:     int size = 0;
 423:     if ((base & (base - 1)) == 0)
 424:       {
 425:         // The base is a power of 2.  Read the input string from
 426:         // least to most significant character/digit.  */
 427: 
 428:         int next_bitpos = 0;
 429:         int bits_per_indigit = 0;
 430:         for (int i = base; (i >>= 1) != 0; ) bits_per_indigit++;
 431:         int res_digit = 0;
 432: 
 433:         for (int i = str_len;  --i >= 0; )
 434:           {
 435:             int inp_digit = str[i];
 436:             res_digit |= inp_digit << next_bitpos;
 437:             next_bitpos += bits_per_indigit;
 438:             if (next_bitpos >= 32)
 439:               {
 440:                 dest[size++] = res_digit;
 441:                 next_bitpos -= 32;
 442:                 res_digit = inp_digit >> (bits_per_indigit - next_bitpos);
 443:               }
 444:           }
 445: 
 446:         if (res_digit != 0)
 447:           dest[size++] = res_digit;
 448:       }
 449:     else
 450:       {
 451:         // General case.  The base is not a power of 2.
 452:         int indigits_per_limb = MPN.chars_per_word (base);
 453:         int str_pos = 0;
 454: 
 455:         while (str_pos < str_len)
 456:           {
 457:             int chunk = str_len - str_pos;
 458:             if (chunk > indigits_per_limb)
 459:               chunk = indigits_per_limb;
 460:             int res_digit = str[str_pos++];
 461:             int big_base = base;
 462: 
 463:             while (--chunk > 0)
 464:               {
 465:                 res_digit = res_digit * base + str[str_pos++];
 466:                 big_base *= base;
 467:               }
 468: 
 469:             int cy_limb;
 470:             if (size == 0)
 471:               cy_limb = res_digit;
 472:             else
 473:               {
 474:                 cy_limb = MPN.mul_1 (dest, dest, size, big_base);
 475:                 cy_limb += MPN.add_1 (dest, dest, size, res_digit);
 476:               }
 477:             if (cy_limb != 0)
 478:               dest[size++] = cy_limb;
 479:           }
 480:        }
 481:     return size;
 482:   }
 483: 
 484:   /** Compare x[0:size-1] with y[0:size-1], treating them as unsigned integers.
 485:    * @result -1, 0, or 1 depending on if x&lt;y, x==y, or x&gt;y.
 486:    * This is basically the same as gmp's mpn_cmp function.
 487:    */
 488:   public static int cmp (int[] x, int[] y, int size)
 489:   {
 490:     while (--size >= 0)
 491:       {
 492:         int x_word = x[size];
 493:         int y_word = y[size];
 494:         if (x_word != y_word)
 495:           {
 496:             // Invert the high-order bit, because:
 497:             // (unsigned) X > (unsigned) Y iff
 498:             // (int) (X^0x80000000) > (int) (Y^0x80000000).
 499:             return (x_word ^ 0x80000000) > (y_word ^0x80000000) ? 1 : -1;
 500:           }
 501:       }
 502:     return 0;
 503:   }
 504: 
 505:   /**
 506:    * Compare x[0:xlen-1] with y[0:ylen-1], treating them as unsigned integers.
 507:    *
 508:    * @return -1, 0, or 1 depending on if x&lt;y, x==y, or x&gt;y.
 509:    */
 510:   public static int cmp (int[] x, int xlen, int[] y, int ylen)
 511:   {
 512:     return xlen > ylen ? 1 : xlen < ylen ? -1 : cmp (x, y, xlen);
 513:   }
 514: 
 515:   /**
 516:    * Shift x[x_start:x_start+len-1] count bits to the "right"
 517:    * (i.e. divide by 2**count).
 518:    * Store the len least significant words of the result at dest.
 519:    * The bits shifted out to the right are returned.
 520:    * OK if dest==x.
 521:    * Assumes: 0 &lt; count &lt; 32
 522:    */
 523:   public static int rshift (int[] dest, int[] x, int x_start,
 524:                             int len, int count)
 525:   {
 526:     int count_2 = 32 - count;
 527:     int low_word = x[x_start];
 528:     int retval = low_word << count_2;
 529:     int i = 1;
 530:     for (; i < len;  i++)
 531:       {
 532:         int high_word = x[x_start+i];
 533:         dest[i-1] = (low_word >>> count) | (high_word << count_2);
 534:         low_word = high_word;
 535:       }
 536:     dest[i-1] = low_word >>> count;
 537:     return retval;
 538:   }
 539: 
 540:   /**
 541:    * Shift x[x_start:x_start+len-1] count bits to the "right"
 542:    * (i.e. divide by 2**count).
 543:    * Store the len least significant words of the result at dest.
 544:    * OK if dest==x.
 545:    * Assumes: 0 &lt;= count &lt; 32
 546:    * Same as rshift, but handles count==0 (and has no return value).
 547:    */
 548:   public static void rshift0 (int[] dest, int[] x, int x_start,
 549:                               int len, int count)
 550:   {
 551:     if (count > 0)
 552:       rshift(dest, x, x_start, len, count);
 553:     else
 554:       for (int i = 0;  i < len;  i++)
 555:         dest[i] = x[i + x_start];
 556:   }
 557: 
 558:   /** Return the long-truncated value of right shifting.
 559:   * @param x a two's-complement "bignum"
 560:   * @param len the number of significant words in x
 561:   * @param count the shift count
 562:   * @return (long)(x[0..len-1] &gt;&gt; count).
 563:   */
 564:   public static long rshift_long (int[] x, int len, int count)
 565:   {
 566:     int wordno = count >> 5;
 567:     count &= 31;
 568:     int sign = x[len-1] < 0 ? -1 : 0;
 569:     int w0 = wordno >= len ? sign : x[wordno];
 570:     wordno++;
 571:     int w1 = wordno >= len ? sign : x[wordno];
 572:     if (count != 0)
 573:       {
 574:         wordno++;
 575:         int w2 = wordno >= len ? sign : x[wordno];
 576:         w0 = (w0 >>> count) | (w1 << (32-count));
 577:         w1 = (w1 >>> count) | (w2 << (32-count));
 578:       }
 579:     return ((long)w1 << 32) | ((long)w0 & 0xffffffffL);
 580:   }
 581: 
 582:   /* Shift x[0:len-1] left by count bits, and store the len least
 583:    * significant words of the result in dest[d_offset:d_offset+len-1].
 584:    * Return the bits shifted out from the most significant digit.
 585:    * Assumes 0 &lt; count &lt; 32.
 586:    * OK if dest==x.
 587:    */
 588: 
 589:   public static int lshift (int[] dest, int d_offset,
 590:                             int[] x, int len, int count)
 591:   {
 592:     int count_2 = 32 - count;
 593:     int i = len - 1;
 594:     int high_word = x[i];
 595:     int retval = high_word >>> count_2;
 596:     d_offset++;
 597:     while (--i >= 0)
 598:       {
 599:         int low_word = x[i];
 600:         dest[d_offset+i] = (high_word << count) | (low_word >>> count_2);
 601:         high_word = low_word;
 602:       }
 603:     dest[d_offset+i] = high_word << count;
 604:     return retval;
 605:   }
 606: 
 607:   /** Return least i such that word &amp; (1&lt;&lt;i). Assumes word!=0. */
 608: 
 609:   public static int findLowestBit (int word)
 610:   {
 611:     int i = 0;
 612:     while ((word & 0xF) == 0)
 613:       {
 614:         word >>= 4;
 615:         i += 4;
 616:       }
 617:     if ((word & 3) == 0)
 618:       {
 619:         word >>= 2;
 620:         i += 2;
 621:       }
 622:     if ((word & 1) == 0)
 623:       i += 1;
 624:     return i;
 625:   }
 626: 
 627:   /** Return least i such that words &amp; (1&lt;&lt;i). Assumes there is such an i. */
 628: 
 629:   public static int findLowestBit (int[] words)
 630:   {
 631:     for (int i = 0;  ; i++)
 632:       {
 633:         if (words[i] != 0)
 634:           return 32 * i + findLowestBit (words[i]);
 635:       }
 636:   }
 637: 
 638:   /** Calculate Greatest Common Divisior of x[0:len-1] and y[0:len-1].
 639:     * Assumes both arguments are non-zero.
 640:     * Leaves result in x, and returns len of result.
 641:     * Also destroys y (actually sets it to a copy of the result). */
 642: 
 643:   public static int gcd (int[] x, int[] y, int len)
 644:   {
 645:     int i, word;
 646:     // Find sh such that both x and y are divisible by 2**sh.
 647:     for (i = 0; ; i++)
 648:       {
 649:         word = x[i] | y[i];
 650:         if (word != 0)
 651:           {
 652:             // Must terminate, since x and y are non-zero.
 653:             break;
 654:           }
 655:       }
 656:     int initShiftWords = i;
 657:     int initShiftBits = findLowestBit (word);
 658:     // Logically: sh = initShiftWords * 32 + initShiftBits
 659: 
 660:     // Temporarily devide both x and y by 2**sh.
 661:     len -= initShiftWords;
 662:     MPN.rshift0 (x, x, initShiftWords, len, initShiftBits);
 663:     MPN.rshift0 (y, y, initShiftWords, len, initShiftBits);
 664: 
 665:     int[] odd_arg; /* One of x or y which is odd. */
 666:     int[] other_arg; /* The other one can be even or odd. */
 667:     if ((x[0] & 1) != 0)
 668:       {
 669:         odd_arg = x;
 670:         other_arg = y;
 671:       }
 672:     else
 673:       {
 674:         odd_arg = y;
 675:         other_arg = x;
 676:       }
 677: 
 678:     for (;;)
 679:       {
 680:         // Shift other_arg until it is odd; this doesn't
 681:         // affect the gcd, since we divide by 2**k, which does not
 682:         // divide odd_arg.
 683:         for (i = 0; other_arg[i] == 0; ) i++;
 684:         if (i > 0)
 685:           {
 686:             int j;
 687:             for (j = 0; j < len-i; j++)
 688:                 other_arg[j] = other_arg[j+i];
 689:             for ( ; j < len; j++)
 690:               other_arg[j] = 0;
 691:           }
 692:         i = findLowestBit(other_arg[0]);
 693:         if (i > 0)
 694:           MPN.rshift (other_arg, other_arg, 0, len, i);
 695: 
 696:         // Now both odd_arg and other_arg are odd.
 697: 
 698:         // Subtract the smaller from the larger.
 699:         // This does not change the result, since gcd(a-b,b)==gcd(a,b).
 700:         i = MPN.cmp(odd_arg, other_arg, len);
 701:         if (i == 0)
 702:             break;
 703:         if (i > 0)
 704:           { // odd_arg > other_arg
 705:             MPN.sub_n (odd_arg, odd_arg, other_arg, len);
 706:             // Now odd_arg is even, so swap with other_arg;
 707:             int[] tmp = odd_arg; odd_arg = other_arg; other_arg = tmp;
 708:           }
 709:         else
 710:           { // other_arg > odd_arg
 711:             MPN.sub_n (other_arg, other_arg, odd_arg, len);
 712:         }
 713:         while (odd_arg[len-1] == 0 && other_arg[len-1] == 0)
 714:           len--;
 715:     }
 716:     if (initShiftWords + initShiftBits > 0)
 717:       {
 718:         if (initShiftBits > 0)
 719:           {
 720:             int sh_out = MPN.lshift (x, initShiftWords, x, len, initShiftBits);
 721:             if (sh_out != 0)
 722:               x[(len++)+initShiftWords] = sh_out;
 723:           }
 724:         else
 725:           {
 726:             for (i = len; --i >= 0;)
 727:               x[i+initShiftWords] = x[i];
 728:           }
 729:         for (i = initShiftWords;  --i >= 0; )
 730:           x[i] = 0;
 731:         len += initShiftWords;
 732:       }
 733:     return len;
 734:   }
 735: 
 736:   public static int intLength (int i)
 737:   {
 738:     return 32 - count_leading_zeros (i < 0 ? ~i : i);
 739:   }
 740: 
 741:   /** Calcaulte the Common Lisp "integer-length" function.
 742:    * Assumes input is canonicalized:  len==BigInteger.wordsNeeded(words,len) */
 743:   public static int intLength (int[] words, int len)
 744:   {
 745:     len--;
 746:     return intLength (words[len]) + 32 * len;
 747:   }
 748: 
 749:   /* DEBUGGING:
 750:   public static void dprint (BigInteger x)
 751:   {
 752:     if (x.words == null)
 753:       System.err.print(Long.toString((long) x.ival & 0xffffffffL, 16));
 754:     else
 755:       dprint (System.err, x.words, x.ival);
 756:   }
 757:   public static void dprint (int[] x) { dprint (System.err, x, x.length); }
 758:   public static void dprint (int[] x, int len) { dprint (System.err, x, len); }
 759:   public static void dprint (java.io.PrintStream ps, int[] x, int len)
 760:   {
 761:     ps.print('(');
 762:     for (int i = 0;  i < len; i++)
 763:       {
 764:         if (i > 0)
 765:           ps.print (' ');
 766:         ps.print ("#x" + Long.toString ((long) x[i] & 0xffffffffL, 16));
 767:       }
 768:     ps.print(')');
 769:   }
 770:   */
 771: }