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 }