1244
1245 // Clear the quotient
1246 quotient.offset = quotient.intLen = 0;
1247
1248 if (r < s) {
1249 return this;
1250 } else {
1251 // Unlike Knuth division, we don't check for common powers of two here because
1252 // BZ already runs faster if both numbers contain powers of two and cancelling them has no
1253 // additional benefit.
1254
1255 // step 1: let m = min{2^k | (2^k)*BURNIKEL_ZIEGLER_THRESHOLD > s}
1256 int m = 1 << (32-Integer.numberOfLeadingZeros(s/BigInteger.BURNIKEL_ZIEGLER_THRESHOLD));
1257
1258 int j = (s+m-1) / m; // step 2a: j = ceil(s/m)
1259 int n = j * m; // step 2b: block length in 32-bit units
1260 long n32 = 32L * n; // block length in bits
1261 int sigma = (int) Math.max(0, n32 - b.bitLength()); // step 3: sigma = max{T | (2^T)*B < beta^n}
1262 MutableBigInteger bShifted = new MutableBigInteger(b);
1263 bShifted.safeLeftShift(sigma); // step 4a: shift b so its length is a multiple of n
1264 safeLeftShift(sigma); // step 4b: shift this by the same amount
1265
1266 // step 5: t is the number of blocks needed to accommodate this plus one additional bit
1267 int t = (int) ((bitLength()+n32) / n32);
1268 if (t < 2) {
1269 t = 2;
1270 }
1271
1272 // step 6: conceptually split this into blocks a[t-1], ..., a[0]
1273 MutableBigInteger a1 = getBlock(t-1, t, n); // the most significant block of this
1274
1275 // step 7: z[t-2] = [a[t-1], a[t-2]]
1276 MutableBigInteger z = getBlock(t-2, t, n); // the second to most significant block
1277 z.addDisjoint(a1, n); // z[t-2]
1278
1279 // do schoolbook division on blocks, dividing 2-block numbers by 1-block numbers
1280 MutableBigInteger qi = new MutableBigInteger();
1281 MutableBigInteger ri;
1282 for (int i=t-2; i > 0; i--) {
1283 // step 8a: compute (qi,ri) such that z=b*qi+ri
1284 ri = z.divide2n1n(bShifted, qi);
1285
1286 // step 8b: z = [ri, a[i-1]]
1287 z = getBlock(i-1, t, n); // a[i-1]
1288 z.addDisjoint(ri, n);
1289 quotient.addShifted(qi, i*n); // update q (part of step 9)
1290 }
1291 // final iteration of step 8: do the loop one more time for i=0 but leave z unchanged
1292 ri = z.divide2n1n(bShifted, qi);
1293 quotient.add(qi);
1294
1295 ri.rightShift(sigma); // step 9: this and b were shifted, so shift back
1296 return ri;
1297 }
1298 }
1299
1300 /**
1301 * This method implements algorithm 1 from pg. 4 of the Burnikel-Ziegler paper.
1302 * It divides a 2n-digit number by a n-digit number.<br/>
1303 * The parameter beta is 2<sup>32</sup> so all shifts are multiples of 32 bits.
1304 * <br/>
1305 * {@code this} must be a nonnegative number such that {@code this.bitLength() <= 2*b.bitLength()}
1306 * @param b a positive number such that {@code b.bitLength()} is even
1307 * @param quotient output parameter for {@code this/b}
1308 * @return {@code this%b}
1309 */
1310 private MutableBigInteger divide2n1n(MutableBigInteger b, MutableBigInteger quotient) {
1311 int n = b.intLen;
1312
1313 // step 1: base case
1314 if (n%2 != 0 || n < BigInteger.BURNIKEL_ZIEGLER_THRESHOLD) {
1315 return divideKnuth(b, quotient);
|
1244
1245 // Clear the quotient
1246 quotient.offset = quotient.intLen = 0;
1247
1248 if (r < s) {
1249 return this;
1250 } else {
1251 // Unlike Knuth division, we don't check for common powers of two here because
1252 // BZ already runs faster if both numbers contain powers of two and cancelling them has no
1253 // additional benefit.
1254
1255 // step 1: let m = min{2^k | (2^k)*BURNIKEL_ZIEGLER_THRESHOLD > s}
1256 int m = 1 << (32-Integer.numberOfLeadingZeros(s/BigInteger.BURNIKEL_ZIEGLER_THRESHOLD));
1257
1258 int j = (s+m-1) / m; // step 2a: j = ceil(s/m)
1259 int n = j * m; // step 2b: block length in 32-bit units
1260 long n32 = 32L * n; // block length in bits
1261 int sigma = (int) Math.max(0, n32 - b.bitLength()); // step 3: sigma = max{T | (2^T)*B < beta^n}
1262 MutableBigInteger bShifted = new MutableBigInteger(b);
1263 bShifted.safeLeftShift(sigma); // step 4a: shift b so its length is a multiple of n
1264 MutableBigInteger aShifted = new MutableBigInteger (this);
1265 aShifted.safeLeftShift(sigma); // step 4b: shift a by the same amount
1266
1267 // step 5: t is the number of blocks needed to accommodate a plus one additional bit
1268 int t = (int) ((aShifted.bitLength()+n32) / n32);
1269 if (t < 2) {
1270 t = 2;
1271 }
1272
1273 // step 6: conceptually split a into blocks a[t-1], ..., a[0]
1274 MutableBigInteger a1 = aShifted.getBlock(t-1, t, n); // the most significant block of a
1275
1276 // step 7: z[t-2] = [a[t-1], a[t-2]]
1277 MutableBigInteger z = aShifted.getBlock(t-2, t, n); // the second to most significant block
1278 z.addDisjoint(a1, n); // z[t-2]
1279
1280 // do schoolbook division on blocks, dividing 2-block numbers by 1-block numbers
1281 MutableBigInteger qi = new MutableBigInteger();
1282 MutableBigInteger ri;
1283 for (int i=t-2; i > 0; i--) {
1284 // step 8a: compute (qi,ri) such that z=b*qi+ri
1285 ri = z.divide2n1n(bShifted, qi);
1286
1287 // step 8b: z = [ri, a[i-1]]
1288 z = aShifted.getBlock(i-1, t, n); // a[i-1]
1289 z.addDisjoint(ri, n);
1290 quotient.addShifted(qi, i*n); // update q (part of step 9)
1291 }
1292 // final iteration of step 8: do the loop one more time for i=0 but leave z unchanged
1293 ri = z.divide2n1n(bShifted, qi);
1294 quotient.add(qi);
1295
1296 ri.rightShift(sigma); // step 9: a and b were shifted, so shift back
1297 return ri;
1298 }
1299 }
1300
1301 /**
1302 * This method implements algorithm 1 from pg. 4 of the Burnikel-Ziegler paper.
1303 * It divides a 2n-digit number by a n-digit number.<br/>
1304 * The parameter beta is 2<sup>32</sup> so all shifts are multiples of 32 bits.
1305 * <br/>
1306 * {@code this} must be a nonnegative number such that {@code this.bitLength() <= 2*b.bitLength()}
1307 * @param b a positive number such that {@code b.bitLength()} is even
1308 * @param quotient output parameter for {@code this/b}
1309 * @return {@code this%b}
1310 */
1311 private MutableBigInteger divide2n1n(MutableBigInteger b, MutableBigInteger quotient) {
1312 int n = b.intLen;
1313
1314 // step 1: base case
1315 if (n%2 != 0 || n < BigInteger.BURNIKEL_ZIEGLER_THRESHOLD) {
1316 return divideKnuth(b, quotient);
|