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