Java > Open Source Codes > java > math > MutableBigInteger


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

A to Z: JavaDoc & Examples

---

Daily Java News & Articles

---

Open Source Projects

---

Open Source Codes

---

Free Computer Books

---

Remove Frame

---

Popular Tags