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