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 static final 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}, but the next bit should be 1

   388      * {@code (highestWord & 0x08000000) != 0}.

   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, {@code 0 <= q < 10}.

   550      *

   551      * @param S 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 The multiplication result.

   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>{@code

   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>{@code

   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>{@code

  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.