Save This Page
Home » openjdk-7 » java » math » [javadoc | source]
    1   /*
    2    * Copyright (c) 1999, 2007, Oracle and/or its affiliates. All rights reserved.
    3    * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
    4    *
    5    * This code is free software; you can redistribute it and/or modify it
    6    * under the terms of the GNU General Public License version 2 only, as
    7    * published by the Free Software Foundation.  Oracle designates this
    8    * particular file as subject to the "Classpath" exception as provided
    9    * by Oracle in the LICENSE file that accompanied this code.
   10    *
   11    * This code is distributed in the hope that it will be useful, but WITHOUT
   12    * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
   13    * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
   14    * version 2 for more details (a copy is included in the LICENSE file that
   15    * accompanied this code).
   16    *
   17    * You should have received a copy of the GNU General Public License version
   18    * 2 along with this work; if not, write to the Free Software Foundation,
   19    * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
   20    *
   21    * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA
   22    * or visit www.oracle.com if you need additional information or have any
   23    * questions.
   24    */
   25   
   26   package java.math;
   27   
   28   /**
   29    * A class used to represent multiprecision integers that makes efficient
   30    * use of allocated space by allowing a number to occupy only part of
   31    * an array so that the arrays do not have to be reallocated as often.
   32    * When performing an operation with many iterations the array used to
   33    * hold a number is only reallocated when necessary and does not have to
   34    * be the same size as the number it represents. A mutable number allows
   35    * calculations to occur on the same number without having to create
   36    * a new number for every step of the calculation as occurs with
   37    * BigIntegers.
   38    *
   39    * @see     BigInteger
   40    * @author  Michael McCloskey
   41    * @since   1.3
   42    */
   43   
   44   import java.util.Arrays;
   45   
   46   import static java.math.BigInteger.LONG_MASK;
   47   import static java.math.BigDecimal.INFLATED;
   48   
   49   class MutableBigInteger {
   50       /**
   51        * Holds the magnitude of this MutableBigInteger in big endian order.
   52        * The magnitude may start at an offset into the value array, and it may
   53        * end before the length of the value array.
   54        */
   55       int[] value;
   56   
   57       /**
   58        * The number of ints of the value array that are currently used
   59        * to hold the magnitude of this MutableBigInteger. The magnitude starts
   60        * at an offset and offset + intLen may be less than value.length.
   61        */
   62       int intLen;
   63   
   64       /**
   65        * The offset into the value array where the magnitude of this
   66        * MutableBigInteger begins.
   67        */
   68       int offset = 0;
   69   
   70       // Constants
   71       /**
   72        * MutableBigInteger with one element value array with the value 1. Used by
   73        * BigDecimal divideAndRound to increment the quotient. Use this constant
   74        * only when the method is not going to modify this object.
   75        */
   76       static final MutableBigInteger ONE = new MutableBigInteger(1);
   77   
   78       // Constructors
   79   
   80       /**
   81        * The default constructor. An empty MutableBigInteger is created with
   82        * a one word capacity.
   83        */
   84       MutableBigInteger() {
   85           value = new int[1];
   86           intLen = 0;
   87       }
   88   
   89       /**
   90        * Construct a new MutableBigInteger with a magnitude specified by
   91        * the int val.
   92        */
   93       MutableBigInteger(int val) {
   94           value = new int[1];
   95           intLen = 1;
   96           value[0] = val;
   97       }
   98   
   99       /**
  100        * Construct a new MutableBigInteger with the specified value array
  101        * up to the length of the array supplied.
  102        */
  103       MutableBigInteger(int[] val) {
  104           value = val;
  105           intLen = val.length;
  106       }
  107   
  108       /**
  109        * Construct a new MutableBigInteger with a magnitude equal to the
  110        * specified BigInteger.
  111        */
  112       MutableBigInteger(BigInteger b) {
  113           intLen = b.mag.length;
  114           value = Arrays.copyOf(b.mag, intLen);
  115       }
  116   
  117       /**
  118        * Construct a new MutableBigInteger with a magnitude equal to the
  119        * specified MutableBigInteger.
  120        */
  121       MutableBigInteger(MutableBigInteger val) {
  122           intLen = val.intLen;
  123           value = Arrays.copyOfRange(val.value, val.offset, val.offset + intLen);
  124       }
  125   
  126       /**
  127        * Internal helper method to return the magnitude array. The caller is not
  128        * supposed to modify the returned array.
  129        */
  130       private int[] getMagnitudeArray() {
  131           if (offset > 0 || value.length != intLen)
  132               return Arrays.copyOfRange(value, offset, offset + intLen);
  133           return value;
  134       }
  135   
  136       /**
  137        * Convert this MutableBigInteger to a long value. The caller has to make
  138        * sure this MutableBigInteger can be fit into long.
  139        */
  140       private long toLong() {
  141           assert (intLen <= 2) : "this MutableBigInteger exceeds the range of long";
  142           if (intLen == 0)
  143               return 0;
  144           long d = value[offset] & LONG_MASK;
  145           return (intLen == 2) ? d << 32 | (value[offset + 1] & LONG_MASK) : d;
  146       }
  147   
  148       /**
  149        * Convert this MutableBigInteger to a BigInteger object.
  150        */
  151       BigInteger toBigInteger(int sign) {
  152           if (intLen == 0 || sign == 0)
  153               return BigInteger.ZERO;
  154           return new BigInteger(getMagnitudeArray(), sign);
  155       }
  156   
  157       /**
  158        * Convert this MutableBigInteger to BigDecimal object with the specified sign
  159        * and scale.
  160        */
  161       BigDecimal toBigDecimal(int sign, int scale) {
  162           if (intLen == 0 || sign == 0)
  163               return BigDecimal.valueOf(0, scale);
  164           int[] mag = getMagnitudeArray();
  165           int len = mag.length;
  166           int d = mag[0];
  167           // If this MutableBigInteger can't be fit into long, we need to
  168           // make a BigInteger object for the resultant BigDecimal object.
  169           if (len > 2 || (d < 0 && len == 2))
  170               return new BigDecimal(new BigInteger(mag, sign), INFLATED, scale, 0);
  171           long v = (len == 2) ?
  172               ((mag[1] & LONG_MASK) | (d & LONG_MASK) << 32) :
  173               d & LONG_MASK;
  174           return new BigDecimal(null, sign == -1 ? -v : v, scale, 0);
  175       }
  176   
  177       /**
  178        * Clear out a MutableBigInteger for reuse.
  179        */
  180       void clear() {
  181           offset = intLen = 0;
  182           for (int index=0, n=value.length; index < n; index++)
  183               value[index] = 0;
  184       }
  185   
  186       /**
  187        * Set a MutableBigInteger to zero, removing its offset.
  188        */
  189       void reset() {
  190           offset = intLen = 0;
  191       }
  192   
  193       /**
  194        * Compare the magnitude of two MutableBigIntegers. Returns -1, 0 or 1
  195        * as this MutableBigInteger is numerically less than, equal to, or
  196        * greater than <tt>b</tt>.
  197        */
  198       final int compare(MutableBigInteger b) {
  199           int blen = b.intLen;
  200           if (intLen < blen)
  201               return -1;
  202           if (intLen > blen)
  203              return 1;
  204   
  205           // Add Integer.MIN_VALUE to make the comparison act as unsigned integer
  206           // comparison.
  207           int[] bval = b.value;
  208           for (int i = offset, j = b.offset; i < intLen + offset; i++, j++) {
  209               int b1 = value[i] + 0x80000000;
  210               int b2 = bval[j]  + 0x80000000;
  211               if (b1 < b2)
  212                   return -1;
  213               if (b1 > b2)
  214                   return 1;
  215           }
  216           return 0;
  217       }
  218   
  219       /**
  220        * Compare this against half of a MutableBigInteger object (Needed for
  221        * remainder tests).
  222        * Assumes no leading unnecessary zeros, which holds for results
  223        * from divide().
  224        */
  225       final int compareHalf(MutableBigInteger b) {
  226           int blen = b.intLen;
  227           int len = intLen;
  228           if (len <= 0)
  229               return blen <=0 ? 0 : -1;
  230           if (len > blen)
  231               return 1;
  232           if (len < blen - 1)
  233               return -1;
  234           int[] bval = b.value;
  235           int bstart = 0;
  236           int carry = 0;
  237           // Only 2 cases left:len == blen or len == blen - 1
  238           if (len != blen) { // len == blen - 1
  239               if (bval[bstart] == 1) {
  240                   ++bstart;
  241                   carry = 0x80000000;
  242               } else
  243                   return -1;
  244           }
  245           // compare values with right-shifted values of b,
  246           // carrying shifted-out bits across words
  247           int[] val = value;
  248           for (int i = offset, j = bstart; i < len + offset;) {
  249               int bv = bval[j++];
  250               long hb = ((bv >>> 1) + carry) & LONG_MASK;
  251               long v = val[i++] & LONG_MASK;
  252               if (v != hb)
  253                   return v < hb ? -1 : 1;
  254               carry = (bv & 1) << 31; // carray will be either 0x80000000 or 0
  255           }
  256           return carry == 0? 0 : -1;
  257       }
  258   
  259       /**
  260        * Return the index of the lowest set bit in this MutableBigInteger. If the
  261        * magnitude of this MutableBigInteger is zero, -1 is returned.
  262        */
  263       private final int getLowestSetBit() {
  264           if (intLen == 0)
  265               return -1;
  266           int j, b;
  267           for (j=intLen-1; (j>0) && (value[j+offset]==0); j--)
  268               ;
  269           b = value[j+offset];
  270           if (b==0)
  271               return -1;
  272           return ((intLen-1-j)<<5) + Integer.numberOfTrailingZeros(b);
  273       }
  274   
  275       /**
  276        * Return the int in use in this MutableBigInteger at the specified
  277        * index. This method is not used because it is not inlined on all
  278        * platforms.
  279        */
  280       private final int getInt(int index) {
  281           return value[offset+index];
  282       }
  283   
  284       /**
  285        * Return a long which is equal to the unsigned value of the int in
  286        * use in this MutableBigInteger at the specified index. This method is
  287        * not used because it is not inlined on all platforms.
  288        */
  289       private final long getLong(int index) {
  290           return value[offset+index] & LONG_MASK;
  291       }
  292   
  293       /**
  294        * Ensure that the MutableBigInteger is in normal form, specifically
  295        * making sure that there are no leading zeros, and that if the
  296        * magnitude is zero, then intLen is zero.
  297        */
  298       final void normalize() {
  299           if (intLen == 0) {
  300               offset = 0;
  301               return;
  302           }
  303   
  304           int index = offset;
  305           if (value[index] != 0)
  306               return;
  307   
  308           int indexBound = index+intLen;
  309           do {
  310               index++;
  311           } while(index < indexBound && value[index]==0);
  312   
  313           int numZeros = index - offset;
  314           intLen -= numZeros;
  315           offset = (intLen==0 ?  0 : offset+numZeros);
  316       }
  317   
  318       /**
  319        * If this MutableBigInteger cannot hold len words, increase the size
  320        * of the value array to len words.
  321        */
  322       private final void ensureCapacity(int len) {
  323           if (value.length < len) {
  324               value = new int[len];
  325               offset = 0;
  326               intLen = len;
  327           }
  328       }
  329   
  330       /**
  331        * Convert this MutableBigInteger into an int array with no leading
  332        * zeros, of a length that is equal to this MutableBigInteger's intLen.
  333        */
  334       int[] toIntArray() {
  335           int[] result = new int[intLen];
  336           for(int i=0; i<intLen; i++)
  337               result[i] = value[offset+i];
  338           return result;
  339       }
  340   
  341       /**
  342        * Sets the int at index+offset in this MutableBigInteger to val.
  343        * This does not get inlined on all platforms so it is not used
  344        * as often as originally intended.
  345        */
  346       void setInt(int index, int val) {
  347           value[offset + index] = val;
  348       }
  349   
  350       /**
  351        * Sets this MutableBigInteger's value array to the specified array.
  352        * The intLen is set to the specified length.
  353        */
  354       void setValue(int[] val, int length) {
  355           value = val;
  356           intLen = length;
  357           offset = 0;
  358       }
  359   
  360       /**
  361        * Sets this MutableBigInteger's value array to a copy of the specified
  362        * array. The intLen is set to the length of the new array.
  363        */
  364       void copyValue(MutableBigInteger src) {
  365           int len = src.intLen;
  366           if (value.length < len)
  367               value = new int[len];
  368           System.arraycopy(src.value, src.offset, value, 0, len);
  369           intLen = len;
  370           offset = 0;
  371       }
  372   
  373       /**
  374        * Sets this MutableBigInteger's value array to a copy of the specified
  375        * array. The intLen is set to the length of the specified array.
  376        */
  377       void copyValue(int[] val) {
  378           int len = val.length;
  379           if (value.length < len)
  380               value = new int[len];
  381           System.arraycopy(val, 0, value, 0, len);
  382           intLen = len;
  383           offset = 0;
  384       }
  385   
  386       /**
  387        * Returns true iff this MutableBigInteger has a value of one.
  388        */
  389       boolean isOne() {
  390           return (intLen == 1) && (value[offset] == 1);
  391       }
  392   
  393       /**
  394        * Returns true iff this MutableBigInteger has a value of zero.
  395        */
  396       boolean isZero() {
  397           return (intLen == 0);
  398       }
  399   
  400       /**
  401        * Returns true iff this MutableBigInteger is even.
  402        */
  403       boolean isEven() {
  404           return (intLen == 0) || ((value[offset + intLen - 1] & 1) == 0);
  405       }
  406   
  407       /**
  408        * Returns true iff this MutableBigInteger is odd.
  409        */
  410       boolean isOdd() {
  411           return isZero() ? false : ((value[offset + intLen - 1] & 1) == 1);
  412       }
  413   
  414       /**
  415        * Returns true iff this MutableBigInteger is in normal form. A
  416        * MutableBigInteger is in normal form if it has no leading zeros
  417        * after the offset, and intLen + offset <= value.length.
  418        */
  419       boolean isNormal() {
  420           if (intLen + offset > value.length)
  421               return false;
  422           if (intLen ==0)
  423               return true;
  424           return (value[offset] != 0);
  425       }
  426   
  427       /**
  428        * Returns a String representation of this MutableBigInteger in radix 10.
  429        */
  430       public String toString() {
  431           BigInteger b = toBigInteger(1);
  432           return b.toString();
  433       }
  434   
  435       /**
  436        * Right shift this MutableBigInteger n bits. The MutableBigInteger is left
  437        * in normal form.
  438        */
  439       void rightShift(int n) {
  440           if (intLen == 0)
  441               return;
  442           int nInts = n >>> 5;
  443           int nBits = n & 0x1F;
  444           this.intLen -= nInts;
  445           if (nBits == 0)
  446               return;
  447           int bitsInHighWord = BigInteger.bitLengthForInt(value[offset]);
  448           if (nBits >= bitsInHighWord) {
  449               this.primitiveLeftShift(32 - nBits);
  450               this.intLen--;
  451           } else {
  452               primitiveRightShift(nBits);
  453           }
  454       }
  455   
  456       /**
  457        * Left shift this MutableBigInteger n bits.
  458        */
  459       void leftShift(int n) {
  460           /*
  461            * If there is enough storage space in this MutableBigInteger already
  462            * the available space will be used. Space to the right of the used
  463            * ints in the value array is faster to utilize, so the extra space
  464            * will be taken from the right if possible.
  465            */
  466           if (intLen == 0)
  467              return;
  468           int nInts = n >>> 5;
  469           int nBits = n&0x1F;
  470           int bitsInHighWord = BigInteger.bitLengthForInt(value[offset]);
  471   
  472           // If shift can be done without moving words, do so
  473           if (n <= (32-bitsInHighWord)) {
  474               primitiveLeftShift(nBits);
  475               return;
  476           }
  477   
  478           int newLen = intLen + nInts +1;
  479           if (nBits <= (32-bitsInHighWord))
  480               newLen--;
  481           if (value.length < newLen) {
  482               // The array must grow
  483               int[] result = new int[newLen];
  484               for (int i=0; i<intLen; i++)
  485                   result[i] = value[offset+i];
  486               setValue(result, newLen);
  487           } else if (value.length - offset >= newLen) {
  488               // Use space on right
  489               for(int i=0; i<newLen - intLen; i++)
  490                   value[offset+intLen+i] = 0;
  491           } else {
  492               // Must use space on left
  493               for (int i=0; i<intLen; i++)
  494                   value[i] = value[offset+i];
  495               for (int i=intLen; i<newLen; i++)
  496                   value[i] = 0;
  497               offset = 0;
  498           }
  499           intLen = newLen;
  500           if (nBits == 0)
  501               return;
  502           if (nBits <= (32-bitsInHighWord))
  503               primitiveLeftShift(nBits);
  504           else
  505               primitiveRightShift(32 -nBits);
  506       }
  507   
  508       /**
  509        * A primitive used for division. This method adds in one multiple of the
  510        * divisor a back to the dividend result at a specified offset. It is used
  511        * when qhat was estimated too large, and must be adjusted.
  512        */
  513       private int divadd(int[] a, int[] result, int offset) {
  514           long carry = 0;
  515   
  516           for (int j=a.length-1; j >= 0; j--) {
  517               long sum = (a[j] & LONG_MASK) +
  518                          (result[j+offset] & LONG_MASK) + carry;
  519               result[j+offset] = (int)sum;
  520               carry = sum >>> 32;
  521           }
  522           return (int)carry;
  523       }
  524   
  525       /**
  526        * This method is used for division. It multiplies an n word input a by one
  527        * word input x, and subtracts the n word product from q. This is needed
  528        * when subtracting qhat*divisor from dividend.
  529        */
  530       private int mulsub(int[] q, int[] a, int x, int len, int offset) {
  531           long xLong = x & LONG_MASK;
  532           long carry = 0;
  533           offset += len;
  534   
  535           for (int j=len-1; j >= 0; j--) {
  536               long product = (a[j] & LONG_MASK) * xLong + carry;
  537               long difference = q[offset] - product;
  538               q[offset--] = (int)difference;
  539               carry = (product >>> 32)
  540                        + (((difference & LONG_MASK) >
  541                            (((~(int)product) & LONG_MASK))) ? 1:0);
  542           }
  543           return (int)carry;
  544       }
  545   
  546       /**
  547        * Right shift this MutableBigInteger n bits, where n is
  548        * less than 32.
  549        * Assumes that intLen > 0, n > 0 for speed
  550        */
  551       private final void primitiveRightShift(int n) {
  552           int[] val = value;
  553           int n2 = 32 - n;
  554           for (int i=offset+intLen-1, c=val[i]; i>offset; i--) {
  555               int b = c;
  556               c = val[i-1];
  557               val[i] = (c << n2) | (b >>> n);
  558           }
  559           val[offset] >>>= n;
  560       }
  561   
  562       /**
  563        * Left shift this MutableBigInteger n bits, where n is
  564        * less than 32.
  565        * Assumes that intLen > 0, n > 0 for speed
  566        */
  567       private final void primitiveLeftShift(int n) {
  568           int[] val = value;
  569           int n2 = 32 - n;
  570           for (int i=offset, c=val[i], m=i+intLen-1; i<m; i++) {
  571               int b = c;
  572               c = val[i+1];
  573               val[i] = (b << n) | (c >>> n2);
  574           }
  575           val[offset+intLen-1] <<= n;
  576       }
  577   
  578       /**
  579        * Adds the contents of two MutableBigInteger objects.The result
  580        * is placed within this MutableBigInteger.
  581        * The contents of the addend are not changed.
  582        */
  583       void add(MutableBigInteger addend) {
  584           int x = intLen;
  585           int y = addend.intLen;
  586           int resultLen = (intLen > addend.intLen ? intLen : addend.intLen);
  587           int[] result = (value.length < resultLen ? new int[resultLen] : value);
  588   
  589           int rstart = result.length-1;
  590           long sum;
  591           long carry = 0;
  592   
  593           // Add common parts of both numbers
  594           while(x>0 && y>0) {
  595               x--; y--;
  596               sum = (value[x+offset] & LONG_MASK) +
  597                   (addend.value[y+addend.offset] & LONG_MASK) + carry;
  598               result[rstart--] = (int)sum;
  599               carry = sum >>> 32;
  600           }
  601   
  602           // Add remainder of the longer number
  603           while(x>0) {
  604               x--;
  605               if (carry == 0 && result == value && rstart == (x + offset))
  606                   return;
  607               sum = (value[x+offset] & LONG_MASK) + carry;
  608               result[rstart--] = (int)sum;
  609               carry = sum >>> 32;
  610           }
  611           while(y>0) {
  612               y--;
  613               sum = (addend.value[y+addend.offset] & LONG_MASK) + carry;
  614               result[rstart--] = (int)sum;
  615               carry = sum >>> 32;
  616           }
  617   
  618           if (carry > 0) { // Result must grow in length
  619               resultLen++;
  620               if (result.length < resultLen) {
  621                   int temp[] = new int[resultLen];
  622                   // Result one word longer from carry-out; copy low-order
  623                   // bits into new result.
  624                   System.arraycopy(result, 0, temp, 1, result.length);
  625                   temp[0] = 1;
  626                   result = temp;
  627               } else {
  628                   result[rstart--] = 1;
  629               }
  630           }
  631   
  632           value = result;
  633           intLen = resultLen;
  634           offset = result.length - resultLen;
  635       }
  636   
  637   
  638       /**
  639        * Subtracts the smaller of this and b from the larger and places the
  640        * result into this MutableBigInteger.
  641        */
  642       int subtract(MutableBigInteger b) {
  643           MutableBigInteger a = this;
  644   
  645           int[] result = value;
  646           int sign = a.compare(b);
  647   
  648           if (sign == 0) {
  649               reset();
  650               return 0;
  651           }
  652           if (sign < 0) {
  653               MutableBigInteger tmp = a;
  654               a = b;
  655               b = tmp;
  656           }
  657   
  658           int resultLen = a.intLen;
  659           if (result.length < resultLen)
  660               result = new int[resultLen];
  661   
  662           long diff = 0;
  663           int x = a.intLen;
  664           int y = b.intLen;
  665           int rstart = result.length - 1;
  666   
  667           // Subtract common parts of both numbers
  668           while (y>0) {
  669               x--; y--;
  670   
  671               diff = (a.value[x+a.offset] & LONG_MASK) -
  672                      (b.value[y+b.offset] & LONG_MASK) - ((int)-(diff>>32));
  673               result[rstart--] = (int)diff;
  674           }
  675           // Subtract remainder of longer number
  676           while (x>0) {
  677               x--;
  678               diff = (a.value[x+a.offset] & LONG_MASK) - ((int)-(diff>>32));
  679               result[rstart--] = (int)diff;
  680           }
  681   
  682           value = result;
  683           intLen = resultLen;
  684           offset = value.length - resultLen;
  685           normalize();
  686           return sign;
  687       }
  688   
  689       /**
  690        * Subtracts the smaller of a and b from the larger and places the result
  691        * into the larger. Returns 1 if the answer is in a, -1 if in b, 0 if no
  692        * operation was performed.
  693        */
  694       private int difference(MutableBigInteger b) {
  695           MutableBigInteger a = this;
  696           int sign = a.compare(b);
  697           if (sign ==0)
  698               return 0;
  699           if (sign < 0) {
  700               MutableBigInteger tmp = a;
  701               a = b;
  702               b = tmp;
  703           }
  704   
  705           long diff = 0;
  706           int x = a.intLen;
  707           int y = b.intLen;
  708   
  709           // Subtract common parts of both numbers
  710           while (y>0) {
  711               x--; y--;
  712               diff = (a.value[a.offset+ x] & LONG_MASK) -
  713                   (b.value[b.offset+ y] & LONG_MASK) - ((int)-(diff>>32));
  714               a.value[a.offset+x] = (int)diff;
  715           }
  716           // Subtract remainder of longer number
  717           while (x>0) {
  718               x--;
  719               diff = (a.value[a.offset+ x] & LONG_MASK) - ((int)-(diff>>32));
  720               a.value[a.offset+x] = (int)diff;
  721           }
  722   
  723           a.normalize();
  724           return sign;
  725       }
  726   
  727       /**
  728        * Multiply the contents of two MutableBigInteger objects. The result is
  729        * placed into MutableBigInteger z. The contents of y are not changed.
  730        */
  731       void multiply(MutableBigInteger y, MutableBigInteger z) {
  732           int xLen = intLen;
  733           int yLen = y.intLen;
  734           int newLen = xLen + yLen;
  735   
  736           // Put z into an appropriate state to receive product
  737           if (z.value.length < newLen)
  738               z.value = new int[newLen];
  739           z.offset = 0;
  740           z.intLen = newLen;
  741   
  742           // The first iteration is hoisted out of the loop to avoid extra add
  743           long carry = 0;
  744           for (int j=yLen-1, k=yLen+xLen-1; j >= 0; j--, k--) {
  745                   long product = (y.value[j+y.offset] & LONG_MASK) *
  746                                  (value[xLen-1+offset] & LONG_MASK) + carry;
  747                   z.value[k] = (int)product;
  748                   carry = product >>> 32;
  749           }
  750           z.value[xLen-1] = (int)carry;
  751   
  752           // Perform the multiplication word by word
  753           for (int i = xLen-2; i >= 0; i--) {
  754               carry = 0;
  755               for (int j=yLen-1, k=yLen+i; j >= 0; j--, k--) {
  756                   long product = (y.value[j+y.offset] & LONG_MASK) *
  757                                  (value[i+offset] & LONG_MASK) +
  758                                  (z.value[k] & LONG_MASK) + carry;
  759                   z.value[k] = (int)product;
  760                   carry = product >>> 32;
  761               }
  762               z.value[i] = (int)carry;
  763           }
  764   
  765           // Remove leading zeros from product
  766           z.normalize();
  767       }
  768   
  769       /**
  770        * Multiply the contents of this MutableBigInteger by the word y. The
  771        * result is placed into z.
  772        */
  773       void mul(int y, MutableBigInteger z) {
  774           if (y == 1) {
  775               z.copyValue(this);
  776               return;
  777           }
  778   
  779           if (y == 0) {
  780               z.clear();
  781               return;
  782           }
  783   
  784           // Perform the multiplication word by word
  785           long ylong = y & LONG_MASK;
  786           int[] zval = (z.value.length<intLen+1 ? new int[intLen + 1]
  787                                                 : z.value);
  788           long carry = 0;
  789           for (int i = intLen-1; i >= 0; i--) {
  790               long product = ylong * (value[i+offset] & LONG_MASK) + carry;
  791               zval[i+1] = (int)product;
  792               carry = product >>> 32;
  793           }
  794   
  795           if (carry == 0) {
  796               z.offset = 1;
  797               z.intLen = intLen;
  798           } else {
  799               z.offset = 0;
  800               z.intLen = intLen + 1;
  801               zval[0] = (int)carry;
  802           }
  803           z.value = zval;
  804       }
  805   
  806        /**
  807        * This method is used for division of an n word dividend by a one word
  808        * divisor. The quotient is placed into quotient. The one word divisor is
  809        * specified by divisor.
  810        *
  811        * @return the remainder of the division is returned.
  812        *
  813        */
  814       int divideOneWord(int divisor, MutableBigInteger quotient) {
  815           long divisorLong = divisor & LONG_MASK;
  816   
  817           // Special case of one word dividend
  818           if (intLen == 1) {
  819               long dividendValue = value[offset] & LONG_MASK;
  820               int q = (int) (dividendValue / divisorLong);
  821               int r = (int) (dividendValue - q * divisorLong);
  822               quotient.value[0] = q;
  823               quotient.intLen = (q == 0) ? 0 : 1;
  824               quotient.offset = 0;
  825               return r;
  826           }
  827   
  828           if (quotient.value.length < intLen)
  829               quotient.value = new int[intLen];
  830           quotient.offset = 0;
  831           quotient.intLen = intLen;
  832   
  833           // Normalize the divisor
  834           int shift = Integer.numberOfLeadingZeros(divisor);
  835   
  836           int rem = value[offset];
  837           long remLong = rem & LONG_MASK;
  838           if (remLong < divisorLong) {
  839               quotient.value[0] = 0;
  840           } else {
  841               quotient.value[0] = (int)(remLong / divisorLong);
  842               rem = (int) (remLong - (quotient.value[0] * divisorLong));
  843               remLong = rem & LONG_MASK;
  844           }
  845   
  846           int xlen = intLen;
  847           int[] qWord = new int[2];
  848           while (--xlen > 0) {
  849               long dividendEstimate = (remLong<<32) |
  850                   (value[offset + intLen - xlen] & LONG_MASK);
  851               if (dividendEstimate >= 0) {
  852                   qWord[0] = (int) (dividendEstimate / divisorLong);
  853                   qWord[1] = (int) (dividendEstimate - qWord[0] * divisorLong);
  854               } else {
  855                   divWord(qWord, dividendEstimate, divisor);
  856               }
  857               quotient.value[intLen - xlen] = qWord[0];
  858               rem = qWord[1];
  859               remLong = rem & LONG_MASK;
  860           }
  861   
  862           quotient.normalize();
  863           // Unnormalize
  864           if (shift > 0)
  865               return rem % divisor;
  866           else
  867               return rem;
  868       }
  869   
  870       /**
  871        * Calculates the quotient of this div b and places the quotient in the
  872        * provided MutableBigInteger objects and the remainder object is returned.
  873        *
  874        * Uses Algorithm D in Knuth section 4.3.1.
  875        * Many optimizations to that algorithm have been adapted from the Colin
  876        * Plumb C library.
  877        * It special cases one word divisors for speed. The content of b is not
  878        * changed.
  879        *
  880        */
  881       MutableBigInteger divide(MutableBigInteger b, MutableBigInteger quotient) {
  882           if (b.intLen == 0)
  883               throw new ArithmeticException("BigInteger divide by zero");
  884   
  885           // Dividend is zero
  886           if (intLen == 0) {
  887               quotient.intLen = quotient.offset;
  888               return new MutableBigInteger();
  889           }
  890   
  891           int cmp = compare(b);
  892           // Dividend less than divisor
  893           if (cmp < 0) {
  894               quotient.intLen = quotient.offset = 0;
  895               return new MutableBigInteger(this);
  896           }
  897           // Dividend equal to divisor
  898           if (cmp == 0) {
  899               quotient.value[0] = quotient.intLen = 1;
  900               quotient.offset = 0;
  901               return new MutableBigInteger();
  902           }
  903   
  904           quotient.clear();
  905           // Special case one word divisor
  906           if (b.intLen == 1) {
  907               int r = divideOneWord(b.value[b.offset], quotient);
  908               if (r == 0)
  909                   return new MutableBigInteger();
  910               return new MutableBigInteger(r);
  911           }
  912   
  913           // Copy divisor value to protect divisor
  914           int[] div = Arrays.copyOfRange(b.value, b.offset, b.offset + b.intLen);
  915           return divideMagnitude(div, quotient);
  916       }
  917   
  918       /**
  919        * Internally used  to calculate the quotient of this div v and places the
  920        * quotient in the provided MutableBigInteger object and the remainder is
  921        * returned.
  922        *
  923        * @return the remainder of the division will be returned.
  924        */
  925       long divide(long v, MutableBigInteger quotient) {
  926           if (v == 0)
  927               throw new ArithmeticException("BigInteger divide by zero");
  928   
  929           // Dividend is zero
  930           if (intLen == 0) {
  931               quotient.intLen = quotient.offset = 0;
  932               return 0;
  933           }
  934           if (v < 0)
  935               v = -v;
  936   
  937           int d = (int)(v >>> 32);
  938           quotient.clear();
  939           // Special case on word divisor
  940           if (d == 0)
  941               return divideOneWord((int)v, quotient) & LONG_MASK;
  942           else {
  943               int[] div = new int[]{ d, (int)(v & LONG_MASK) };
  944               return divideMagnitude(div, quotient).toLong();
  945           }
  946       }
  947   
  948       /**
  949        * Divide this MutableBigInteger by the divisor represented by its magnitude
  950        * array. The quotient will be placed into the provided quotient object &
  951        * the remainder object is returned.
  952        */
  953       private MutableBigInteger divideMagnitude(int[] divisor,
  954                                                 MutableBigInteger quotient) {
  955   
  956           // Remainder starts as dividend with space for a leading zero
  957           MutableBigInteger rem = new MutableBigInteger(new int[intLen + 1]);
  958           System.arraycopy(value, offset, rem.value, 1, intLen);
  959           rem.intLen = intLen;
  960           rem.offset = 1;
  961   
  962           int nlen = rem.intLen;
  963   
  964           // Set the quotient size
  965           int dlen = divisor.length;
  966           int limit = nlen - dlen + 1;
  967           if (quotient.value.length < limit) {
  968               quotient.value = new int[limit];
  969               quotient.offset = 0;
  970           }
  971           quotient.intLen = limit;
  972           int[] q = quotient.value;
  973   
  974           // D1 normalize the divisor
  975           int shift = Integer.numberOfLeadingZeros(divisor[0]);
  976           if (shift > 0) {
  977               // First shift will not grow array
  978               BigInteger.primitiveLeftShift(divisor, dlen, shift);
  979               // But this one might
  980               rem.leftShift(shift);
  981           }
  982   
  983           // Must insert leading 0 in rem if its length did not change
  984           if (rem.intLen == nlen) {
  985               rem.offset = 0;
  986               rem.value[0] = 0;
  987               rem.intLen++;
  988           }
  989   
  990           int dh = divisor[0];
  991           long dhLong = dh & LONG_MASK;
  992           int dl = divisor[1];
  993           int[] qWord = new int[2];
  994   
  995           // D2 Initialize j
  996           for(int j=0; j<limit; j++) {
  997               // D3 Calculate qhat
  998               // estimate qhat
  999               int qhat = 0;
 1000               int qrem = 0;
 1001               boolean skipCorrection = false;
 1002               int nh = rem.value[j+rem.offset];
 1003               int nh2 = nh + 0x80000000;
 1004               int nm = rem.value[j+1+rem.offset];
 1005   
 1006               if (nh == dh) {
 1007                   qhat = ~0;
 1008                   qrem = nh + nm;
 1009                   skipCorrection = qrem + 0x80000000 < nh2;
 1010               } else {
 1011                   long nChunk = (((long)nh) << 32) | (nm & LONG_MASK);
 1012                   if (nChunk >= 0) {
 1013                       qhat = (int) (nChunk / dhLong);
 1014                       qrem = (int) (nChunk - (qhat * dhLong));
 1015                   } else {
 1016                       divWord(qWord, nChunk, dh);
 1017                       qhat = qWord[0];
 1018                       qrem = qWord[1];
 1019                   }
 1020               }
 1021   
 1022               if (qhat == 0)
 1023                   continue;
 1024   
 1025               if (!skipCorrection) { // Correct qhat
 1026                   long nl = rem.value[j+2+rem.offset] & LONG_MASK;
 1027                   long rs = ((qrem & LONG_MASK) << 32) | nl;
 1028                   long estProduct = (dl & LONG_MASK) * (qhat & LONG_MASK);
 1029   
 1030                   if (unsignedLongCompare(estProduct, rs)) {
 1031                       qhat--;
 1032                       qrem = (int)((qrem & LONG_MASK) + dhLong);
 1033                       if ((qrem & LONG_MASK) >=  dhLong) {
 1034                           estProduct -= (dl & LONG_MASK);
 1035                           rs = ((qrem & LONG_MASK) << 32) | nl;
 1036                           if (unsignedLongCompare(estProduct, rs))
 1037                               qhat--;
 1038                       }
 1039                   }
 1040               }
 1041   
 1042               // D4 Multiply and subtract
 1043               rem.value[j+rem.offset] = 0;
 1044               int borrow = mulsub(rem.value, divisor, qhat, dlen, j+rem.offset);
 1045   
 1046               // D5 Test remainder
 1047               if (borrow + 0x80000000 > nh2) {
 1048                   // D6 Add back
 1049                   divadd(divisor, rem.value, j+1+rem.offset);
 1050                   qhat--;
 1051               }
 1052   
 1053               // Store the quotient digit
 1054               q[j] = qhat;
 1055           } // D7 loop on j
 1056   
 1057           // D8 Unnormalize
 1058           if (shift > 0)
 1059               rem.rightShift(shift);
 1060   
 1061           quotient.normalize();
 1062           rem.normalize();
 1063           return rem;
 1064       }
 1065   
 1066       /**
 1067        * Compare two longs as if they were unsigned.
 1068        * Returns true iff one is bigger than two.
 1069        */
 1070       private boolean unsignedLongCompare(long one, long two) {
 1071           return (one+Long.MIN_VALUE) > (two+Long.MIN_VALUE);
 1072       }
 1073   
 1074       /**
 1075        * This method divides a long quantity by an int to estimate
 1076        * qhat for two multi precision numbers. It is used when
 1077        * the signed value of n is less than zero.
 1078        */
 1079       private void divWord(int[] result, long n, int d) {
 1080           long dLong = d & LONG_MASK;
 1081   
 1082           if (dLong == 1) {
 1083               result[0] = (int)n;
 1084               result[1] = 0;
 1085               return;
 1086           }
 1087   
 1088           // Approximate the quotient and remainder
 1089           long q = (n >>> 1) / (dLong >>> 1);
 1090           long r = n - q*dLong;
 1091   
 1092           // Correct the approximation
 1093           while (r < 0) {
 1094               r += dLong;
 1095               q--;
 1096           }
 1097           while (r >= dLong) {
 1098               r -= dLong;
 1099               q++;
 1100           }
 1101   
 1102           // n - q*dlong == r && 0 <= r <dLong, hence we're done.
 1103           result[0] = (int)q;
 1104           result[1] = (int)r;
 1105       }
 1106   
 1107       /**
 1108        * Calculate GCD of this and b. This and b are changed by the computation.
 1109        */
 1110       MutableBigInteger hybridGCD(MutableBigInteger b) {
 1111           // Use Euclid's algorithm until the numbers are approximately the
 1112           // same length, then use the binary GCD algorithm to find the GCD.
 1113           MutableBigInteger a = this;
 1114           MutableBigInteger q = new MutableBigInteger();
 1115   
 1116           while (b.intLen != 0) {
 1117               if (Math.abs(a.intLen - b.intLen) < 2)
 1118                   return a.binaryGCD(b);
 1119   
 1120               MutableBigInteger r = a.divide(b, q);
 1121               a = b;
 1122               b = r;
 1123           }
 1124           return a;
 1125       }
 1126   
 1127       /**
 1128        * Calculate GCD of this and v.
 1129        * Assumes that this and v are not zero.
 1130        */
 1131       private MutableBigInteger binaryGCD(MutableBigInteger v) {
 1132           // Algorithm B from Knuth section 4.5.2
 1133           MutableBigInteger u = this;
 1134           MutableBigInteger r = new MutableBigInteger();
 1135   
 1136           // step B1
 1137           int s1 = u.getLowestSetBit();
 1138           int s2 = v.getLowestSetBit();
 1139           int k = (s1 < s2) ? s1 : s2;
 1140           if (k != 0) {
 1141               u.rightShift(k);
 1142               v.rightShift(k);
 1143           }
 1144   
 1145           // step B2
 1146           boolean uOdd = (k==s1);
 1147           MutableBigInteger t = uOdd ? v: u;
 1148           int tsign = uOdd ? -1 : 1;
 1149   
 1150           int lb;
 1151           while ((lb = t.getLowestSetBit()) >= 0) {
 1152               // steps B3 and B4
 1153               t.rightShift(lb);
 1154               // step B5
 1155               if (tsign > 0)
 1156                   u = t;
 1157               else
 1158                   v = t;
 1159   
 1160               // Special case one word numbers
 1161               if (u.intLen < 2 && v.intLen < 2) {
 1162                   int x = u.value[u.offset];
 1163                   int y = v.value[v.offset];
 1164                   x  = binaryGcd(x, y);
 1165                   r.value[0] = x;
 1166                   r.intLen = 1;
 1167                   r.offset = 0;
 1168                   if (k > 0)
 1169                       r.leftShift(k);
 1170                   return r;
 1171               }
 1172   
 1173               // step B6
 1174               if ((tsign = u.difference(v)) == 0)
 1175                   break;
 1176               t = (tsign >= 0) ? u : v;
 1177           }
 1178   
 1179           if (k > 0)
 1180               u.leftShift(k);
 1181           return u;
 1182       }
 1183   
 1184       /**
 1185        * Calculate GCD of a and b interpreted as unsigned integers.
 1186        */
 1187       static int binaryGcd(int a, int b) {
 1188           if (b==0)
 1189               return a;
 1190           if (a==0)
 1191               return b;
 1192   
 1193           // Right shift a & b till their last bits equal to 1.
 1194           int aZeros = Integer.numberOfTrailingZeros(a);
 1195           int bZeros = Integer.numberOfTrailingZeros(b);
 1196           a >>>= aZeros;
 1197           b >>>= bZeros;
 1198   
 1199           int t = (aZeros < bZeros ? aZeros : bZeros);
 1200   
 1201           while (a != b) {
 1202               if ((a+0x80000000) > (b+0x80000000)) {  // a > b as unsigned
 1203                   a -= b;
 1204                   a >>>= Integer.numberOfTrailingZeros(a);
 1205               } else {
 1206                   b -= a;
 1207                   b >>>= Integer.numberOfTrailingZeros(b);
 1208               }
 1209           }
 1210           return a<<t;
 1211       }
 1212   
 1213       /**
 1214        * Returns the modInverse of this mod p. This and p are not affected by
 1215        * the operation.
 1216        */
 1217       MutableBigInteger mutableModInverse(MutableBigInteger p) {
 1218           // Modulus is odd, use Schroeppel's algorithm
 1219           if (p.isOdd())
 1220               return modInverse(p);
 1221   
 1222           // Base and modulus are even, throw exception
 1223           if (isEven())
 1224               throw new ArithmeticException("BigInteger not invertible.");
 1225   
 1226           // Get even part of modulus expressed as a power of 2
 1227           int powersOf2 = p.getLowestSetBit();
 1228   
 1229           // Construct odd part of modulus
 1230           MutableBigInteger oddMod = new MutableBigInteger(p);
 1231           oddMod.rightShift(powersOf2);
 1232   
 1233           if (oddMod.isOne())
 1234               return modInverseMP2(powersOf2);
 1235   
 1236           // Calculate 1/a mod oddMod
 1237           MutableBigInteger oddPart = modInverse(oddMod);
 1238   
 1239           // Calculate 1/a mod evenMod
 1240           MutableBigInteger evenPart = modInverseMP2(powersOf2);
 1241   
 1242           // Combine the results using Chinese Remainder Theorem
 1243           MutableBigInteger y1 = modInverseBP2(oddMod, powersOf2);
 1244           MutableBigInteger y2 = oddMod.modInverseMP2(powersOf2);
 1245   
 1246           MutableBigInteger temp1 = new MutableBigInteger();
 1247           MutableBigInteger temp2 = new MutableBigInteger();
 1248           MutableBigInteger result = new MutableBigInteger();
 1249   
 1250           oddPart.leftShift(powersOf2);
 1251           oddPart.multiply(y1, result);
 1252   
 1253           evenPart.multiply(oddMod, temp1);
 1254           temp1.multiply(y2, temp2);
 1255   
 1256           result.add(temp2);
 1257           return result.divide(p, temp1);
 1258       }
 1259   
 1260       /*
 1261        * Calculate the multiplicative inverse of this mod 2^k.
 1262        */
 1263       MutableBigInteger modInverseMP2(int k) {
 1264           if (isEven())
 1265               throw new ArithmeticException("Non-invertible. (GCD != 1)");
 1266   
 1267           if (k > 64)
 1268               return euclidModInverse(k);
 1269   
 1270           int t = inverseMod32(value[offset+intLen-1]);
 1271   
 1272           if (k < 33) {
 1273               t = (k == 32 ? t : t & ((1 << k) - 1));
 1274               return new MutableBigInteger(t);
 1275           }
 1276   
 1277           long pLong = (value[offset+intLen-1] & LONG_MASK);
 1278           if (intLen > 1)
 1279               pLong |=  ((long)value[offset+intLen-2] << 32);
 1280           long tLong = t & LONG_MASK;
 1281           tLong = tLong * (2 - pLong * tLong);  // 1 more Newton iter step
 1282           tLong = (k == 64 ? tLong : tLong & ((1L << k) - 1));
 1283   
 1284           MutableBigInteger result = new MutableBigInteger(new int[2]);
 1285           result.value[0] = (int)(tLong >>> 32);
 1286           result.value[1] = (int)tLong;
 1287           result.intLen = 2;
 1288           result.normalize();
 1289           return result;
 1290       }
 1291   
 1292       /*
 1293        * Returns the multiplicative inverse of val mod 2^32.  Assumes val is odd.
 1294        */
 1295       static int inverseMod32(int val) {
 1296           // Newton's iteration!
 1297           int t = val;
 1298           t *= 2 - val*t;
 1299           t *= 2 - val*t;
 1300           t *= 2 - val*t;
 1301           t *= 2 - val*t;
 1302           return t;
 1303       }
 1304   
 1305       /*
 1306        * Calculate the multiplicative inverse of 2^k mod mod, where mod is odd.
 1307        */
 1308       static MutableBigInteger modInverseBP2(MutableBigInteger mod, int k) {
 1309           // Copy the mod to protect original
 1310           return fixup(new MutableBigInteger(1), new MutableBigInteger(mod), k);
 1311       }
 1312   
 1313       /**
 1314        * Calculate the multiplicative inverse of this mod mod, where mod is odd.
 1315        * This and mod are not changed by the calculation.
 1316        *
 1317        * This method implements an algorithm due to Richard Schroeppel, that uses
 1318        * the same intermediate representation as Montgomery Reduction
 1319        * ("Montgomery Form").  The algorithm is described in an unpublished
 1320        * manuscript entitled "Fast Modular Reciprocals."
 1321        */
 1322       private MutableBigInteger modInverse(MutableBigInteger mod) {
 1323           MutableBigInteger p = new MutableBigInteger(mod);
 1324           MutableBigInteger f = new MutableBigInteger(this);
 1325           MutableBigInteger g = new MutableBigInteger(p);
 1326           SignedMutableBigInteger c = new SignedMutableBigInteger(1);
 1327           SignedMutableBigInteger d = new SignedMutableBigInteger();
 1328           MutableBigInteger temp = null;
 1329           SignedMutableBigInteger sTemp = null;
 1330   
 1331           int k = 0;
 1332           // Right shift f k times until odd, left shift d k times
 1333           if (f.isEven()) {
 1334               int trailingZeros = f.getLowestSetBit();
 1335               f.rightShift(trailingZeros);
 1336               d.leftShift(trailingZeros);
 1337               k = trailingZeros;
 1338           }
 1339   
 1340           // The Almost Inverse Algorithm
 1341           while(!f.isOne()) {
 1342               // If gcd(f, g) != 1, number is not invertible modulo mod
 1343               if (f.isZero())
 1344                   throw new ArithmeticException("BigInteger not invertible.");
 1345   
 1346               // If f < g exchange f, g and c, d
 1347               if (f.compare(g) < 0) {
 1348                   temp = f; f = g; g = temp;
 1349                   sTemp = d; d = c; c = sTemp;
 1350               }
 1351   
 1352               // If f == g (mod 4)
 1353               if (((f.value[f.offset + f.intLen - 1] ^
 1354                    g.value[g.offset + g.intLen - 1]) & 3) == 0) {
 1355                   f.subtract(g);
 1356                   c.signedSubtract(d);
 1357               } else { // If f != g (mod 4)
 1358                   f.add(g);
 1359                   c.signedAdd(d);
 1360               }
 1361   
 1362               // Right shift f k times until odd, left shift d k times
 1363               int trailingZeros = f.getLowestSetBit();
 1364               f.rightShift(trailingZeros);
 1365               d.leftShift(trailingZeros);
 1366               k += trailingZeros;
 1367           }
 1368   
 1369           while (c.sign < 0)
 1370              c.signedAdd(p);
 1371   
 1372           return fixup(c, p, k);
 1373       }
 1374   
 1375       /*
 1376        * The Fixup Algorithm
 1377        * Calculates X such that X = C * 2^(-k) (mod P)
 1378        * Assumes C<P and P is odd.
 1379        */
 1380       static MutableBigInteger fixup(MutableBigInteger c, MutableBigInteger p,
 1381                                                                         int k) {
 1382           MutableBigInteger temp = new MutableBigInteger();
 1383           // Set r to the multiplicative inverse of p mod 2^32
 1384           int r = -inverseMod32(p.value[p.offset+p.intLen-1]);
 1385   
 1386           for(int i=0, numWords = k >> 5; i<numWords; i++) {
 1387               // V = R * c (mod 2^j)
 1388               int  v = r * c.value[c.offset + c.intLen-1];
 1389               // c = c + (v * p)
 1390               p.mul(v, temp);
 1391               c.add(temp);
 1392               // c = c / 2^j
 1393               c.intLen--;
 1394           }
 1395           int numBits = k & 0x1f;
 1396           if (numBits != 0) {
 1397               // V = R * c (mod 2^j)
 1398               int v = r * c.value[c.offset + c.intLen-1];
 1399               v &= ((1<<numBits) - 1);
 1400               // c = c + (v * p)
 1401               p.mul(v, temp);
 1402               c.add(temp);
 1403               // c = c / 2^j
 1404               c.rightShift(numBits);
 1405           }
 1406   
 1407           // In theory, c may be greater than p at this point (Very rare!)
 1408           while (c.compare(p) >= 0)
 1409               c.subtract(p);
 1410   
 1411           return c;
 1412       }
 1413   
 1414       /**
 1415        * Uses the extended Euclidean algorithm to compute the modInverse of base
 1416        * mod a modulus that is a power of 2. The modulus is 2^k.
 1417        */
 1418       MutableBigInteger euclidModInverse(int k) {
 1419           MutableBigInteger b = new MutableBigInteger(1);
 1420           b.leftShift(k);
 1421           MutableBigInteger mod = new MutableBigInteger(b);
 1422   
 1423           MutableBigInteger a = new MutableBigInteger(this);
 1424           MutableBigInteger q = new MutableBigInteger();
 1425           MutableBigInteger r = b.divide(a, q);
 1426   
 1427           MutableBigInteger swapper = b;
 1428           // swap b & r
 1429           b = r;
 1430           r = swapper;
 1431   
 1432           MutableBigInteger t1 = new MutableBigInteger(q);
 1433           MutableBigInteger t0 = new MutableBigInteger(1);
 1434           MutableBigInteger temp = new MutableBigInteger();
 1435   
 1436           while (!b.isOne()) {
 1437               r = a.divide(b, q);
 1438   
 1439               if (r.intLen == 0)
 1440                   throw new ArithmeticException("BigInteger not invertible.");
 1441   
 1442               swapper = r;
 1443               a = swapper;
 1444   
 1445               if (q.intLen == 1)
 1446                   t1.mul(q.value[q.offset], temp);
 1447               else
 1448                   q.multiply(t1, temp);
 1449               swapper = q;
 1450               q = temp;
 1451               temp = swapper;
 1452               t0.add(q);
 1453   
 1454               if (a.isOne())
 1455                   return t0;
 1456   
 1457               r = b.divide(a, q);
 1458   
 1459               if (r.intLen == 0)
 1460                   throw new ArithmeticException("BigInteger not invertible.");
 1461   
 1462               swapper = b;
 1463               b =  r;
 1464   
 1465               if (q.intLen == 1)
 1466                   t0.mul(q.value[q.offset], temp);
 1467               else
 1468                   q.multiply(t0, temp);
 1469               swapper = q; q = temp; temp = swapper;
 1470   
 1471               t1.add(q);
 1472           }
 1473           mod.subtract(t1);
 1474           return mod;
 1475       }
 1476   
 1477   }

Save This Page
Home » openjdk-7 » java » math » [javadoc | source]