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