Robust/src/ClassLibrary/MGC/gnu/MPN.java

   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] += (int)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 }