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