1 /*
2 * Copyright (c) 2013, Oracle and/or its affiliates. All rights reserved.
3 * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
4 *
5 * This code is free software; you can redistribute it and/or modify it
6 * under the terms of the GNU General Public License version 2 only, as
7 * published by the Free Software Foundation. Oracle designates this
8 * particular file as subject to the "Classpath" exception as provided
9 * by Oracle in the LICENSE file that accompanied this code.
10 *
11 * This code is distributed in the hope that it will be useful, but WITHOUT
12 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
14 * version 2 for more details (a copy is included in the LICENSE file that
15 * accompanied this code).
16 *
17 * You should have received a copy of the GNU General Public License version
18 * 2 along with this work; if not, write to the Free Software Foundation,
19 * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
20 *
21 * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA
22 * or visit www.oracle.com if you need additional information or have any
23 * questions.
24 */
25 package jdk.internal.math;
26
27 import jdk.internal.misc.VM;
28
29 import java.math.BigInteger;
30 import java.util.Arrays;
31 //@ model import org.jmlspecs.models.JMLMath;
32
33 /**
34 * A simple big integer package specifically for floating point base conversion.
35 */
36 public /*@ spec_bigint_math @*/ class FDBigInteger {
37
38 //
39 // This class contains many comments that start with "/*@" mark.
40 // They are behavourial specification in
41 // the Java Modelling Language (JML):
42 // http://www.eecs.ucf.edu/~leavens/JML//index.shtml
43 //
44
45 /*@
46 @ public pure model static \bigint UNSIGNED(int v) {
47 @ return v >= 0 ? v : v + (((\bigint)1) << 32);
48 @ }
49 @
50 @ public pure model static \bigint UNSIGNED(long v) {
51 @ return v >= 0 ? v : v + (((\bigint)1) << 64);
52 @ }
53 @
54 @ public pure model static \bigint AP(int[] data, int len) {
55 @ return (\sum int i; 0 <= 0 && i < len; UNSIGNED(data[i]) << (i*32));
56 @ }
57 @
58 @ public pure model static \bigint pow52(int p5, int p2) {
59 @ ghost \bigint v = 1;
60 @ for (int i = 0; i < p5; i++) v *= 5;
61 @ return v << p2;
62 @ }
63 @
64 @ public pure model static \bigint pow10(int p10) {
65 @ return pow52(p10, p10);
66 @ }
67 @*/
68
69 static final int[] SMALL_5_POW;
70
71 static final long[] LONG_5_POW;
72
73 // Maximum size of cache of powers of 5 as FDBigIntegers.
74 private static final int MAX_FIVE_POW = 340;
75
76 // Cache of big powers of 5 as FDBigIntegers.
77 private static final FDBigInteger POW_5_CACHE[];
78
79 // Archive proxy
80 private static Object[] archivedCaches;
81
82 // Initialize FDBigInteger cache of powers of 5.
83 static {
84 VM.initializeFromArchive(FDBigInteger.class);
85 if (archivedCaches == null) {
86 long[] long5pow = {
87 1L,
88 5L,
89 5L * 5,
90 5L * 5 * 5,
91 5L * 5 * 5 * 5,
92 5L * 5 * 5 * 5 * 5,
93 5L * 5 * 5 * 5 * 5 * 5,
94 5L * 5 * 5 * 5 * 5 * 5 * 5,
95 5L * 5 * 5 * 5 * 5 * 5 * 5 * 5,
96 5L * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5,
97 5L * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5,
98 5L * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5,
99 5L * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5,
100 5L * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5,
101 5L * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5,
102 5L * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5,
103 5L * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5,
104 5L * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5,
105 5L * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5,
106 5L * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5,
107 5L * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5,
108 5L * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5,
109 5L * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5,
110 5L * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5,
111 5L * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5,
112 5L * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5,
113 5L * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5,
114 };
115 int[] small5pow = {
116 1,
117 5,
118 5 * 5,
119 5 * 5 * 5,
120 5 * 5 * 5 * 5,
121 5 * 5 * 5 * 5 * 5,
122 5 * 5 * 5 * 5 * 5 * 5,
123 5 * 5 * 5 * 5 * 5 * 5 * 5,
124 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5,
125 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5,
126 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5,
127 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5,
128 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5,
129 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5
130 };
131 FDBigInteger[] pow5cache = new FDBigInteger[MAX_FIVE_POW];
132 int i = 0;
133 while (i < small5pow.length) {
134 FDBigInteger pow5 = new FDBigInteger(new int[] { small5pow[i] }, 0);
135 pow5.makeImmutable();
136 pow5cache[i] = pow5;
137 i++;
138 }
139 FDBigInteger prev = pow5cache[i - 1];
140 while (i < MAX_FIVE_POW) {
141 pow5cache[i] = prev = prev.mult(5);
142 prev.makeImmutable();
143 i++;
144 }
145 FDBigInteger zero = new FDBigInteger(new int[0], 0);
146 zero.makeImmutable();
147 archivedCaches = new Object[] {small5pow, long5pow, pow5cache, zero};
148 }
149 SMALL_5_POW = (int[])archivedCaches[0];
150 LONG_5_POW = (long[])archivedCaches[1];
151 POW_5_CACHE = (FDBigInteger[])archivedCaches[2];
152 ZERO = (FDBigInteger)archivedCaches[3];
153 }
154
155 // Zero as an FDBigInteger.
156 public static final FDBigInteger ZERO;
157
158 // Constant for casting an int to a long via bitwise AND.
159 private static final long LONG_MASK = 0xffffffffL;
160
161 //@ spec_public non_null;
162 private int data[]; // value: data[0] is least significant
163 //@ spec_public;
164 private int offset; // number of least significant zero padding ints
165 //@ spec_public;
166 private int nWords; // data[nWords-1]!=0, all values above are zero
167 // if nWords==0 -> this FDBigInteger is zero
168 //@ spec_public;
169 private boolean isImmutable = false;
170
171 /*@
172 @ public invariant 0 <= nWords && nWords <= data.length && offset >= 0;
173 @ public invariant nWords == 0 ==> offset == 0;
174 @ public invariant nWords > 0 ==> data[nWords - 1] != 0;
175 @ public invariant (\forall int i; nWords <= i && i < data.length; data[i] == 0);
176 @ public pure model \bigint value() {
177 @ return AP(data, nWords) << (offset*32);
178 @ }
179 @*/
180
181 /**
182 * Constructs an <code>FDBigInteger</code> from data and padding. The
183 * <code>data</code> parameter has the least significant <code>int</code> at
184 * the zeroth index. The <code>offset</code> parameter gives the number of
185 * zero <code>int</code>s to be inferred below the least significant element
186 * of <code>data</code>.
187 *
188 * @param data An array containing all non-zero <code>int</code>s of the value.
189 * @param offset An offset indicating the number of zero <code>int</code>s to pad
190 * below the least significant element of <code>data</code>.
191 */
192 /*@
193 @ requires data != null && offset >= 0;
194 @ ensures this.value() == \old(AP(data, data.length) << (offset*32));
195 @ ensures this.data == \old(data);
196 @*/
197 private FDBigInteger(int[] data, int offset) {
198 this.data = data;
199 this.offset = offset;
200 this.nWords = data.length;
201 trimLeadingZeros();
202 }
203
204 /**
205 * Constructs an <code>FDBigInteger</code> from a starting value and some
206 * decimal digits.
207 *
208 * @param lValue The starting value.
209 * @param digits The decimal digits.
210 * @param kDigits The initial index into <code>digits</code>.
211 * @param nDigits The final index into <code>digits</code>.
212 */
213 /*@
214 @ requires digits != null;
215 @ requires 0 <= kDigits && kDigits <= nDigits && nDigits <= digits.length;
216 @ requires (\forall int i; 0 <= i && i < nDigits; '0' <= digits[i] && digits[i] <= '9');
217 @ ensures this.value() == \old(lValue * pow10(nDigits - kDigits) + (\sum int i; kDigits <= i && i < nDigits; (digits[i] - '0') * pow10(nDigits - i - 1)));
218 @*/
219 public FDBigInteger(long lValue, char[] digits, int kDigits, int nDigits) {
220 int n = Math.max((nDigits + 8) / 9, 2); // estimate size needed.
221 data = new int[n]; // allocate enough space
222 data[0] = (int) lValue; // starting value
223 data[1] = (int) (lValue >>> 32);
224 offset = 0;
225 nWords = 2;
226 int i = kDigits;
227 int limit = nDigits - 5; // slurp digits 5 at a time.
228 int v;
229 while (i < limit) {
230 int ilim = i + 5;
231 v = (int) digits[i++] - (int) '0';
232 while (i < ilim) {
233 v = 10 * v + (int) digits[i++] - (int) '0';
234 }
235 multAddMe(100000, v); // ... where 100000 is 10^5.
236 }
237 int factor = 1;
238 v = 0;
239 while (i < nDigits) {
240 v = 10 * v + (int) digits[i++] - (int) '0';
241 factor *= 10;
242 }
243 if (factor != 1) {
244 multAddMe(factor, v);
245 }
246 trimLeadingZeros();
247 }
248
249 /**
250 * Returns an <code>FDBigInteger</code> with the numerical value
251 * <code>5<sup>p5</sup> * 2<sup>p2</sup></code>.
252 *
253 * @param p5 The exponent of the power-of-five factor.
254 * @param p2 The exponent of the power-of-two factor.
255 * @return <code>5<sup>p5</sup> * 2<sup>p2</sup></code>
256 */
257 /*@
258 @ requires p5 >= 0 && p2 >= 0;
259 @ assignable \nothing;
260 @ ensures \result.value() == \old(pow52(p5, p2));
261 @*/
262 public static FDBigInteger valueOfPow52(int p5, int p2) {
263 if (p5 != 0) {
264 if (p2 == 0) {
265 return big5pow(p5);
266 } else if (p5 < SMALL_5_POW.length) {
267 int pow5 = SMALL_5_POW[p5];
268 int wordcount = p2 >> 5;
269 int bitcount = p2 & 0x1f;
270 if (bitcount == 0) {
271 return new FDBigInteger(new int[]{pow5}, wordcount);
272 } else {
273 return new FDBigInteger(new int[]{
274 pow5 << bitcount,
275 pow5 >>> (32 - bitcount)
276 }, wordcount);
277 }
278 } else {
279 return big5pow(p5).leftShift(p2);
280 }
281 } else {
282 return valueOfPow2(p2);
283 }
284 }
285
286 /**
287 * Returns an <code>FDBigInteger</code> with the numerical value
288 * <code>value * 5<sup>p5</sup> * 2<sup>p2</sup></code>.
289 *
290 * @param value The constant factor.
291 * @param p5 The exponent of the power-of-five factor.
292 * @param p2 The exponent of the power-of-two factor.
293 * @return <code>value * 5<sup>p5</sup> * 2<sup>p2</sup></code>
294 */
295 /*@
296 @ requires p5 >= 0 && p2 >= 0;
297 @ assignable \nothing;
298 @ ensures \result.value() == \old(UNSIGNED(value) * pow52(p5, p2));
299 @*/
300 public static FDBigInteger valueOfMulPow52(long value, int p5, int p2) {
301 assert p5 >= 0 : p5;
302 assert p2 >= 0 : p2;
303 int v0 = (int) value;
304 int v1 = (int) (value >>> 32);
305 int wordcount = p2 >> 5;
306 int bitcount = p2 & 0x1f;
307 if (p5 != 0) {
308 if (p5 < SMALL_5_POW.length) {
309 long pow5 = SMALL_5_POW[p5] & LONG_MASK;
310 long carry = (v0 & LONG_MASK) * pow5;
311 v0 = (int) carry;
312 carry >>>= 32;
313 carry = (v1 & LONG_MASK) * pow5 + carry;
314 v1 = (int) carry;
315 int v2 = (int) (carry >>> 32);
316 if (bitcount == 0) {
317 return new FDBigInteger(new int[]{v0, v1, v2}, wordcount);
318 } else {
319 return new FDBigInteger(new int[]{
320 v0 << bitcount,
321 (v1 << bitcount) | (v0 >>> (32 - bitcount)),
322 (v2 << bitcount) | (v1 >>> (32 - bitcount)),
323 v2 >>> (32 - bitcount)
324 }, wordcount);
325 }
326 } else {
327 FDBigInteger pow5 = big5pow(p5);
328 int[] r;
329 if (v1 == 0) {
330 r = new int[pow5.nWords + 1 + ((p2 != 0) ? 1 : 0)];
331 mult(pow5.data, pow5.nWords, v0, r);
332 } else {
333 r = new int[pow5.nWords + 2 + ((p2 != 0) ? 1 : 0)];
334 mult(pow5.data, pow5.nWords, v0, v1, r);
335 }
336 return (new FDBigInteger(r, pow5.offset)).leftShift(p2);
337 }
338 } else if (p2 != 0) {
339 if (bitcount == 0) {
340 return new FDBigInteger(new int[]{v0, v1}, wordcount);
341 } else {
342 return new FDBigInteger(new int[]{
343 v0 << bitcount,
344 (v1 << bitcount) | (v0 >>> (32 - bitcount)),
345 v1 >>> (32 - bitcount)
346 }, wordcount);
347 }
348 }
349 return new FDBigInteger(new int[]{v0, v1}, 0);
350 }
351
352 /**
353 * Returns an <code>FDBigInteger</code> with the numerical value
354 * <code>2<sup>p2</sup></code>.
355 *
356 * @param p2 The exponent of 2.
357 * @return <code>2<sup>p2</sup></code>
358 */
359 /*@
360 @ requires p2 >= 0;
361 @ assignable \nothing;
362 @ ensures \result.value() == pow52(0, p2);
363 @*/
364 private static FDBigInteger valueOfPow2(int p2) {
365 int wordcount = p2 >> 5;
366 int bitcount = p2 & 0x1f;
367 return new FDBigInteger(new int[]{1 << bitcount}, wordcount);
368 }
369
370 /**
371 * Removes all leading zeros from this <code>FDBigInteger</code> adjusting
372 * the offset and number of non-zero leading words accordingly.
373 */
374 /*@
375 @ requires data != null;
376 @ requires 0 <= nWords && nWords <= data.length && offset >= 0;
377 @ requires nWords == 0 ==> offset == 0;
378 @ ensures nWords == 0 ==> offset == 0;
379 @ ensures nWords > 0 ==> data[nWords - 1] != 0;
380 @*/
381 private /*@ helper @*/ void trimLeadingZeros() {
382 int i = nWords;
383 if (i > 0 && (data[--i] == 0)) {
384 //for (; i > 0 && data[i - 1] == 0; i--) ;
385 while(i > 0 && data[i - 1] == 0) {
386 i--;
387 }
388 this.nWords = i;
389 if (i == 0) { // all words are zero
390 this.offset = 0;
391 }
392 }
393 }
394
395 /**
396 * Retrieves the normalization bias of the <code>FDBigIntger</code>. The
397 * normalization bias is a left shift such that after it the highest word
398 * of the value will have the 4 highest bits equal to zero:
399 * {@code (highestWord & 0xf0000000) == 0}, but the next bit should be 1
400 * {@code (highestWord & 0x08000000) != 0}.
401 *
402 * @return The normalization bias.
403 */
404 /*@
405 @ requires this.value() > 0;
406 @*/
407 public /*@ pure @*/ int getNormalizationBias() {
408 if (nWords == 0) {
409 throw new IllegalArgumentException("Zero value cannot be normalized");
410 }
411 int zeros = Integer.numberOfLeadingZeros(data[nWords - 1]);
412 return (zeros < 4) ? 28 + zeros : zeros - 4;
413 }
414
415 // TODO: Why is anticount param needed if it is always 32 - bitcount?
416 /**
417 * Left shifts the contents of one int array into another.
418 *
419 * @param src The source array.
420 * @param idx The initial index of the source array.
421 * @param result The destination array.
422 * @param bitcount The left shift.
423 * @param anticount The left anti-shift, e.g., <code>32-bitcount</code>.
424 * @param prev The prior source value.
425 */
426 /*@
427 @ requires 0 < bitcount && bitcount < 32 && anticount == 32 - bitcount;
428 @ requires src.length >= idx && result.length > idx;
429 @ assignable result[*];
430 @ ensures AP(result, \old(idx + 1)) == \old((AP(src, idx) + UNSIGNED(prev) << (idx*32)) << bitcount);
431 @*/
432 private static void leftShift(int[] src, int idx, int result[], int bitcount, int anticount, int prev){
433 for (; idx > 0; idx--) {
434 int v = (prev << bitcount);
435 prev = src[idx - 1];
436 v |= (prev >>> anticount);
437 result[idx] = v;
438 }
439 int v = prev << bitcount;
440 result[0] = v;
441 }
442
443 /**
444 * Shifts this <code>FDBigInteger</code> to the left. The shift is performed
445 * in-place unless the <code>FDBigInteger</code> is immutable in which case
446 * a new instance of <code>FDBigInteger</code> is returned.
447 *
448 * @param shift The number of bits to shift left.
449 * @return The shifted <code>FDBigInteger</code>.
450 */
451 /*@
452 @ requires this.value() == 0 || shift == 0;
453 @ assignable \nothing;
454 @ ensures \result == this;
455 @
456 @ also
457 @
458 @ requires this.value() > 0 && shift > 0 && this.isImmutable;
459 @ assignable \nothing;
460 @ ensures \result.value() == \old(this.value() << shift);
461 @
462 @ also
463 @
464 @ requires this.value() > 0 && shift > 0 && this.isImmutable;
465 @ assignable \nothing;
466 @ ensures \result == this;
467 @ ensures \result.value() == \old(this.value() << shift);
468 @*/
469 public FDBigInteger leftShift(int shift) {
470 if (shift == 0 || nWords == 0) {
471 return this;
472 }
473 int wordcount = shift >> 5;
474 int bitcount = shift & 0x1f;
475 if (this.isImmutable) {
476 if (bitcount == 0) {
477 return new FDBigInteger(Arrays.copyOf(data, nWords), offset + wordcount);
478 } else {
479 int anticount = 32 - bitcount;
480 int idx = nWords - 1;
481 int prev = data[idx];
482 int hi = prev >>> anticount;
483 int[] result;
484 if (hi != 0) {
485 result = new int[nWords + 1];
486 result[nWords] = hi;
487 } else {
488 result = new int[nWords];
489 }
490 leftShift(data,idx,result,bitcount,anticount,prev);
491 return new FDBigInteger(result, offset + wordcount);
492 }
493 } else {
494 if (bitcount != 0) {
495 int anticount = 32 - bitcount;
496 if ((data[0] << bitcount) == 0) {
497 int idx = 0;
498 int prev = data[idx];
499 for (; idx < nWords - 1; idx++) {
500 int v = (prev >>> anticount);
501 prev = data[idx + 1];
502 v |= (prev << bitcount);
503 data[idx] = v;
504 }
505 int v = prev >>> anticount;
506 data[idx] = v;
507 if(v==0) {
508 nWords--;
509 }
510 offset++;
511 } else {
512 int idx = nWords - 1;
513 int prev = data[idx];
514 int hi = prev >>> anticount;
515 int[] result = data;
516 int[] src = data;
517 if (hi != 0) {
518 if(nWords == data.length) {
519 data = result = new int[nWords + 1];
520 }
521 result[nWords++] = hi;
522 }
523 leftShift(src,idx,result,bitcount,anticount,prev);
524 }
525 }
526 offset += wordcount;
527 return this;
528 }
529 }
530
531 /**
532 * Returns the number of <code>int</code>s this <code>FDBigInteger</code> represents.
533 *
534 * @return Number of <code>int</code>s required to represent this <code>FDBigInteger</code>.
535 */
536 /*@
537 @ requires this.value() == 0;
538 @ ensures \result == 0;
539 @
540 @ also
541 @
542 @ requires this.value() > 0;
543 @ ensures ((\bigint)1) << (\result - 1) <= this.value() && this.value() <= ((\bigint)1) << \result;
544 @*/
545 private /*@ pure @*/ int size() {
546 return nWords + offset;
547 }
548
549
550 /**
551 * Computes
552 * <pre>
553 * q = (int)( this / S )
554 * this = 10 * ( this mod S )
555 * Return q.
556 * </pre>
557 * This is the iteration step of digit development for output.
558 * We assume that S has been normalized, as above, and that
559 * "this" has been left-shifted accordingly.
560 * Also assumed, of course, is that the result, q, can be expressed
561 * as an integer, {@code 0 <= q < 10}.
562 *
563 * @param S The divisor of this <code>FDBigInteger</code>.
564 * @return <code>q = (int)(this / S)</code>.
565 */
566 /*@
567 @ requires !this.isImmutable;
568 @ requires this.size() <= S.size();
569 @ requires this.data.length + this.offset >= S.size();
570 @ requires S.value() >= ((\bigint)1) << (S.size()*32 - 4);
571 @ assignable this.nWords, this.offset, this.data, this.data[*];
572 @ ensures \result == \old(this.value() / S.value());
573 @ ensures this.value() == \old(10 * (this.value() % S.value()));
574 @*/
575 public int quoRemIteration(FDBigInteger S) throws IllegalArgumentException {
576 assert !this.isImmutable : "cannot modify immutable value";
577 // ensure that this and S have the same number of
578 // digits. If S is properly normalized and q < 10 then
579 // this must be so.
580 int thSize = this.size();
581 int sSize = S.size();
582 if (thSize < sSize) {
583 // this value is significantly less than S, result of division is zero.
584 // just mult this by 10.
585 int p = multAndCarryBy10(this.data, this.nWords, this.data);
586 if(p!=0) {
587 this.data[nWords++] = p;
588 } else {
589 trimLeadingZeros();
590 }
591 return 0;
592 } else if (thSize > sSize) {
593 throw new IllegalArgumentException("disparate values");
594 }
595 // estimate q the obvious way. We will usually be
596 // right. If not, then we're only off by a little and
597 // will re-add.
598 long q = (this.data[this.nWords - 1] & LONG_MASK) / (S.data[S.nWords - 1] & LONG_MASK);
599 long diff = multDiffMe(q, S);
600 if (diff != 0L) {
601 //@ assert q != 0;
602 //@ assert this.offset == \old(Math.min(this.offset, S.offset));
603 //@ assert this.offset <= S.offset;
604
605 // q is too big.
606 // add S back in until this turns +. This should
607 // not be very many times!
608 long sum = 0L;
609 int tStart = S.offset - this.offset;
610 //@ assert tStart >= 0;
611 int[] sd = S.data;
612 int[] td = this.data;
613 while (sum == 0L) {
614 for (int sIndex = 0, tIndex = tStart; tIndex < this.nWords; sIndex++, tIndex++) {
615 sum += (td[tIndex] & LONG_MASK) + (sd[sIndex] & LONG_MASK);
616 td[tIndex] = (int) sum;
617 sum >>>= 32; // Signed or unsigned, answer is 0 or 1
618 }
619 //
620 // Originally the following line read
621 // "if ( sum !=0 && sum != -1 )"
622 // but that would be wrong, because of the
623 // treatment of the two values as entirely unsigned,
624 // it would be impossible for a carry-out to be interpreted
625 // as -1 -- it would have to be a single-bit carry-out, or +1.
626 //
627 assert sum == 0 || sum == 1 : sum; // carry out of division correction
628 q -= 1;
629 }
630 }
631 // finally, we can multiply this by 10.
632 // it cannot overflow, right, as the high-order word has
633 // at least 4 high-order zeros!
634 int p = multAndCarryBy10(this.data, this.nWords, this.data);
635 assert p == 0 : p; // Carry out of *10
636 trimLeadingZeros();
637 return (int) q;
638 }
639
640 /**
641 * Multiplies this <code>FDBigInteger</code> by 10. The operation will be
642 * performed in place unless the <code>FDBigInteger</code> is immutable in
643 * which case a new <code>FDBigInteger</code> will be returned.
644 *
645 * @return The <code>FDBigInteger</code> multiplied by 10.
646 */
647 /*@
648 @ requires this.value() == 0;
649 @ assignable \nothing;
650 @ ensures \result == this;
651 @
652 @ also
653 @
654 @ requires this.value() > 0 && this.isImmutable;
655 @ assignable \nothing;
656 @ ensures \result.value() == \old(this.value() * 10);
657 @
658 @ also
659 @
660 @ requires this.value() > 0 && !this.isImmutable;
661 @ assignable this.nWords, this.data, this.data[*];
662 @ ensures \result == this;
663 @ ensures \result.value() == \old(this.value() * 10);
664 @*/
665 public FDBigInteger multBy10() {
666 if (nWords == 0) {
667 return this;
668 }
669 if (isImmutable) {
670 int[] res = new int[nWords + 1];
671 res[nWords] = multAndCarryBy10(data, nWords, res);
672 return new FDBigInteger(res, offset);
673 } else {
674 int p = multAndCarryBy10(this.data, this.nWords, this.data);
675 if (p != 0) {
676 if (nWords == data.length) {
677 if (data[0] == 0) {
678 System.arraycopy(data, 1, data, 0, --nWords);
679 offset++;
680 } else {
681 data = Arrays.copyOf(data, data.length + 1);
682 }
683 }
684 data[nWords++] = p;
685 } else {
686 trimLeadingZeros();
687 }
688 return this;
689 }
690 }
691
692 /**
693 * Multiplies this <code>FDBigInteger</code> by
694 * <code>5<sup>p5</sup> * 2<sup>p2</sup></code>. The operation will be
695 * performed in place if possible, otherwise a new <code>FDBigInteger</code>
696 * will be returned.
697 *
698 * @param p5 The exponent of the power-of-five factor.
699 * @param p2 The exponent of the power-of-two factor.
700 * @return The multiplication result.
701 */
702 /*@
703 @ requires this.value() == 0 || p5 == 0 && p2 == 0;
704 @ assignable \nothing;
705 @ ensures \result == this;
706 @
707 @ also
708 @
709 @ requires this.value() > 0 && (p5 > 0 && p2 >= 0 || p5 == 0 && p2 > 0 && this.isImmutable);
710 @ assignable \nothing;
711 @ ensures \result.value() == \old(this.value() * pow52(p5, p2));
712 @
713 @ also
714 @
715 @ requires this.value() > 0 && p5 == 0 && p2 > 0 && !this.isImmutable;
716 @ assignable this.nWords, this.data, this.data[*];
717 @ ensures \result == this;
718 @ ensures \result.value() == \old(this.value() * pow52(p5, p2));
719 @*/
720 public FDBigInteger multByPow52(int p5, int p2) {
721 if (this.nWords == 0) {
722 return this;
723 }
724 FDBigInteger res = this;
725 if (p5 != 0) {
726 int[] r;
727 int extraSize = (p2 != 0) ? 1 : 0;
728 if (p5 < SMALL_5_POW.length) {
729 r = new int[this.nWords + 1 + extraSize];
730 mult(this.data, this.nWords, SMALL_5_POW[p5], r);
731 res = new FDBigInteger(r, this.offset);
732 } else {
733 FDBigInteger pow5 = big5pow(p5);
734 r = new int[this.nWords + pow5.size() + extraSize];
735 mult(this.data, this.nWords, pow5.data, pow5.nWords, r);
736 res = new FDBigInteger(r, this.offset + pow5.offset);
737 }
738 }
739 return res.leftShift(p2);
740 }
741
742 /**
743 * Multiplies two big integers represented as int arrays.
744 *
745 * @param s1 The first array factor.
746 * @param s1Len The number of elements of <code>s1</code> to use.
747 * @param s2 The second array factor.
748 * @param s2Len The number of elements of <code>s2</code> to use.
749 * @param dst The product array.
750 */
751 /*@
752 @ requires s1 != dst && s2 != dst;
753 @ requires s1.length >= s1Len && s2.length >= s2Len && dst.length >= s1Len + s2Len;
754 @ assignable dst[0 .. s1Len + s2Len - 1];
755 @ ensures AP(dst, s1Len + s2Len) == \old(AP(s1, s1Len) * AP(s2, s2Len));
756 @*/
757 private static void mult(int[] s1, int s1Len, int[] s2, int s2Len, int[] dst) {
758 for (int i = 0; i < s1Len; i++) {
759 long v = s1[i] & LONG_MASK;
760 long p = 0L;
761 for (int j = 0; j < s2Len; j++) {
762 p += (dst[i + j] & LONG_MASK) + v * (s2[j] & LONG_MASK);
763 dst[i + j] = (int) p;
764 p >>>= 32;
765 }
766 dst[i + s2Len] = (int) p;
767 }
768 }
769
770 /**
771 * Subtracts the supplied <code>FDBigInteger</code> subtrahend from this
772 * <code>FDBigInteger</code>. Assert that the result is positive.
773 * If the subtrahend is immutable, store the result in this(minuend).
774 * If this(minuend) is immutable a new <code>FDBigInteger</code> is created.
775 *
776 * @param subtrahend The <code>FDBigInteger</code> to be subtracted.
777 * @return This <code>FDBigInteger</code> less the subtrahend.
778 */
779 /*@
780 @ requires this.isImmutable;
781 @ requires this.value() >= subtrahend.value();
782 @ assignable \nothing;
783 @ ensures \result.value() == \old(this.value() - subtrahend.value());
784 @
785 @ also
786 @
787 @ requires !subtrahend.isImmutable;
788 @ requires this.value() >= subtrahend.value();
789 @ assignable this.nWords, this.offset, this.data, this.data[*];
790 @ ensures \result == this;
791 @ ensures \result.value() == \old(this.value() - subtrahend.value());
792 @*/
793 public FDBigInteger leftInplaceSub(FDBigInteger subtrahend) {
794 assert this.size() >= subtrahend.size() : "result should be positive";
795 FDBigInteger minuend;
796 if (this.isImmutable) {
797 minuend = new FDBigInteger(this.data.clone(), this.offset);
798 } else {
799 minuend = this;
800 }
801 int offsetDiff = subtrahend.offset - minuend.offset;
802 int[] sData = subtrahend.data;
803 int[] mData = minuend.data;
804 int subLen = subtrahend.nWords;
805 int minLen = minuend.nWords;
806 if (offsetDiff < 0) {
807 // need to expand minuend
808 int rLen = minLen - offsetDiff;
809 if (rLen < mData.length) {
810 System.arraycopy(mData, 0, mData, -offsetDiff, minLen);
811 Arrays.fill(mData, 0, -offsetDiff, 0);
812 } else {
813 int[] r = new int[rLen];
814 System.arraycopy(mData, 0, r, -offsetDiff, minLen);
815 minuend.data = mData = r;
816 }
817 minuend.offset = subtrahend.offset;
818 minuend.nWords = minLen = rLen;
819 offsetDiff = 0;
820 }
821 long borrow = 0L;
822 int mIndex = offsetDiff;
823 for (int sIndex = 0; sIndex < subLen && mIndex < minLen; sIndex++, mIndex++) {
824 long diff = (mData[mIndex] & LONG_MASK) - (sData[sIndex] & LONG_MASK) + borrow;
825 mData[mIndex] = (int) diff;
826 borrow = diff >> 32; // signed shift
827 }
828 for (; borrow != 0 && mIndex < minLen; mIndex++) {
829 long diff = (mData[mIndex] & LONG_MASK) + borrow;
830 mData[mIndex] = (int) diff;
831 borrow = diff >> 32; // signed shift
832 }
833 assert borrow == 0L : borrow; // borrow out of subtract,
834 // result should be positive
835 minuend.trimLeadingZeros();
836 return minuend;
837 }
838
839 /**
840 * Subtracts the supplied <code>FDBigInteger</code> subtrahend from this
841 * <code>FDBigInteger</code>. Assert that the result is positive.
842 * If the this(minuend) is immutable, store the result in subtrahend.
843 * If subtrahend is immutable a new <code>FDBigInteger</code> is created.
844 *
845 * @param subtrahend The <code>FDBigInteger</code> to be subtracted.
846 * @return This <code>FDBigInteger</code> less the subtrahend.
847 */
848 /*@
849 @ requires subtrahend.isImmutable;
850 @ requires this.value() >= subtrahend.value();
851 @ assignable \nothing;
852 @ ensures \result.value() == \old(this.value() - subtrahend.value());
853 @
854 @ also
855 @
856 @ requires !subtrahend.isImmutable;
857 @ requires this.value() >= subtrahend.value();
858 @ assignable subtrahend.nWords, subtrahend.offset, subtrahend.data, subtrahend.data[*];
859 @ ensures \result == subtrahend;
860 @ ensures \result.value() == \old(this.value() - subtrahend.value());
861 @*/
862 public FDBigInteger rightInplaceSub(FDBigInteger subtrahend) {
863 assert this.size() >= subtrahend.size() : "result should be positive";
864 FDBigInteger minuend = this;
865 if (subtrahend.isImmutable) {
866 subtrahend = new FDBigInteger(subtrahend.data.clone(), subtrahend.offset);
867 }
868 int offsetDiff = minuend.offset - subtrahend.offset;
869 int[] sData = subtrahend.data;
870 int[] mData = minuend.data;
871 int subLen = subtrahend.nWords;
872 int minLen = minuend.nWords;
873 if (offsetDiff < 0) {
874 int rLen = minLen;
875 if (rLen < sData.length) {
876 System.arraycopy(sData, 0, sData, -offsetDiff, subLen);
877 Arrays.fill(sData, 0, -offsetDiff, 0);
878 } else {
879 int[] r = new int[rLen];
880 System.arraycopy(sData, 0, r, -offsetDiff, subLen);
881 subtrahend.data = sData = r;
882 }
883 subtrahend.offset = minuend.offset;
884 subLen -= offsetDiff;
885 offsetDiff = 0;
886 } else {
887 int rLen = minLen + offsetDiff;
888 if (rLen >= sData.length) {
889 subtrahend.data = sData = Arrays.copyOf(sData, rLen);
890 }
891 }
892 //@ assert minuend == this && minuend.value() == \old(this.value());
893 //@ assert mData == minuend.data && minLen == minuend.nWords;
894 //@ assert subtrahend.offset + subtrahend.data.length >= minuend.size();
895 //@ assert sData == subtrahend.data;
896 //@ assert AP(subtrahend.data, subtrahend.data.length) << subtrahend.offset == \old(subtrahend.value());
897 //@ assert subtrahend.offset == Math.min(\old(this.offset), minuend.offset);
898 //@ assert offsetDiff == minuend.offset - subtrahend.offset;
899 //@ assert 0 <= offsetDiff && offsetDiff + minLen <= sData.length;
900 int sIndex = 0;
901 long borrow = 0L;
902 for (; sIndex < offsetDiff; sIndex++) {
903 long diff = 0L - (sData[sIndex] & LONG_MASK) + borrow;
904 sData[sIndex] = (int) diff;
905 borrow = diff >> 32; // signed shift
906 }
907 //@ assert sIndex == offsetDiff;
908 for (int mIndex = 0; mIndex < minLen; sIndex++, mIndex++) {
909 //@ assert sIndex == offsetDiff + mIndex;
910 long diff = (mData[mIndex] & LONG_MASK) - (sData[sIndex] & LONG_MASK) + borrow;
911 sData[sIndex] = (int) diff;
912 borrow = diff >> 32; // signed shift
913 }
914 assert borrow == 0L : borrow; // borrow out of subtract,
915 // result should be positive
916 subtrahend.nWords = sIndex;
917 subtrahend.trimLeadingZeros();
918 return subtrahend;
919
920 }
921
922 /**
923 * Determines whether all elements of an array are zero for all indices less
924 * than a given index.
925 *
926 * @param a The array to be examined.
927 * @param from The index strictly below which elements are to be examined.
928 * @return Zero if all elements in range are zero, 1 otherwise.
929 */
930 /*@
931 @ requires 0 <= from && from <= a.length;
932 @ ensures \result == (AP(a, from) == 0 ? 0 : 1);
933 @*/
934 private /*@ pure @*/ static int checkZeroTail(int[] a, int from) {
935 while (from > 0) {
936 if (a[--from] != 0) {
937 return 1;
938 }
939 }
940 return 0;
941 }
942
943 /**
944 * Compares the parameter with this <code>FDBigInteger</code>. Returns an
945 * integer accordingly as:
946 * <pre>{@code
947 * > 0: this > other
948 * 0: this == other
949 * < 0: this < other
950 * }</pre>
951 *
952 * @param other The <code>FDBigInteger</code> to compare.
953 * @return A negative value, zero, or a positive value according to the
954 * result of the comparison.
955 */
956 /*@
957 @ ensures \result == (this.value() < other.value() ? -1 : this.value() > other.value() ? +1 : 0);
958 @*/
959 public /*@ pure @*/ int cmp(FDBigInteger other) {
960 int aSize = nWords + offset;
961 int bSize = other.nWords + other.offset;
962 if (aSize > bSize) {
963 return 1;
964 } else if (aSize < bSize) {
965 return -1;
966 }
967 int aLen = nWords;
968 int bLen = other.nWords;
969 while (aLen > 0 && bLen > 0) {
970 int a = data[--aLen];
971 int b = other.data[--bLen];
972 if (a != b) {
973 return ((a & LONG_MASK) < (b & LONG_MASK)) ? -1 : 1;
974 }
975 }
976 if (aLen > 0) {
977 return checkZeroTail(data, aLen);
978 }
979 if (bLen > 0) {
980 return -checkZeroTail(other.data, bLen);
981 }
982 return 0;
983 }
984
985 /**
986 * Compares this <code>FDBigInteger</code> with
987 * <code>5<sup>p5</sup> * 2<sup>p2</sup></code>.
988 * Returns an integer accordingly as:
989 * <pre>{@code
990 * > 0: this > other
991 * 0: this == other
992 * < 0: this < other
993 * }</pre>
994 * @param p5 The exponent of the power-of-five factor.
995 * @param p2 The exponent of the power-of-two factor.
996 * @return A negative value, zero, or a positive value according to the
997 * result of the comparison.
998 */
999 /*@
1000 @ requires p5 >= 0 && p2 >= 0;
1001 @ ensures \result == (this.value() < pow52(p5, p2) ? -1 : this.value() > pow52(p5, p2) ? +1 : 0);
1002 @*/
1003 public /*@ pure @*/ int cmpPow52(int p5, int p2) {
1004 if (p5 == 0) {
1005 int wordcount = p2 >> 5;
1006 int bitcount = p2 & 0x1f;
1007 int size = this.nWords + this.offset;
1008 if (size > wordcount + 1) {
1009 return 1;
1010 } else if (size < wordcount + 1) {
1011 return -1;
1012 }
1013 int a = this.data[this.nWords -1];
1014 int b = 1 << bitcount;
1015 if (a != b) {
1016 return ( (a & LONG_MASK) < (b & LONG_MASK)) ? -1 : 1;
1017 }
1018 return checkZeroTail(this.data, this.nWords - 1);
1019 }
1020 return this.cmp(big5pow(p5).leftShift(p2));
1021 }
1022
1023 /**
1024 * Compares this <code>FDBigInteger</code> with <code>x + y</code>. Returns a
1025 * value according to the comparison as:
1026 * <pre>{@code
1027 * -1: this < x + y
1028 * 0: this == x + y
1029 * 1: this > x + y
1030 * }</pre>
1031 * @param x The first addend of the sum to compare.
1032 * @param y The second addend of the sum to compare.
1033 * @return -1, 0, or 1 according to the result of the comparison.
1034 */
1035 /*@
1036 @ ensures \result == (this.value() < x.value() + y.value() ? -1 : this.value() > x.value() + y.value() ? +1 : 0);
1037 @*/
1038 public /*@ pure @*/ int addAndCmp(FDBigInteger x, FDBigInteger y) {
1039 FDBigInteger big;
1040 FDBigInteger small;
1041 int xSize = x.size();
1042 int ySize = y.size();
1043 int bSize;
1044 int sSize;
1045 if (xSize >= ySize) {
1046 big = x;
1047 small = y;
1048 bSize = xSize;
1049 sSize = ySize;
1050 } else {
1051 big = y;
1052 small = x;
1053 bSize = ySize;
1054 sSize = xSize;
1055 }
1056 int thSize = this.size();
1057 if (bSize == 0) {
1058 return thSize == 0 ? 0 : 1;
1059 }
1060 if (sSize == 0) {
1061 return this.cmp(big);
1062 }
1063 if (bSize > thSize) {
1064 return -1;
1065 }
1066 if (bSize + 1 < thSize) {
1067 return 1;
1068 }
1069 long top = (big.data[big.nWords - 1] & LONG_MASK);
1070 if (sSize == bSize) {
1071 top += (small.data[small.nWords - 1] & LONG_MASK);
1072 }
1073 if ((top >>> 32) == 0) {
1074 if (((top + 1) >>> 32) == 0) {
1075 // good case - no carry extension
1076 if (bSize < thSize) {
1077 return 1;
1078 }
1079 // here sum.nWords == this.nWords
1080 long v = (this.data[this.nWords - 1] & LONG_MASK);
1081 if (v < top) {
1082 return -1;
1083 }
1084 if (v > top + 1) {
1085 return 1;
1086 }
1087 }
1088 } else { // (top>>>32)!=0 guaranteed carry extension
1089 if (bSize + 1 > thSize) {
1090 return -1;
1091 }
1092 // here sum.nWords == this.nWords
1093 top >>>= 32;
1094 long v = (this.data[this.nWords - 1] & LONG_MASK);
1095 if (v < top) {
1096 return -1;
1097 }
1098 if (v > top + 1) {
1099 return 1;
1100 }
1101 }
1102 return this.cmp(big.add(small));
1103 }
1104
1105 /**
1106 * Makes this <code>FDBigInteger</code> immutable.
1107 */
1108 /*@
1109 @ assignable this.isImmutable;
1110 @ ensures this.isImmutable;
1111 @*/
1112 public void makeImmutable() {
1113 this.isImmutable = true;
1114 }
1115
1116 /**
1117 * Multiplies this <code>FDBigInteger</code> by an integer.
1118 *
1119 * @param i The factor by which to multiply this <code>FDBigInteger</code>.
1120 * @return This <code>FDBigInteger</code> multiplied by an integer.
1121 */
1122 /*@
1123 @ requires this.value() == 0;
1124 @ assignable \nothing;
1125 @ ensures \result == this;
1126 @
1127 @ also
1128 @
1129 @ requires this.value() != 0;
1130 @ assignable \nothing;
1131 @ ensures \result.value() == \old(this.value() * UNSIGNED(i));
1132 @*/
1133 private FDBigInteger mult(int i) {
1134 if (this.nWords == 0) {
1135 return this;
1136 }
1137 int[] r = new int[nWords + 1];
1138 mult(data, nWords, i, r);
1139 return new FDBigInteger(r, offset);
1140 }
1141
1142 /**
1143 * Multiplies this <code>FDBigInteger</code> by another <code>FDBigInteger</code>.
1144 *
1145 * @param other The <code>FDBigInteger</code> factor by which to multiply.
1146 * @return The product of this and the parameter <code>FDBigInteger</code>s.
1147 */
1148 /*@
1149 @ requires this.value() == 0;
1150 @ assignable \nothing;
1151 @ ensures \result == this;
1152 @
1153 @ also
1154 @
1155 @ requires this.value() != 0 && other.value() == 0;
1156 @ assignable \nothing;
1157 @ ensures \result == other;
1158 @
1159 @ also
1160 @
1161 @ requires this.value() != 0 && other.value() != 0;
1162 @ assignable \nothing;
1163 @ ensures \result.value() == \old(this.value() * other.value());
1164 @*/
1165 private FDBigInteger mult(FDBigInteger other) {
1166 if (this.nWords == 0) {
1167 return this;
1168 }
1169 if (this.size() == 1) {
1170 return other.mult(data[0]);
1171 }
1172 if (other.nWords == 0) {
1173 return other;
1174 }
1175 if (other.size() == 1) {
1176 return this.mult(other.data[0]);
1177 }
1178 int[] r = new int[nWords + other.nWords];
1179 mult(this.data, this.nWords, other.data, other.nWords, r);
1180 return new FDBigInteger(r, this.offset + other.offset);
1181 }
1182
1183 /**
1184 * Adds another <code>FDBigInteger</code> to this <code>FDBigInteger</code>.
1185 *
1186 * @param other The <code>FDBigInteger</code> to add.
1187 * @return The sum of the <code>FDBigInteger</code>s.
1188 */
1189 /*@
1190 @ assignable \nothing;
1191 @ ensures \result.value() == \old(this.value() + other.value());
1192 @*/
1193 private FDBigInteger add(FDBigInteger other) {
1194 FDBigInteger big, small;
1195 int bigLen, smallLen;
1196 int tSize = this.size();
1197 int oSize = other.size();
1198 if (tSize >= oSize) {
1199 big = this;
1200 bigLen = tSize;
1201 small = other;
1202 smallLen = oSize;
1203 } else {
1204 big = other;
1205 bigLen = oSize;
1206 small = this;
1207 smallLen = tSize;
1208 }
1209 int[] r = new int[bigLen + 1];
1210 int i = 0;
1211 long carry = 0L;
1212 for (; i < smallLen; i++) {
1213 carry += (i < big.offset ? 0L : (big.data[i - big.offset] & LONG_MASK) )
1214 + ((i < small.offset ? 0L : (small.data[i - small.offset] & LONG_MASK)));
1215 r[i] = (int) carry;
1216 carry >>= 32; // signed shift.
1217 }
1218 for (; i < bigLen; i++) {
1219 carry += (i < big.offset ? 0L : (big.data[i - big.offset] & LONG_MASK) );
1220 r[i] = (int) carry;
1221 carry >>= 32; // signed shift.
1222 }
1223 r[bigLen] = (int) carry;
1224 return new FDBigInteger(r, 0);
1225 }
1226
1227
1228 /**
1229 * Multiplies a <code>FDBigInteger</code> by an int and adds another int. The
1230 * result is computed in place. This method is intended only to be invoked
1231 * from
1232 * <code>
1233 * FDBigInteger(long lValue, char[] digits, int kDigits, int nDigits)
1234 * </code>.
1235 *
1236 * @param iv The factor by which to multiply this <code>FDBigInteger</code>.
1237 * @param addend The value to add to the product of this
1238 * <code>FDBigInteger</code> and <code>iv</code>.
1239 */
1240 /*@
1241 @ requires this.value()*UNSIGNED(iv) + UNSIGNED(addend) < ((\bigint)1) << ((this.data.length + this.offset)*32);
1242 @ assignable this.data[*];
1243 @ ensures this.value() == \old(this.value()*UNSIGNED(iv) + UNSIGNED(addend));
1244 @*/
1245 private /*@ helper @*/ void multAddMe(int iv, int addend) {
1246 long v = iv & LONG_MASK;
1247 // unroll 0th iteration, doing addition.
1248 long p = v * (data[0] & LONG_MASK) + (addend & LONG_MASK);
1249 data[0] = (int) p;
1250 p >>>= 32;
1251 for (int i = 1; i < nWords; i++) {
1252 p += v * (data[i] & LONG_MASK);
1253 data[i] = (int) p;
1254 p >>>= 32;
1255 }
1256 if (p != 0L) {
1257 data[nWords++] = (int) p; // will fail noisily if illegal!
1258 }
1259 }
1260
1261 //
1262 // original doc:
1263 //
1264 // do this -=q*S
1265 // returns borrow
1266 //
1267 /**
1268 * Multiplies the parameters and subtracts them from this
1269 * <code>FDBigInteger</code>.
1270 *
1271 * @param q The integer parameter.
1272 * @param S The <code>FDBigInteger</code> parameter.
1273 * @return <code>this - q*S</code>.
1274 */
1275 /*@
1276 @ ensures nWords == 0 ==> offset == 0;
1277 @ ensures nWords > 0 ==> data[nWords - 1] != 0;
1278 @*/
1279 /*@
1280 @ requires 0 < q && q <= (1L << 31);
1281 @ requires data != null;
1282 @ requires 0 <= nWords && nWords <= data.length && offset >= 0;
1283 @ requires !this.isImmutable;
1284 @ requires this.size() == S.size();
1285 @ requires this != S;
1286 @ assignable this.nWords, this.offset, this.data, this.data[*];
1287 @ ensures -q <= \result && \result <= 0;
1288 @ ensures this.size() == \old(this.size());
1289 @ ensures this.value() + (\result << (this.size()*32)) == \old(this.value() - q*S.value());
1290 @ ensures this.offset == \old(Math.min(this.offset, S.offset));
1291 @ ensures \old(this.offset <= S.offset) ==> this.nWords == \old(this.nWords);
1292 @ ensures \old(this.offset <= S.offset) ==> this.offset == \old(this.offset);
1293 @ ensures \old(this.offset <= S.offset) ==> this.data == \old(this.data);
1294 @
1295 @ also
1296 @
1297 @ requires q == 0;
1298 @ assignable \nothing;
1299 @ ensures \result == 0;
1300 @*/
1301 private /*@ helper @*/ long multDiffMe(long q, FDBigInteger S) {
1302 long diff = 0L;
1303 if (q != 0) {
1304 int deltaSize = S.offset - this.offset;
1305 if (deltaSize >= 0) {
1306 int[] sd = S.data;
1307 int[] td = this.data;
1308 for (int sIndex = 0, tIndex = deltaSize; sIndex < S.nWords; sIndex++, tIndex++) {
1309 diff += (td[tIndex] & LONG_MASK) - q * (sd[sIndex] & LONG_MASK);
1310 td[tIndex] = (int) diff;
1311 diff >>= 32; // N.B. SIGNED shift.
1312 }
1313 } else {
1314 deltaSize = -deltaSize;
1315 int[] rd = new int[nWords + deltaSize];
1316 int sIndex = 0;
1317 int rIndex = 0;
1318 int[] sd = S.data;
1319 for (; rIndex < deltaSize && sIndex < S.nWords; sIndex++, rIndex++) {
1320 diff -= q * (sd[sIndex] & LONG_MASK);
1321 rd[rIndex] = (int) diff;
1322 diff >>= 32; // N.B. SIGNED shift.
1323 }
1324 int tIndex = 0;
1325 int[] td = this.data;
1326 for (; sIndex < S.nWords; sIndex++, tIndex++, rIndex++) {
1327 diff += (td[tIndex] & LONG_MASK) - q * (sd[sIndex] & LONG_MASK);
1328 rd[rIndex] = (int) diff;
1329 diff >>= 32; // N.B. SIGNED shift.
1330 }
1331 this.nWords += deltaSize;
1332 this.offset -= deltaSize;
1333 this.data = rd;
1334 }
1335 }
1336 return diff;
1337 }
1338
1339
1340 /**
1341 * Multiplies by 10 a big integer represented as an array. The final carry
1342 * is returned.
1343 *
1344 * @param src The array representation of the big integer.
1345 * @param srcLen The number of elements of <code>src</code> to use.
1346 * @param dst The product array.
1347 * @return The final carry of the multiplication.
1348 */
1349 /*@
1350 @ requires src.length >= srcLen && dst.length >= srcLen;
1351 @ assignable dst[0 .. srcLen - 1];
1352 @ ensures 0 <= \result && \result < 10;
1353 @ ensures AP(dst, srcLen) + (\result << (srcLen*32)) == \old(AP(src, srcLen) * 10);
1354 @*/
1355 private static int multAndCarryBy10(int[] src, int srcLen, int[] dst) {
1356 long carry = 0;
1357 for (int i = 0; i < srcLen; i++) {
1358 long product = (src[i] & LONG_MASK) * 10L + carry;
1359 dst[i] = (int) product;
1360 carry = product >>> 32;
1361 }
1362 return (int) carry;
1363 }
1364
1365 /**
1366 * Multiplies by a constant value a big integer represented as an array.
1367 * The constant factor is an <code>int</code>.
1368 *
1369 * @param src The array representation of the big integer.
1370 * @param srcLen The number of elements of <code>src</code> to use.
1371 * @param value The constant factor by which to multiply.
1372 * @param dst The product array.
1373 */
1374 /*@
1375 @ requires src.length >= srcLen && dst.length >= srcLen + 1;
1376 @ assignable dst[0 .. srcLen];
1377 @ ensures AP(dst, srcLen + 1) == \old(AP(src, srcLen) * UNSIGNED(value));
1378 @*/
1379 private static void mult(int[] src, int srcLen, int value, int[] dst) {
1380 long val = value & LONG_MASK;
1381 long carry = 0;
1382 for (int i = 0; i < srcLen; i++) {
1383 long product = (src[i] & LONG_MASK) * val + carry;
1384 dst[i] = (int) product;
1385 carry = product >>> 32;
1386 }
1387 dst[srcLen] = (int) carry;
1388 }
1389
1390 /**
1391 * Multiplies by a constant value a big integer represented as an array.
1392 * The constant factor is a long represent as two <code>int</code>s.
1393 *
1394 * @param src The array representation of the big integer.
1395 * @param srcLen The number of elements of <code>src</code> to use.
1396 * @param v0 The lower 32 bits of the long factor.
1397 * @param v1 The upper 32 bits of the long factor.
1398 * @param dst The product array.
1399 */
1400 /*@
1401 @ requires src != dst;
1402 @ requires src.length >= srcLen && dst.length >= srcLen + 2;
1403 @ assignable dst[0 .. srcLen + 1];
1404 @ ensures AP(dst, srcLen + 2) == \old(AP(src, srcLen) * (UNSIGNED(v0) + (UNSIGNED(v1) << 32)));
1405 @*/
1406 private static void mult(int[] src, int srcLen, int v0, int v1, int[] dst) {
1407 long v = v0 & LONG_MASK;
1408 long carry = 0;
1409 for (int j = 0; j < srcLen; j++) {
1410 long product = v * (src[j] & LONG_MASK) + carry;
1411 dst[j] = (int) product;
1412 carry = product >>> 32;
1413 }
1414 dst[srcLen] = (int) carry;
1415 v = v1 & LONG_MASK;
1416 carry = 0;
1417 for (int j = 0; j < srcLen; j++) {
1418 long product = (dst[j + 1] & LONG_MASK) + v * (src[j] & LONG_MASK) + carry;
1419 dst[j + 1] = (int) product;
1420 carry = product >>> 32;
1421 }
1422 dst[srcLen + 1] = (int) carry;
1423 }
1424
1425 // Fails assertion for negative exponent.
1426 /**
1427 * Computes <code>5</code> raised to a given power.
1428 *
1429 * @param p The exponent of 5.
1430 * @return <code>5<sup>p</sup></code>.
1431 */
1432 private static FDBigInteger big5pow(int p) {
1433 assert p >= 0 : p; // negative power of 5
1434 if (p < MAX_FIVE_POW) {
1435 return POW_5_CACHE[p];
1436 }
1437 return big5powRec(p);
1438 }
1439
1440 // slow path
1441 /**
1442 * Computes <code>5</code> raised to a given power.
1443 *
1444 * @param p The exponent of 5.
1445 * @return <code>5<sup>p</sup></code>.
1446 */
1447 private static FDBigInteger big5powRec(int p) {
1448 if (p < MAX_FIVE_POW) {
1449 return POW_5_CACHE[p];
1450 }
1451 // construct the value.
1452 // recursively.
1453 int q, r;
1454 // in order to compute 5^p,
1455 // compute its square root, 5^(p/2) and square.
1456 // or, let q = p / 2, r = p -q, then
1457 // 5^p = 5^(q+r) = 5^q * 5^r
1458 q = p >> 1;
1459 r = p - q;
1460 FDBigInteger bigq = big5powRec(q);
1461 if (r < SMALL_5_POW.length) {
1462 return bigq.mult(SMALL_5_POW[r]);
1463 } else {
1464 return bigq.mult(big5powRec(r));
1465 }
1466 }
1467
1468 // for debugging ...
1469 /**
1470 * Converts this <code>FDBigInteger</code> to a hexadecimal string.
1471 *
1472 * @return The hexadecimal string representation.
1473 */
1474 public String toHexString(){
1475 if(nWords ==0) {
1476 return "0";
1477 }
1478 StringBuilder sb = new StringBuilder((nWords +offset)*8);
1479 for(int i= nWords -1; i>=0; i--) {
1480 String subStr = Integer.toHexString(data[i]);
1481 for(int j = subStr.length(); j<8; j++) {
1482 sb.append('0');
1483 }
1484 sb.append(subStr);
1485 }
1486 for(int i=offset; i>0; i--) {
1487 sb.append("00000000");
1488 }
1489 return sb.toString();
1490 }
1491
1492 // for debugging ...
1493 /**
1494 * Converts this <code>FDBigInteger</code> to a <code>BigInteger</code>.
1495 *
1496 * @return The <code>BigInteger</code> representation.
1497 */
1498 public BigInteger toBigInteger() {
1499 byte[] magnitude = new byte[nWords * 4 + 1];
1500 for (int i = 0; i < nWords; i++) {
1501 int w = data[i];
1502 magnitude[magnitude.length - 4 * i - 1] = (byte) w;
1503 magnitude[magnitude.length - 4 * i - 2] = (byte) (w >> 8);
1504 magnitude[magnitude.length - 4 * i - 3] = (byte) (w >> 16);
1505 magnitude[magnitude.length - 4 * i - 4] = (byte) (w >> 24);
1506 }
1507 return new BigInteger(magnitude).shiftLeft(offset * 32);
1508 }
1509
1510 // for debugging ...
1511 /**
1512 * Converts this <code>FDBigInteger</code> to a string.
1513 *
1514 * @return The string representation.
1515 */
1516 @Override
1517 public String toString(){
1518 return toBigInteger().toString();
1519 }
1520 }