1 /* 2 * Copyright (c) 2013, 2020, 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 java.math.BigInteger; 28 import java.util.Arrays; 29 //@ model import org.jmlspecs.models.JMLMath; 30 31 /** 32 * A simple big integer package specifically for floating point base conversion. 33 */ 34 public /*@ spec_bigint_math @*/ class FDBigInteger { 35 36 // 37 // This class contains many comments that start with "/*@" mark. 38 // They are behavourial specification in 39 // the Java Modelling Language (JML): 40 // http://www.eecs.ucf.edu/~leavens/JML//index.shtml 41 // 42 43 /*@ 44 @ public pure model static \bigint UNSIGNED(int v) { 45 @ return v >= 0 ? v : v + (((\bigint)1) << 32); 46 @ } 47 @ 48 @ public pure model static \bigint UNSIGNED(long v) { 49 @ return v >= 0 ? v : v + (((\bigint)1) << 64); 50 @ } 51 @ 52 @ public pure model static \bigint AP(int[] data, int len) { 53 @ return (\sum int i; 0 <= 0 && i < len; UNSIGNED(data[i]) << (i*32)); 54 @ } 55 @ 56 @ public pure model static \bigint pow52(int p5, int p2) { 57 @ ghost \bigint v = 1; 58 @ for (int i = 0; i < p5; i++) v *= 5; 59 @ return v << p2; 60 @ } 61 @ 62 @ public pure model static \bigint pow10(int p10) { 63 @ return pow52(p10, p10); 64 @ } 65 @*/ 66 67 static final int[] SMALL_5_POW; 68 69 static final long[] LONG_5_POW; 70 71 // Maximum size of cache of powers of 5 as FDBigIntegers. 72 private static final int MAX_FIVE_POW = 340; 73 74 // Cache of big powers of 5 as FDBigIntegers. 75 private static final FDBigInteger POW_5_CACHE[]; 76 77 // Zero as an FDBigInteger. 78 public static final FDBigInteger ZERO; 79 80 // Archive proxy 81 private static Object[] archivedCaches; 82 83 // Initialize FDBigInteger cache of powers of 5. 84 static { 85 // Android-removed: CDS is not available. 86 // CDS.initializeFromArchive(FDBigInteger.class); 87 Object[] caches = archivedCaches; 88 if (caches == null) { 89 long[] long5pow = { 90 1L, 91 5L, 92 5L * 5, 93 5L * 5 * 5, 94 5L * 5 * 5 * 5, 95 5L * 5 * 5 * 5 * 5, 96 5L * 5 * 5 * 5 * 5 * 5, 97 5L * 5 * 5 * 5 * 5 * 5 * 5, 98 5L * 5 * 5 * 5 * 5 * 5 * 5 * 5, 99 5L * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5, 100 5L * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5, 101 5L * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5, 102 5L * 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, 104 5L * 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, 106 5L * 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, 108 5L * 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, 110 5L * 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, 112 5L * 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, 114 5L * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5, 115 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, 116 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, 117 }; 118 int[] small5pow = { 119 1, 120 5, 121 5 * 5, 122 5 * 5 * 5, 123 5 * 5 * 5 * 5, 124 5 * 5 * 5 * 5 * 5, 125 5 * 5 * 5 * 5 * 5 * 5, 126 5 * 5 * 5 * 5 * 5 * 5 * 5, 127 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5, 128 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5, 129 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5, 130 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5, 131 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5, 132 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 * 5 133 }; 134 FDBigInteger[] pow5cache = new FDBigInteger[MAX_FIVE_POW]; 135 int i = 0; 136 while (i < small5pow.length) { 137 FDBigInteger pow5 = new FDBigInteger(new int[] { small5pow[i] }, 0); pow5.makeImmutable()138 pow5.makeImmutable(); 139 pow5cache[i] = pow5; 140 i++; 141 } 142 FDBigInteger prev = pow5cache[i - 1]; 143 while (i < MAX_FIVE_POW) { 144 pow5cache[i] = prev = prev.mult(5); prev.makeImmutable()145 prev.makeImmutable(); 146 i++; 147 } 148 FDBigInteger zero = new FDBigInteger(new int[0], 0); zero.makeImmutable()149 zero.makeImmutable(); 150 archivedCaches = caches = new Object[] {small5pow, long5pow, pow5cache, zero}; 151 } 152 SMALL_5_POW = (int[])caches[0]; 153 LONG_5_POW = (long[])caches[1]; 154 POW_5_CACHE = (FDBigInteger[])caches[2]; 155 ZERO = (FDBigInteger)caches[3]; 156 } 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 @*/ FDBigInteger(int[] data, int offset)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 @*/ FDBigInteger(long lValue, char[] digits, int kDigits, int nDigits)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 @*/ valueOfPow52(int p5, int p2)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 @*/ valueOfMulPow52(long value, int p5, int p2)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 @*/ valueOfPow2(int p2)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 @*/ trimLeadingZeros()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 @*/ getNormalizationBias()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 @*/ leftShift(int[] src, int idx, int result[], int bitcount, int anticount, int prev)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 @*/ leftShift(int shift)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 @*/ size()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 @*/ quoRemIteration(FDBigInteger S)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 @*/ multBy10()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 @*/ multByPow52(int p5, int p2)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 @*/ mult(int[] s1, int s1Len, int[] s2, int s2Len, int[] dst)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 @*/ leftInplaceSub(FDBigInteger subtrahend)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 @*/ rightInplaceSub(FDBigInteger subtrahend)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 @*/ checkZeroTail(int[] a, int from)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 @*/ cmp(FDBigInteger other)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 @*/ cmpPow52(int p5, int p2)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 @*/ addAndCmp(FDBigInteger x, FDBigInteger y)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 @*/ makeImmutable()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 @*/ mult(int i)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 @*/ mult(FDBigInteger other)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 @*/ add(FDBigInteger other)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 @*/ multAddMe(int iv, int addend)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 @*/ multDiffMe(long q, FDBigInteger S)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 @*/ multAndCarryBy10(int[] src, int srcLen, int[] dst)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 @*/ mult(int[] src, int srcLen, int value, int[] dst)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 @*/ mult(int[] src, int srcLen, int v0, int v1, int[] dst)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 */ big5pow(int p)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 */ big5powRec(int p)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 */ toHexString()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 */ toBigInteger()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 toString()1517 public String toString(){ 1518 return toBigInteger().toString(); 1519 } 1520 } 1521