|
|
@@ -1,3200 +0,0 @@
|
|
|
-#include <vector>
|
|
|
-#include <algorithm>
|
|
|
-#include <cstring>
|
|
|
-#include <random>
|
|
|
-#include <iterator>
|
|
|
-#include <stdexcept>
|
|
|
-#include <cstdint>
|
|
|
-#include <cassert>
|
|
|
-using std::uint64_t;
|
|
|
-static uint64_t bitsPerDigit[] = { 0, 0,
|
|
|
- 1024, 1624, 2048, 2378, 2648, 2875, 3072, 3247, 3402, 3543, 3672,
|
|
|
- 3790, 3899, 4001, 4096, 4186, 4271, 4350, 4426, 4498, 4567, 4633,
|
|
|
- 4696, 4756, 4814, 4870, 4923, 4975, 5025, 5074, 5120, 5166, 5210,
|
|
|
- 5253, 5295};
|
|
|
-class BigInteger{
|
|
|
- int signum;
|
|
|
- std::vector<int> mag;
|
|
|
- int bitCount;
|
|
|
- int bitLength;
|
|
|
- int lowestSetBit;
|
|
|
- int firstNonzeroIntNum;
|
|
|
- const static uint64_t LONG_MASK = 0xffffffffL;
|
|
|
- static const int MAX_MAG_LENGTH = (1 << 26);
|
|
|
- static const int PRIME_SEARCH_BIT_LENGTH_LIMIT = 500000000;
|
|
|
- static const int KARATSUBA_THRESHOLD = 80;
|
|
|
- static const int TOOM_COOK_THRESHOLD = 240;
|
|
|
- static const int KARATSUBA_SQUARE_THRESHOLD = 128;
|
|
|
- static const int TOOM_COOK_SQUARE_THRESHOLD = 216;
|
|
|
- static const int BURNIKEL_ZIEGLER_THRESHOLD = 80;
|
|
|
- static const int BURNIKEL_ZIEGLER_OFFSET = 40;
|
|
|
- static const int SCHOENHAGE_BASE_CONVERSION_THRESHOLD = 20;
|
|
|
- static const int MULTIPLY_SQUARE_THRESHOLD = 20;
|
|
|
- static const int MONTGOMERY_INTRINSIC_THRESHOLD = 512;
|
|
|
- BigInteger(std::vector<char> val) {
|
|
|
- assert(val.size() != 0);
|
|
|
- if (val[0] < 0) {
|
|
|
- mag = makePositive(val);
|
|
|
- signum = -1;
|
|
|
- } else {
|
|
|
- mag = stripLeadingZeroBytes(val);
|
|
|
- signum = (mag.size() == 0 ? 0 : 1);
|
|
|
- }
|
|
|
- if (mag.size() >= MAX_MAG_LENGTH) {
|
|
|
- checkRange();
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- BigInteger(std::vector<int> val) {
|
|
|
- assert(val.size() != 0);
|
|
|
-
|
|
|
- if (val[0] < 0) {
|
|
|
- mag = makePositive(val);
|
|
|
- signum = -1;
|
|
|
- } else {
|
|
|
- mag = trustedStripLeadingZeroInts(val);
|
|
|
- signum = (mag.size() == 0 ? 0 : 1);
|
|
|
- }
|
|
|
- if (mag.size() >= MAX_MAG_LENGTH) {
|
|
|
- checkRange();
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- BigInteger(int signum, std::vector<char> magnitude) {
|
|
|
- mag = stripLeadingZeroBytes(magnitude);
|
|
|
-
|
|
|
- assert(!(signum < -1 || signum > 1));
|
|
|
-
|
|
|
- if (mag.size() == 0) {
|
|
|
- signum = 0;
|
|
|
- } else {
|
|
|
- assert(signum != 0);
|
|
|
- signum = signum;
|
|
|
- }
|
|
|
- if (mag.size() >= MAX_MAG_LENGTH) {
|
|
|
- checkRange();
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
- BigInteger():BigInteger((int)0){}
|
|
|
- BigInteger(int signum, std::vector<int> magnitude) {
|
|
|
- mag = stripLeadingZeroInts(magnitude);
|
|
|
-
|
|
|
- assert(!(signum < -1 || signum > 1));
|
|
|
-
|
|
|
- if (this.mag.size() == 0) {
|
|
|
- signum = 0;
|
|
|
- } else {
|
|
|
- assert(signum != 0);
|
|
|
- signum = signum;
|
|
|
- }
|
|
|
- if (mag.size() >= MAX_MAG_LENGTH) {
|
|
|
- checkRange();
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
-
|
|
|
- /*
|
|
|
- * Constructs a new BigInteger using a char array with radix=10.
|
|
|
- * Sign is precalculated outside and not allowed in the val.
|
|
|
- */
|
|
|
- BigInteger(std::vector<char> val, int sign, int len) {
|
|
|
- int cursor = 0, numDigits;
|
|
|
-
|
|
|
- // Skip leading zeros and compute number of digits in magnitude
|
|
|
- while (cursor < len && Character.digit(val[cursor], 10) == 0) {
|
|
|
- cursor++;
|
|
|
- }
|
|
|
- if (cursor == len) {
|
|
|
- signum = 0;
|
|
|
- mag = ZERO.mag;
|
|
|
- return;
|
|
|
- }
|
|
|
-
|
|
|
- numDigits = len - cursor;
|
|
|
- signum = sign;
|
|
|
- // Pre-allocate array of expected size
|
|
|
- unsigned int numWords;
|
|
|
- if (len < 10) {
|
|
|
- numWords = 1;
|
|
|
- } else {
|
|
|
- uint64_t numBits = ((numDigits * bitsPerDigit[10]) >>> 10) + 1;
|
|
|
- if (numBits + 31 >= (1L << 32)) {
|
|
|
- reportOverflow();
|
|
|
- }
|
|
|
- numWords = (int) (numBits + 31) >>> 5;
|
|
|
- }
|
|
|
- std::vector<int> magnitude(numBits);
|
|
|
-
|
|
|
- // Process first (potentially short) digit group
|
|
|
- int firstGroupLen = numDigits % digitsPerInt[10];
|
|
|
- if (firstGroupLen == 0)
|
|
|
- firstGroupLen = digitsPerInt[10];
|
|
|
- magnitude[numWords - 1] = parseInt(val, cursor, cursor += firstGroupLen);
|
|
|
-
|
|
|
- // Process remaining digit groups
|
|
|
- while (cursor < len) {
|
|
|
- int groupVal = parseInt(val, cursor, cursor += digitsPerInt[10]);
|
|
|
- destructiveMulAdd(magnitude, intRadix[10], groupVal);
|
|
|
- }
|
|
|
- mag = trustedStripLeadingZeroInts(magnitude);
|
|
|
- if (mag.size() >= MAX_MAG_LENGTH) {
|
|
|
- checkRange();
|
|
|
- }
|
|
|
- }
|
|
|
- int digit(char a){
|
|
|
- assert((int)(a - '0') < 10 && (int)(a - '0') >= 0);
|
|
|
- return (int)(a - '0');
|
|
|
- }
|
|
|
- // Create an integer with the digits between the two indexes
|
|
|
- // Assumes start < end. The result may be negative, but it
|
|
|
- // is to be treated as an unsigned value.
|
|
|
- int parseInt(const std::vector<char>& source, int start, int end) {
|
|
|
- int result = digit(source[start++]);
|
|
|
-
|
|
|
- for (int index = start; index < end; index++) {
|
|
|
- int nextVal = digit(source[index]);
|
|
|
- result = 10 * result + nextVal;
|
|
|
- }
|
|
|
-
|
|
|
- return result;
|
|
|
- }
|
|
|
-
|
|
|
- // bitsPerDigit in the given radix times 1024
|
|
|
- // Rounded up to avoid underallocation.
|
|
|
-
|
|
|
-
|
|
|
- // Multiply x array times word y in place, and add word z
|
|
|
- static void destructiveMulAdd(std::vector<int>& x, int y, int z) {
|
|
|
- // Perform the multiplication word by word
|
|
|
- uint64_t ylong = y & LONG_MASK;
|
|
|
- uint64_t zlong = z & LONG_MASK;
|
|
|
- int len = x.size();
|
|
|
-
|
|
|
- uint64_t product = 0;
|
|
|
- uint64_t carry = 0;
|
|
|
- for (int i = len-1; i >= 0; i--) {
|
|
|
- product = ylong * (x[i] & LONG_MASK) + carry;
|
|
|
- x[i] = (int)product;
|
|
|
- carry = product >>> 32;
|
|
|
- }
|
|
|
-
|
|
|
- // Perform the addition
|
|
|
- uint64_t sum = (x[len-1] & LONG_MASK) + zlong;
|
|
|
- x[len-1] = (int)sum;
|
|
|
- carry = sum >> 32;
|
|
|
- for (int i = len-2; i >= 0; i--) {
|
|
|
- sum = (x[i] & LONG_MASK) + carry;
|
|
|
- x[i] = (int)sum;
|
|
|
- carry = sum >> 32;
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- BigInteger(std::string val) : BigInteger(val, 10) {
|
|
|
-
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- BigInteger(int numBits, std::mt19937_64& rnd) : BigInteger(1, randomBits(numBits, rnd)) {
|
|
|
-
|
|
|
- }
|
|
|
-
|
|
|
- static std::vector<char> randomBits(unsigned int numBits, std::mt19937_64& rnd) {
|
|
|
- unsigned int numBytes = (unsigned int)(((uint64_t)numBits+7)/8); // avoid overflow
|
|
|
- std::vector<char> randomBits(numBytes);
|
|
|
- std::uniform_int_distribution<char> dis(-128,127);
|
|
|
- // Generate random bytes and mask out any excess bits
|
|
|
- if (numBytes > 0) {
|
|
|
- std::generate(randomBits.begin(), randomBits.end(), [&](){return dis(rnd)});
|
|
|
- int excessBits = 8*numBytes - numBits;
|
|
|
- randomBits[0] &= (1 << (8-excessBits)) - 1;
|
|
|
- }
|
|
|
- return randomBits;
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- BigInteger(int bitLength, int certainty, std::mt19937_64& rnd) {
|
|
|
- BigInteger prime;
|
|
|
-
|
|
|
- assert(bitLength >= 2);
|
|
|
- prime = (bitLength < SMALL_PRIME_THRESHOLD
|
|
|
- ? smallPrime(bitLength, certainty, rnd)
|
|
|
- : largePrime(bitLength, certainty, rnd));
|
|
|
- signum = 1;
|
|
|
- mag = prime.mag;
|
|
|
- }
|
|
|
-
|
|
|
- // Minimum size in bits that the requested prime number has
|
|
|
- // before we use the large prime number generating algorithms.
|
|
|
- // The cutoff of 95 was chosen empirically for best performance.
|
|
|
- static const int SMALL_PRIME_THRESHOLD = 95;
|
|
|
-
|
|
|
- // Certainty required to meet the spec of probablePrime
|
|
|
- static const int DEFAULT_PRIME_CERTAINTY = 100;
|
|
|
-
|
|
|
-
|
|
|
- /*static BigInteger probablePrime(int bitLength, Random rnd) {
|
|
|
- if (bitLength < 2)
|
|
|
- throw new ArithmeticException("bitLength < 2");
|
|
|
-
|
|
|
- return (bitLength < SMALL_PRIME_THRESHOLD ?
|
|
|
- smallPrime(bitLength, DEFAULT_PRIME_CERTAINTY, rnd) :
|
|
|
- largePrime(bitLength, DEFAULT_PRIME_CERTAINTY, rnd));
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- static BigInteger smallPrime(int bitLength, int certainty, Random rnd) {
|
|
|
- int magLen = (bitLength + 31) >>> 5;
|
|
|
- int temp[] = new int[magLen];
|
|
|
- int highBit = 1 << ((bitLength+31) & 0x1f); // High bit of high int
|
|
|
- int highMask = (highBit << 1) - 1; // Bits to keep in high int
|
|
|
-
|
|
|
- while (true) {
|
|
|
- // Construct a candidate
|
|
|
- for (int i=0; i < magLen; i++)
|
|
|
- temp[i] = rnd.nextInt();
|
|
|
- temp[0] = (temp[0] & highMask) | highBit; // Ensure exact length
|
|
|
- if (bitLength > 2)
|
|
|
- temp[magLen-1] |= 1; // Make odd if bitlen > 2
|
|
|
-
|
|
|
- BigInteger p = new BigInteger(temp, 1);
|
|
|
-
|
|
|
- // Do cheap "pre-test" if applicable
|
|
|
- if (bitLength > 6) {
|
|
|
- long r = p.remainder(SMALL_PRIME_PRODUCT).longValue();
|
|
|
- if ((r%3==0) || (r%5==0) || (r%7==0) || (r%11==0) ||
|
|
|
- (r%13==0) || (r%17==0) || (r%19==0) || (r%23==0) ||
|
|
|
- (r%29==0) || (r%31==0) || (r%37==0) || (r%41==0))
|
|
|
- continue; // Candidate is composite; try another
|
|
|
- }
|
|
|
-
|
|
|
- // All candidates of bitLength 2 and 3 are prime by this point
|
|
|
- if (bitLength < 4)
|
|
|
- return p;
|
|
|
-
|
|
|
- // Do expensive test if we survive pre-test (or it's inapplicable)
|
|
|
- if (p.primeToCertainty(certainty, rnd))
|
|
|
- return p;
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
- static const BigInteger SMALL_PRIME_PRODUCT
|
|
|
- = valueOf(3L*5*7*11*13*17*19*23*29*31*37*41);
|
|
|
-
|
|
|
-
|
|
|
- static BigInteger largePrime(int bitLength, int certainty, Random rnd) {
|
|
|
- BigInteger p;
|
|
|
- p = new BigInteger(bitLength, rnd).setBit(bitLength-1);
|
|
|
- p.mag[p.mag.length-1] &= 0xfffffffe;
|
|
|
-
|
|
|
- // Use a sieve length likely to contain the next prime number
|
|
|
- int searchLen = getPrimeSearchLen(bitLength);
|
|
|
- BitSieve searchSieve = new BitSieve(p, searchLen);
|
|
|
- BigInteger candidate = searchSieve.retrieve(p, certainty, rnd);
|
|
|
-
|
|
|
- while ((candidate == null) || (candidate.bitLength() != bitLength)) {
|
|
|
- p = p.add(BigInteger.valueOf(2*searchLen));
|
|
|
- if (p.bitLength() != bitLength)
|
|
|
- p = new BigInteger(bitLength, rnd).setBit(bitLength-1);
|
|
|
- p.mag[p.mag.length-1] &= 0xfffffffe;
|
|
|
- searchSieve = new BitSieve(p, searchLen);
|
|
|
- candidate = searchSieve.retrieve(p, certainty, rnd);
|
|
|
- }
|
|
|
- return candidate;
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- BigInteger nextProbablePrime() {
|
|
|
- if (this.signum < 0)
|
|
|
- throw new ArithmeticException("start < 0: " + this);
|
|
|
-
|
|
|
- // Handle trivial cases
|
|
|
- if ((this.signum == 0) || this.equals(ONE))
|
|
|
- return TWO;
|
|
|
-
|
|
|
- BigInteger result = this.add(ONE);
|
|
|
-
|
|
|
- // Fastpath for small numbers
|
|
|
- if (result.bitLength() < SMALL_PRIME_THRESHOLD) {
|
|
|
-
|
|
|
- // Ensure an odd number
|
|
|
- if (!result.testBit(0))
|
|
|
- result = result.add(ONE);
|
|
|
-
|
|
|
- while (true) {
|
|
|
- // Do cheap "pre-test" if applicable
|
|
|
- if (result.bitLength() > 6) {
|
|
|
- long r = result.remainder(SMALL_PRIME_PRODUCT).longValue();
|
|
|
- if ((r%3==0) || (r%5==0) || (r%7==0) || (r%11==0) ||
|
|
|
- (r%13==0) || (r%17==0) || (r%19==0) || (r%23==0) ||
|
|
|
- (r%29==0) || (r%31==0) || (r%37==0) || (r%41==0)) {
|
|
|
- result = result.add(TWO);
|
|
|
- continue; // Candidate is composite; try another
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
- // All candidates of bitLength 2 and 3 are prime by this point
|
|
|
- if (result.bitLength() < 4)
|
|
|
- return result;
|
|
|
-
|
|
|
- // The expensive test
|
|
|
- if (result.primeToCertainty(DEFAULT_PRIME_CERTAINTY, null))
|
|
|
- return result;
|
|
|
-
|
|
|
- result = result.add(TWO);
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
- // Start at previous even number
|
|
|
- if (result.testBit(0))
|
|
|
- result = result.subtract(ONE);
|
|
|
-
|
|
|
- // Looking for the next large prime
|
|
|
- int searchLen = getPrimeSearchLen(result.bitLength());
|
|
|
-
|
|
|
- while (true) {
|
|
|
- BitSieve searchSieve = new BitSieve(result, searchLen);
|
|
|
- BigInteger candidate = searchSieve.retrieve(result,
|
|
|
- DEFAULT_PRIME_CERTAINTY, null);
|
|
|
- if (candidate != null)
|
|
|
- return candidate;
|
|
|
- result = result.add(BigInteger.valueOf(2 * searchLen));
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
- static int getPrimeSearchLen(int bitLength) {
|
|
|
- if (bitLength > PRIME_SEARCH_BIT_LENGTH_LIMIT + 1) {
|
|
|
- throw new ArithmeticException("Prime search implementation restriction on bitLength");
|
|
|
- }
|
|
|
- return bitLength / 20 * 64;
|
|
|
- }*/
|
|
|
-
|
|
|
-
|
|
|
- bool primeToCertainty(int certainty, std::mt19937_64& random) {
|
|
|
- int rounds = 0;
|
|
|
- int n = (std::min(certainty, Integer.MAX_VALUE-1)+1)/2;
|
|
|
-
|
|
|
- // The relationship between the certainty and the number of rounds
|
|
|
- // we perform is given in the draft standard ANSI X9.80, "PRIME
|
|
|
- // NUMBER GENERATION, PRIMALITY TESTING, AND PRIMALITY CERTIFICATES".
|
|
|
- int sizeInBits = bitLength();
|
|
|
- if (sizeInBits < 100) {
|
|
|
- rounds = 50;
|
|
|
- rounds = n < rounds ? n : rounds;
|
|
|
- return passesMillerRabin(rounds, random);
|
|
|
- }
|
|
|
-
|
|
|
- if (sizeInBits < 256) {
|
|
|
- rounds = 27;
|
|
|
- } else if (sizeInBits < 512) {
|
|
|
- rounds = 15;
|
|
|
- } else if (sizeInBits < 768) {
|
|
|
- rounds = 8;
|
|
|
- } else if (sizeInBits < 1024) {
|
|
|
- rounds = 4;
|
|
|
- } else {
|
|
|
- rounds = 2;
|
|
|
- }
|
|
|
- rounds = n < rounds ? n : rounds;
|
|
|
-
|
|
|
- return passesMillerRabin(rounds, random) && passesLucasLehmer();
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- bool passesLucasLehmer() {
|
|
|
- BigInteger thisPlusOne = this->add(ONE);
|
|
|
-
|
|
|
- // Step 1
|
|
|
- int d = 5;
|
|
|
- while (jacobiSymbol(d, this) != -1) {
|
|
|
- // 5, -7, 9, -11, ...
|
|
|
- d = (d < 0) ? std::abs(d)+2 : -(d+2);
|
|
|
- }
|
|
|
-
|
|
|
- // Step 2
|
|
|
- BigInteger u = lucasLehmerSequence(d, thisPlusOne, *this);
|
|
|
-
|
|
|
- // Step 3
|
|
|
- return u.mod(*this).equals(ZERO);
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- static int jacobiSymbol(int p, BigInteger n) {
|
|
|
- if (p == 0)
|
|
|
- return 0;
|
|
|
-
|
|
|
- // Algorithm and comments adapted from Colin Plumb's C library.
|
|
|
- int j = 1;
|
|
|
- int u = n.mag[n.mag.size()-1];
|
|
|
-
|
|
|
- // Make p positive
|
|
|
- if (p < 0) {
|
|
|
- p = -p;
|
|
|
- int n8 = u & 7;
|
|
|
- if ((n8 == 3) || (n8 == 7))
|
|
|
- j = -j; // 3 (011) or 7 (111) mod 8
|
|
|
- }
|
|
|
-
|
|
|
- // Get rid of factors of 2 in p
|
|
|
- while ((p & 3) == 0)
|
|
|
- p >>= 2;
|
|
|
- if ((p & 1) == 0) {
|
|
|
- p >>= 1;
|
|
|
- if (((u ^ (u>>1)) & 2) != 0)
|
|
|
- j = -j; // 3 (011) or 5 (101) mod 8
|
|
|
- }
|
|
|
- if (p == 1)
|
|
|
- return j;
|
|
|
- // Then, apply quadratic reciprocity
|
|
|
- if ((p & u & 2) != 0) // p = u = 3 (mod 4)?
|
|
|
- j = -j;
|
|
|
- // And reduce u mod p
|
|
|
- u = n.mod(BigInteger::valueOf(p)).intValue();
|
|
|
-
|
|
|
- // Now compute Jacobi(u,p), u < p
|
|
|
- while (u != 0) {
|
|
|
- while ((u & 3) == 0)
|
|
|
- u >>= 2;
|
|
|
- if ((u & 1) == 0) {
|
|
|
- u >>= 1;
|
|
|
- if (((p ^ (p>>1)) & 2) != 0)
|
|
|
- j = -j; // 3 (011) or 5 (101) mod 8
|
|
|
- }
|
|
|
- if (u == 1)
|
|
|
- return j;
|
|
|
- // Now both u and p are odd, so use quadratic reciprocity
|
|
|
- assert (u < p);
|
|
|
- int t = u; u = p; p = t;
|
|
|
- if ((u & p & 2) != 0) // u = p = 3 (mod 4)?
|
|
|
- j = -j;
|
|
|
- // Now u >= p, so it can be reduced
|
|
|
- u %= p;
|
|
|
- }
|
|
|
- return 0;
|
|
|
- }
|
|
|
-
|
|
|
- static BigInteger lucasLehmerSequence(int z, BigInteger k, BigInteger n) {
|
|
|
- BigInteger d = BigInteger::valueOf(z);
|
|
|
- BigInteger u = ONE; BigInteger u2;
|
|
|
- BigInteger v = ONE; BigInteger v2;
|
|
|
-
|
|
|
- for (int i=k.bitLength()-2; i >= 0; i--) {
|
|
|
- u2 = u.multiply(v).mod(n);
|
|
|
-
|
|
|
- v2 = v.square().add(d.multiply(u.square())).mod(n);
|
|
|
- if (v2.testBit(0))
|
|
|
- v2 = v2.subtract(n);
|
|
|
-
|
|
|
- v2 = v2.shiftRight(1);
|
|
|
-
|
|
|
- u = u2; v = v2;
|
|
|
- if (k.testBit(i)) {
|
|
|
- u2 = u.add(v).mod(n);
|
|
|
- if (u2.testBit(0))
|
|
|
- u2 = u2.subtract(n);
|
|
|
-
|
|
|
- u2 = u2.shiftRight(1);
|
|
|
- v2 = v.add(d.multiply(u)).mod(n);
|
|
|
- if (v2.testBit(0))
|
|
|
- v2 = v2.subtract(n);
|
|
|
- v2 = v2.shiftRight(1);
|
|
|
-
|
|
|
- u = u2; v = v2;
|
|
|
- }
|
|
|
- }
|
|
|
- return u;
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- bool passesMillerRabin(int iterations, std::mt19937_64& rnd) {
|
|
|
- // Find a and m such that m is odd and this == 1 + 2**a * m
|
|
|
- BigInteger thisMinusOne = this->subtract(ONE);
|
|
|
- BigInteger m = thisMinusOne;
|
|
|
- int a = m.getLowestSetBit();
|
|
|
- m = m.shiftRight(a);
|
|
|
-
|
|
|
- // Do the tests
|
|
|
- for (int i=0; i < iterations; i++) {
|
|
|
- // Generate a uniform random on (1, this)
|
|
|
- BigInteger b;
|
|
|
- do {
|
|
|
- b = BigInteger(this->bitLength(), rnd);
|
|
|
- } while (b.compareTo(ONE) <= 0 || b.compareTo(*this) >= 0);
|
|
|
-
|
|
|
- int j = 0;
|
|
|
- BigInteger z = b.modPow(m, *this);
|
|
|
- while (!((j == 0 && z.equals(ONE)) || z.equals(thisMinusOne))) {
|
|
|
- if (j > 0 && z.equals(ONE) || ++j == a)
|
|
|
- return false;
|
|
|
- z = z.modPow(TWO, *this);
|
|
|
- }
|
|
|
- }
|
|
|
- return true;
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- BigInteger(const std::vector<int>& magnitude, int signum) {
|
|
|
- this.signum = (magnitude.size() == 0 ? 0 : signum);
|
|
|
- this.mag = magnitude;
|
|
|
- if (mag.size() >= MAX_MAG_LENGTH) {
|
|
|
- checkRange();
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- BigInteger(const std::vector<int>& magnitude, int signum) {
|
|
|
- signum = (magnitude.size() == 0 ? 0 : signum);
|
|
|
- mag = stripLeadingZeroBytes(magnitude);
|
|
|
- if (mag.size() >= MAX_MAG_LENGTH) {
|
|
|
- checkRange();
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- void checkRange() {
|
|
|
- if (mag.size() > MAX_MAG_LENGTH || mag.size() == MAX_MAG_LENGTH && mag[0] < 0) {
|
|
|
- reportOverflow();
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
- static void reportOverflow() {
|
|
|
- std::cout << "BigInteger would overflow supported range" << std::endl;
|
|
|
- throw 1;
|
|
|
- }
|
|
|
-
|
|
|
- //Static Factory Methods
|
|
|
-
|
|
|
-
|
|
|
- static BigInteger valueOf(std::int64_t val) {
|
|
|
- // If -MAX_CONSTANT < val < MAX_CONSTANT, return stashed constant
|
|
|
- if (val == 0)
|
|
|
- return ZERO;
|
|
|
- if (val > 0 && val <= MAX_CONSTANT)
|
|
|
- return posConst[(int) val];
|
|
|
- else if (val < 0 && val >= -MAX_CONSTANT)
|
|
|
- return negConst[(int) -val];
|
|
|
-
|
|
|
- return BigInteger(val);
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- BigInteger(std::int64_t val) {
|
|
|
- if (val < 0) {
|
|
|
- val = -val;
|
|
|
- signum = -1;
|
|
|
- } else {
|
|
|
- signum = 1;
|
|
|
- }
|
|
|
-
|
|
|
- int highWord = (int)(((uint64_t)val) >> 32);
|
|
|
- if (highWord == 0) {
|
|
|
- mag = std::vector<int>(1);
|
|
|
- mag[0] = (int)val;
|
|
|
- } else {
|
|
|
- mag = std::vector<int>(2);
|
|
|
- mag[0] = highWord;
|
|
|
- mag[1] = (int)val;
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- static BigInteger valueOf(std::vector<int> val) {
|
|
|
- return (val[0] > 0 ? BigInteger(val, 1) : BigInteger(val));
|
|
|
- }
|
|
|
-
|
|
|
- // Constants
|
|
|
-
|
|
|
-
|
|
|
- const static int MAX_CONSTANT = 16;
|
|
|
- static BigInteger posConst[] = new BigInteger[MAX_CONSTANT+1];
|
|
|
- static BigInteger negConst[] = new BigInteger[MAX_CONSTANT+1];
|
|
|
-
|
|
|
-
|
|
|
- static volatile BigInteger[][] powerCache;
|
|
|
-
|
|
|
-
|
|
|
- static const double[] logCache;
|
|
|
-
|
|
|
-
|
|
|
- static const double LOG_TWO = std::log(2.0);
|
|
|
-
|
|
|
- static void initStuff(){
|
|
|
- for (int i = 1; i <= MAX_CONSTANT; i++) {
|
|
|
- std::vector<int> magnitude(1);
|
|
|
- magnitude[0] = i;
|
|
|
- posConst[i] = BigInteger(magnitude, 1);
|
|
|
- negConst[i] = BigInteger(magnitude, -1);
|
|
|
- }
|
|
|
-
|
|
|
- /*
|
|
|
- * Initialize the cache of radix^(2^x) values used for base conversion
|
|
|
- * with just the very first value. Additional values will be created
|
|
|
- * on demand.
|
|
|
- */
|
|
|
- powerCache = new BigInteger[Character.MAX_RADIX+1][];
|
|
|
- logCache = new double[Character.MAX_RADIX+1];
|
|
|
-
|
|
|
- for (int i=Character.MIN_RADIX; i <= Character.MAX_RADIX; i++) {
|
|
|
- powerCache[i] = new BigInteger[] { BigInteger.valueOf(i) };
|
|
|
- logCache[i] = Math.log(i);
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- static const BigInteger ZERO = new BigInteger(std::vector<int>(0), 0);
|
|
|
-
|
|
|
-
|
|
|
- static const BigInteger ONE = valueOf(1);
|
|
|
-
|
|
|
-
|
|
|
- static const BigInteger TWO = valueOf(2);
|
|
|
-
|
|
|
-
|
|
|
- static const BigInteger NEGATIVE_ONE = valueOf(-1);
|
|
|
-
|
|
|
-
|
|
|
- static const BigInteger TEN = valueOf(10);
|
|
|
-
|
|
|
- // Arithmetic Operations
|
|
|
-
|
|
|
-
|
|
|
- BigInteger add(BigInteger val) {
|
|
|
- if (val.signum == 0)
|
|
|
- return *this;
|
|
|
- if (signum == 0)
|
|
|
- return val;
|
|
|
- if (val.signum == signum)
|
|
|
- return BigInteger(add(mag, val.mag), signum);
|
|
|
-
|
|
|
- int cmp = compareMagnitude(val);
|
|
|
- if (cmp == 0)
|
|
|
- return ZERO;
|
|
|
- std::vector<int> resultMag = (cmp > 0 ? subtract(mag, val.mag)
|
|
|
- : subtract(val.mag, mag));
|
|
|
- resultMag = trustedStripLeadingZeroInts(resultMag);
|
|
|
-
|
|
|
- return BigInteger(resultMag, cmp == signum ? 1 : -1);
|
|
|
- }
|
|
|
- static int signum(std::int64_t a){
|
|
|
- if(a == 0)return 0;
|
|
|
- if(a > 0)return 1;
|
|
|
- return -1;
|
|
|
- }
|
|
|
-
|
|
|
- BigInteger add(std::int64_t val) {
|
|
|
- if (val == 0)
|
|
|
- return *this;
|
|
|
- if (signum == 0)
|
|
|
- return valueOf(val);
|
|
|
- if (signum(val) == signum)
|
|
|
- return BigInteger(add(mag, Math.abs(val)), signum);
|
|
|
- int cmp = compareMagnitude(val);
|
|
|
- if (cmp == 0)
|
|
|
- return ZERO;
|
|
|
- std::vector<int> resultMag = (cmp > 0 ? subtract(mag, Math.abs(val)) : subtract(Math.abs(val), mag));
|
|
|
- resultMag = trustedStripLeadingZeroInts(resultMag);
|
|
|
- return BigInteger(resultMag, cmp == signum ? 1 : -1);
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- static std::vector<int> add(std::vector<int> x, std::uint64_t val) {
|
|
|
- std::vector<int> y;
|
|
|
- long sum = 0;
|
|
|
- int xIndex = x.size();
|
|
|
- std::vector<int> result;
|
|
|
- int highWord = (int)(val >> 32);
|
|
|
- if (highWord == 0) {
|
|
|
- result = std::vector<int>(xIndex);
|
|
|
- sum = (x[--xIndex] & LONG_MASK) + val;
|
|
|
- result[xIndex] = (int)sum;
|
|
|
- } else {
|
|
|
- if (xIndex == 1) {
|
|
|
- result = std::vector<int>(2);
|
|
|
- sum = val + (x[0] & LONG_MASK);
|
|
|
- result[1] = (int)sum;
|
|
|
- result[0] = (int)(sum >>> 32);
|
|
|
- return result;
|
|
|
- } else {
|
|
|
- result = std::vector<int>(xIndex);
|
|
|
- sum = (x[--xIndex] & LONG_MASK) + (val & LONG_MASK);
|
|
|
- result[xIndex] = (int)sum;
|
|
|
- sum = (x[--xIndex] & LONG_MASK) + (highWord & LONG_MASK) + (sum >>> 32);
|
|
|
- result[xIndex] = (int)sum;
|
|
|
- }
|
|
|
- }
|
|
|
- // Copy remainder of longer number while carry propagation is required
|
|
|
- bool carry = (sum >>> 32 != 0);
|
|
|
- while (xIndex > 0 && carry)
|
|
|
- carry = ((result[--xIndex] = x[xIndex] + 1) == 0);
|
|
|
- // Copy remainder of longer number
|
|
|
- while (xIndex > 0)
|
|
|
- result[--xIndex] = x[xIndex];
|
|
|
- // Grow result if necessary
|
|
|
- if (carry) {
|
|
|
- std::vector<int> bigger(result.size() + 1);
|
|
|
- //System.arraycopy(result, 0, bigger, 1, result.length);
|
|
|
- std::copy(result.begin(),result.end(), bigger.begin());
|
|
|
- bigger[0] = 0x01;
|
|
|
- return bigger;
|
|
|
- }
|
|
|
- return result;
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- static std::vector<int> add(std::vector<int> x, std::vector<int> y) {
|
|
|
- // If x is shorter, swap the two arrays
|
|
|
- if (x.size() < y.size()) {
|
|
|
- std::swap(x,y);
|
|
|
- /*int[] tmp = x;
|
|
|
- x = y;
|
|
|
- y = tmp;*/
|
|
|
- }
|
|
|
-
|
|
|
- int xIndex = x.size();
|
|
|
- int yIndex = y.size();
|
|
|
- std::vector<int> result(xIndex);// = new int[xIndex];
|
|
|
- long sum = 0;
|
|
|
- if (yIndex == 1) {
|
|
|
- sum = (x[--xIndex] & LONG_MASK) + (y[0] & LONG_MASK) ;
|
|
|
- result[xIndex] = (int)sum;
|
|
|
- } else {
|
|
|
- // Add common parts of both numbers
|
|
|
- while (yIndex > 0) {
|
|
|
- sum = (x[--xIndex] & LONG_MASK) +
|
|
|
- (y[--yIndex] & LONG_MASK) + (sum >>> 32);
|
|
|
- result[xIndex] = (int)sum;
|
|
|
- }
|
|
|
- }
|
|
|
- // Copy remainder of longer number while carry propagation is required
|
|
|
- bool carry = (sum >>> 32 != 0);
|
|
|
- while (xIndex > 0 && carry)
|
|
|
- carry = ((result[--xIndex] = x[xIndex] + 1) == 0);
|
|
|
-
|
|
|
- // Copy remainder of longer number
|
|
|
- while (xIndex > 0)
|
|
|
- result[--xIndex] = x[xIndex];
|
|
|
-
|
|
|
- // Grow result if necessary
|
|
|
- if (carry) {
|
|
|
- std::vector<int> bigger(result.size() + 1);
|
|
|
- std::copy(result.begin(),result.end(),bigger.begin());
|
|
|
- bigger[0] = 0x01;
|
|
|
- return bigger;
|
|
|
- }
|
|
|
- return result;
|
|
|
- }
|
|
|
-
|
|
|
- static std::vector<int> subtract(std::int64_t val, std::vector<int> little) {
|
|
|
- int highWord = (int)(((uint64_t)val) >> 32);
|
|
|
- if (highWord == 0) {
|
|
|
- std::vector<int> result(1);
|
|
|
- result[0] = (int)(val - (little[0] & LONG_MASK));
|
|
|
- return result;
|
|
|
- } else {
|
|
|
- std::vector<int> result(2);
|
|
|
- if (little.size() == 1) {
|
|
|
- std::int64_t difference = ((int)val & LONG_MASK) - (little[0] & LONG_MASK);
|
|
|
- result[1] = (int)difference;
|
|
|
- // Subtract remainder of longer number while borrow propagates
|
|
|
- bool borrow = (difference >> 32 != 0);
|
|
|
- if (borrow) {
|
|
|
- result[0] = highWord - 1;
|
|
|
- } else { // Copy remainder of longer number
|
|
|
- result[0] = highWord;
|
|
|
- }
|
|
|
- return result;
|
|
|
- } else { // little.length == 2
|
|
|
- std::int64_t difference = ((int)val & LONG_MASK) - (little[1] & LONG_MASK);
|
|
|
- result[1] = (int)difference;
|
|
|
- difference = (highWord & LONG_MASK) - (little[0] & LONG_MASK) + (difference >> 32);
|
|
|
- result[0] = (int)difference;
|
|
|
- return result;
|
|
|
- }
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- static std::vector<int> subtract(std::vector<int> big, std::int64_t val) {
|
|
|
- int highWord = (int)(((uint64_t)val) >> 32);
|
|
|
- int bigIndex = big.size();
|
|
|
- std::vector<int> result(bigIndex);
|
|
|
- long difference = 0;
|
|
|
-
|
|
|
- if (highWord == 0) {
|
|
|
- difference = (big[--bigIndex] & LONG_MASK) - val;
|
|
|
- result[bigIndex] = (int)difference;
|
|
|
- } else {
|
|
|
- difference = (big[--bigIndex] & LONG_MASK) - (val & LONG_MASK);
|
|
|
- result[bigIndex] = (int)difference;
|
|
|
- difference = (big[--bigIndex] & LONG_MASK) - (highWord & LONG_MASK) + (difference >> 32);
|
|
|
- result[bigIndex] = (int)difference;
|
|
|
- }
|
|
|
-
|
|
|
- // Subtract remainder of longer number while borrow propagates
|
|
|
- bool borrow = (difference >> 32 != 0);
|
|
|
- while (bigIndex > 0 && borrow)
|
|
|
- borrow = ((result[--bigIndex] = big[bigIndex] - 1) == -1);
|
|
|
-
|
|
|
- // Copy remainder of longer number
|
|
|
- while (bigIndex > 0)
|
|
|
- result[--bigIndex] = big[bigIndex];
|
|
|
-
|
|
|
- return result;
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- BigInteger subtract(BigInteger val) {
|
|
|
- if (val.signum == 0)
|
|
|
- return *this;
|
|
|
- if (signum == 0)
|
|
|
- return val.negate();
|
|
|
- if (val.signum != signum)
|
|
|
- return BigInteger(add(mag, val.mag), signum);
|
|
|
-
|
|
|
- int cmp = compareMagnitude(val);
|
|
|
- if (cmp == 0)
|
|
|
- return ZERO;
|
|
|
- std::vector<int> resultMag = (cmp > 0 ? subtract(mag, val.mag)
|
|
|
- : subtract(val.mag, mag));
|
|
|
- resultMag = trustedStripLeadingZeroInts(resultMag);
|
|
|
- return BigInteger(resultMag, cmp == signum ? 1 : -1);
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- static std::vector<int> subtract(std::vector<int> big, std::vector<int> little) {
|
|
|
- int bigIndex = big.size();
|
|
|
- std::vector<int> result(bigIndex);
|
|
|
- int littleIndex = little.size();
|
|
|
- std::int64_t difference = 0;
|
|
|
-
|
|
|
- // Subtract common parts of both numbers
|
|
|
- while (littleIndex > 0) {
|
|
|
- difference = (big[--bigIndex] & LONG_MASK) -
|
|
|
- (little[--littleIndex] & LONG_MASK) +
|
|
|
- (difference >> 32);
|
|
|
- result[bigIndex] = (int)difference;
|
|
|
- }
|
|
|
-
|
|
|
- // Subtract remainder of longer number while borrow propagates
|
|
|
- bool borrow = (difference >> 32 != 0);
|
|
|
- while (bigIndex > 0 && borrow)
|
|
|
- borrow = ((result[--bigIndex] = big[bigIndex] - 1) == -1);
|
|
|
-
|
|
|
- // Copy remainder of longer number
|
|
|
- while (bigIndex > 0)
|
|
|
- result[--bigIndex] = big[bigIndex];
|
|
|
-
|
|
|
- return result;
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- BigInteger multiply(BigInteger val) {
|
|
|
- if (val.signum == 0 || signum == 0)
|
|
|
- return ZERO;
|
|
|
-
|
|
|
- int xlen = mag.size();
|
|
|
-
|
|
|
- if (val == this && xlen > MULTIPLY_SQUARE_THRESHOLD) {
|
|
|
- return square();
|
|
|
- }
|
|
|
-
|
|
|
- int ylen = val.mag.size();
|
|
|
-
|
|
|
- if ((xlen < KARATSUBA_THRESHOLD) || (ylen < KARATSUBA_THRESHOLD)) {
|
|
|
- int resultSign = signum == val.signum ? 1 : -1;
|
|
|
- if (val.mag.size() == 1) {
|
|
|
- return multiplyByInt(mag,val.mag[0], resultSign);
|
|
|
- }
|
|
|
- if (mag.size() == 1) {
|
|
|
- return multiplyByInt(val.mag,mag[0], resultSign);
|
|
|
- }
|
|
|
- std::vector<int> result = multiplyToLen(mag, xlen,
|
|
|
- val.mag, ylen, null);
|
|
|
- result = trustedStripLeadingZeroInts(result);
|
|
|
- return BigInteger(result, resultSign);
|
|
|
- } else {
|
|
|
- if ((xlen < TOOM_COOK_THRESHOLD) && (ylen < TOOM_COOK_THRESHOLD)) {
|
|
|
- return multiplyKaratsuba(this, val);
|
|
|
- } else {
|
|
|
- return multiplyToomCook3(this, val);
|
|
|
- }
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
- static BigInteger multiplyByInt(std::vector<int> x, int y, int sign) {
|
|
|
- if (__builtin_popcount(y) == 1) {
|
|
|
- return BigInteger(shiftLeft(x,__builtin_ctz(y)), sign);
|
|
|
- }
|
|
|
- int xlen = x.size();
|
|
|
- std::vector<int> rmag(xlen + 1);
|
|
|
- long carry = 0;
|
|
|
- long yl = y & LONG_MASK;
|
|
|
- int rstart = rmag.size() - 1;
|
|
|
- for (int i = xlen - 1; i >= 0; i--) {
|
|
|
- std::uint64_t product = (x[i] & LONG_MASK) * yl + carry;
|
|
|
- rmag[rstart--] = (int)product;
|
|
|
- carry = product >> 32;
|
|
|
- }
|
|
|
- if (carry == 0L) {
|
|
|
- rmag.erase(rmag.begin());
|
|
|
- } else {
|
|
|
- rmag[rstart] = (int)carry;
|
|
|
- }
|
|
|
- return BigInteger(rmag, sign);
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- BigInteger multiply(std::uint64_t v) {
|
|
|
- if (v == 0 || signum == 0)
|
|
|
- return ZERO;
|
|
|
- if (v == std::numeric_limits<std::uint64_t>::max())
|
|
|
- return multiply(BigInteger.valueOf(v));
|
|
|
- int rsign = (v > 0 ? signum : -signum);
|
|
|
- if (v < 0)
|
|
|
- v = -v;
|
|
|
- uint64_t dh = v >> 32; // higher order bits
|
|
|
- uint64_t dl = v & LONG_MASK; // lower order bits
|
|
|
-
|
|
|
- int xlen = mag.size();
|
|
|
- std::vector<int> value = mag;
|
|
|
- std::vector<int> rmag = (dh == 0L) ? (std::vector<int>([xlen + 1])) : (std::vector<int>([xlen + 2]));
|
|
|
- uint64_t carry = 0;
|
|
|
- int rstart = rmag.size() - 1;
|
|
|
- for (int i = xlen - 1; i >= 0; i--) {
|
|
|
- uint64_t product = (value[i] & LONG_MASK) * dl + carry;
|
|
|
- rmag[rstart--] = (int)product;
|
|
|
- carry = product >> 32;
|
|
|
- }
|
|
|
- rmag[rstart] = (int)carry;
|
|
|
- if (dh != 0L) {
|
|
|
- carry = 0;
|
|
|
- rstart = rmag.size() - 2;
|
|
|
- for (int i = xlen - 1; i >= 0; i--) {
|
|
|
- long product = (value[i] & LONG_MASK) * dh +
|
|
|
- (rmag[rstart] & LONG_MASK) + carry;
|
|
|
- rmag[rstart--] = (int)product;
|
|
|
- carry = product >> 32;
|
|
|
- }
|
|
|
- rmag[0] = (int)carry;
|
|
|
- }
|
|
|
- if (carry == 0LL)
|
|
|
- rmag.erase(rmag.begin());
|
|
|
- return BigInteger(rmag, rsign);
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- static std::vector<int> multiplyToLen(std::vector<int> x, int xlen, std::vector<int> y, int ylen, std::vector<int> z) {
|
|
|
- int xstart = xlen - 1;
|
|
|
- int ystart = ylen - 1;
|
|
|
-
|
|
|
- if (z.size() < (xlen+ ylen))
|
|
|
- z = std::vector<int>(xlen+ylen);
|
|
|
-
|
|
|
- uint64_t carry = 0;
|
|
|
- for (int j=ystart, k=ystart+1+xstart; j >= 0; j--, k--) {
|
|
|
- long product = (y[j] & LONG_MASK) *
|
|
|
- (x[xstart] & LONG_MASK) + carry;
|
|
|
- z[k] = (int)product;
|
|
|
- carry = product >> 32;
|
|
|
- }
|
|
|
- z[xstart] = (int)carry;
|
|
|
-
|
|
|
- for (int i = xstart-1; i >= 0; i--) {
|
|
|
- carry = 0;
|
|
|
- for (int j=ystart, k=ystart+1+i; j >= 0; j--, k--) {
|
|
|
- long product = (y[j] & LONG_MASK) *
|
|
|
- (x[i] & LONG_MASK) +
|
|
|
- (z[k] & LONG_MASK) + carry;
|
|
|
- z[k] = (int)product;
|
|
|
- carry = product >> 32;
|
|
|
- }
|
|
|
- z[i] = (int)carry;
|
|
|
- }
|
|
|
- return z;
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- static BigInteger multiplyKaratsuba(BigInteger x, BigInteger y) {
|
|
|
- int xlen = x.mag.size();
|
|
|
- int ylen = y.mag.size();
|
|
|
-
|
|
|
- // The number of ints in each half of the number.
|
|
|
- int half = (std::max(xlen, ylen)+1) / 2;
|
|
|
-
|
|
|
- // xl and yl are the lower halves of x and y respectively,
|
|
|
- // xh and yh are the upper halves.
|
|
|
- BigInteger xl = x.getLower(half);
|
|
|
- BigInteger xh = x.getUpper(half);
|
|
|
- BigInteger yl = y.getLower(half);
|
|
|
- BigInteger yh = y.getUpper(half);
|
|
|
-
|
|
|
- BigInteger p1 = xh.multiply(yh); // p1 = xh*yh
|
|
|
- BigInteger p2 = xl.multiply(yl); // p2 = xl*yl
|
|
|
-
|
|
|
- // p3=(xh+xl)*(yh+yl)
|
|
|
- BigInteger p3 = xh.add(xl).multiply(yh.add(yl));
|
|
|
-
|
|
|
- // result = p1 * 2^(32*2*half) + (p3 - p1 - p2) * 2^(32*half) + p2
|
|
|
- BigInteger result = p1.shiftLeft(32*half).add(p3.subtract(p1).subtract(p2)).shiftLeft(32*half).add(p2);
|
|
|
-
|
|
|
- if (x.signum != y.signum) {
|
|
|
- return result.negate();
|
|
|
- } else {
|
|
|
- return result;
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- static BigInteger multiplyToomCook3(BigInteger a, BigInteger b) {
|
|
|
- int alen = a.mag.size();
|
|
|
- int blen = b.mag.size();
|
|
|
-
|
|
|
- int largest = std::max(alen, blen);
|
|
|
-
|
|
|
- // k is the size (in ints) of the lower-order slices.
|
|
|
- int k = (largest+2)/3; // Equal to ceil(largest/3)
|
|
|
-
|
|
|
- // r is the size (in ints) of the highest-order slice.
|
|
|
- int r = largest - 2*k;
|
|
|
-
|
|
|
- // Obtain slices of the numbers. a2 and b2 are the most significant
|
|
|
- // bits of the numbers a and b, and a0 and b0 the least significant.
|
|
|
- BigInteger a0, a1, a2, b0, b1, b2;
|
|
|
- a2 = a.getToomSlice(k, r, 0, largest);
|
|
|
- a1 = a.getToomSlice(k, r, 1, largest);
|
|
|
- a0 = a.getToomSlice(k, r, 2, largest);
|
|
|
- b2 = b.getToomSlice(k, r, 0, largest);
|
|
|
- b1 = b.getToomSlice(k, r, 1, largest);
|
|
|
- b0 = b.getToomSlice(k, r, 2, largest);
|
|
|
-
|
|
|
- BigInteger v0, v1, v2, vm1, vinf, t1, t2, tm1, da1, db1;
|
|
|
-
|
|
|
- v0 = a0.multiply(b0);
|
|
|
- da1 = a2.add(a0);
|
|
|
- db1 = b2.add(b0);
|
|
|
- vm1 = da1.subtract(a1).multiply(db1.subtract(b1));
|
|
|
- da1 = da1.add(a1);
|
|
|
- db1 = db1.add(b1);
|
|
|
- v1 = da1.multiply(db1);
|
|
|
- v2 = da1.add(a2).shiftLeft(1).subtract(a0).multiply(
|
|
|
- db1.add(b2).shiftLeft(1).subtract(b0));
|
|
|
- vinf = a2.multiply(b2);
|
|
|
-
|
|
|
- // The algorithm requires two divisions by 2 and one by 3.
|
|
|
- // All divisions are known to be exact, that is, they do not produce
|
|
|
- // remainders, and all results are positive. The divisions by 2 are
|
|
|
- // implemented as right shifts which are relatively efficient, leaving
|
|
|
- // only an exact division by 3, which is done by a specialized
|
|
|
- // linear-time algorithm.
|
|
|
- t2 = v2.subtract(vm1).exactDivideBy3();
|
|
|
- tm1 = v1.subtract(vm1).shiftRight(1);
|
|
|
- t1 = v1.subtract(v0);
|
|
|
- t2 = t2.subtract(t1).shiftRight(1);
|
|
|
- t1 = t1.subtract(tm1).subtract(vinf);
|
|
|
- t2 = t2.subtract(vinf.shiftLeft(1));
|
|
|
- tm1 = tm1.subtract(t2);
|
|
|
-
|
|
|
- // Number of bits to shift left.
|
|
|
- int ss = k*32;
|
|
|
-
|
|
|
- BigInteger result = vinf.shiftLeft(ss).add(t2).shiftLeft(ss).add(t1).shiftLeft(ss).add(tm1).shiftLeft(ss).add(v0);
|
|
|
-
|
|
|
- if (a.signum != b.signum) {
|
|
|
- return result.negate();
|
|
|
- } else {
|
|
|
- return result;
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
-
|
|
|
- BigInteger getToomSlice(int lowerSize, int upperSize, int slice,
|
|
|
- int fullsize) {
|
|
|
- int start, end, sliceSize, len, offset;
|
|
|
-
|
|
|
- len = mag.size();
|
|
|
- offset = fullsize - len;
|
|
|
-
|
|
|
- if (slice == 0) {
|
|
|
- start = 0 - offset;
|
|
|
- end = upperSize - 1 - offset;
|
|
|
- } else {
|
|
|
- start = upperSize + (slice-1)*lowerSize - offset;
|
|
|
- end = start + lowerSize - 1;
|
|
|
- }
|
|
|
-
|
|
|
- if (start < 0) {
|
|
|
- start = 0;
|
|
|
- }
|
|
|
- if (end < 0) {
|
|
|
- return ZERO;
|
|
|
- }
|
|
|
-
|
|
|
- sliceSize = (end-start) + 1;
|
|
|
-
|
|
|
- if (sliceSize <= 0) {
|
|
|
- return ZERO;
|
|
|
- }
|
|
|
-
|
|
|
- // While performing Toom-Cook, all slices are positive and
|
|
|
- // the sign is adjusted when the const number is composed.
|
|
|
- if (start == 0 && sliceSize >= len) {
|
|
|
- return this->abs();
|
|
|
- }
|
|
|
-
|
|
|
- std::vector<int> intSlice(sliceSize);
|
|
|
- std::copy(mag.begin() + start,mag.begin() + start + sliceSize, intSlice.begin());
|
|
|
-
|
|
|
- return BigInteger(trustedStripLeadingZeroInts(intSlice), 1);
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- BigInteger exactDivideBy3() {
|
|
|
- int len = mag.size();
|
|
|
- std::vector<int> result(len);
|
|
|
- std::int64_t x, w, q, borrow;
|
|
|
- borrow = 0L;
|
|
|
- for (int i = len-1; i >= 0; i--) {
|
|
|
- x = (mag[i] & LONG_MASK);
|
|
|
- w = x - borrow;
|
|
|
- if (borrow > x) { // Did we make the number go negative?
|
|
|
- borrow = 1LL;
|
|
|
- } else {
|
|
|
- borrow = 0LL;
|
|
|
- }
|
|
|
-
|
|
|
- // 0xAAAAAAAB is the modular inverse of 3 (mod 2^32). Thus,
|
|
|
- // the effect of this is to divide by 3 (mod 2^32).
|
|
|
- // This is much faster than division on most architectures.
|
|
|
- q = (w * 0xAAAAAAABL) & LONG_MASK;
|
|
|
- result[i] = (int) q;
|
|
|
-
|
|
|
- // Now check the borrow. The second check can of course be
|
|
|
- // eliminated if the first fails.
|
|
|
- if (q >= 0x55555556L) {
|
|
|
- borrow++;
|
|
|
- if (q >= 0xAAAAAAABL)
|
|
|
- borrow++;
|
|
|
- }
|
|
|
- }
|
|
|
- result = trustedStripLeadingZeroInts(result);
|
|
|
- return BigInteger(result, signum);
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- BigInteger getLower(int n) {
|
|
|
- int len = mag.size();
|
|
|
-
|
|
|
- if (len <= n) {
|
|
|
- return abs();
|
|
|
- }
|
|
|
-
|
|
|
- std::vector<int>lowerInts(n);
|
|
|
- std::copy(mag.begin() + (len - n), mag.end(), lowerInts.begin());
|
|
|
-
|
|
|
- return BigInteger(trustedStripLeadingZeroInts(lowerInts), 1);
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- BigInteger getUpper(int n) {
|
|
|
- int len = mag.size();
|
|
|
-
|
|
|
- if (len <= n) {
|
|
|
- return ZERO;
|
|
|
- }
|
|
|
-
|
|
|
- int upperLen = len - n;
|
|
|
- std::vector<int> upperInts(upperLen);
|
|
|
- std::copy(mag.begin(), mag.begin() + upperLen, upperInts);
|
|
|
-
|
|
|
- return BigInteger(trustedStripLeadingZeroInts(upperInts), 1);
|
|
|
- }
|
|
|
-
|
|
|
- // Squaring
|
|
|
-
|
|
|
-
|
|
|
- BigInteger square() {
|
|
|
- if (signum == 0) {
|
|
|
- return ZERO;
|
|
|
- }
|
|
|
- int len = mag.size();
|
|
|
-
|
|
|
- if (len < KARATSUBA_SQUARE_THRESHOLD) {
|
|
|
- std::vector<int> z = squareToLen(mag, len, null);
|
|
|
- return BigInteger(trustedStripLeadingZeroInts(z), 1);
|
|
|
- } else {
|
|
|
- if (len < TOOM_COOK_SQUARE_THRESHOLD) {
|
|
|
- return squareKaratsuba();
|
|
|
- } else {
|
|
|
- return squareToomCook3();
|
|
|
- }
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- static std::vector<int> squareToLen(std::vector<int> x, int len, std::vector<int> z) {
|
|
|
- int zlen = len << 1;
|
|
|
- if (z == null || z.size() < zlen)
|
|
|
- z = std::vector<int>(zlen);
|
|
|
-
|
|
|
- // Execute checks before calling intrinsified method.
|
|
|
- implSquareToLenChecks(x, len, z, zlen);
|
|
|
- return implSquareToLen(x, len, z, zlen);
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- static void implSquareToLenChecks(std::vector<int> x, int len, std::vector<int> z, int zlen){
|
|
|
- if (len < 1) {
|
|
|
- throw std::invalid_argument("invalid input length: " + len);
|
|
|
- }
|
|
|
- if (len > x.size()) {
|
|
|
- throw std::invalid_argument("input length out of bound: " +
|
|
|
- len + " > " + x.size());
|
|
|
- }
|
|
|
- if (len * 2 > z.size()) {
|
|
|
- throw std::invalid_argument("input length out of bound: " +
|
|
|
- (len * 2) + " > " + z.size());
|
|
|
- }
|
|
|
- if (zlen < 1) {
|
|
|
- throw std::invalid_argument("invalid input length: " + zlen);
|
|
|
- }
|
|
|
- if (zlen > z.size()) {
|
|
|
- throw std::invalid_argument("input length out of bound: " +
|
|
|
- len + " > " + z.size());
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- static std::vector<int> implSquareToLen(std::vector<int> x, int len, std::vector<int> z, int zlen) {
|
|
|
- /*
|
|
|
- * The algorithm used here is adapted from Colin Plumb's C library.
|
|
|
- * Technique: Consider the partial products in the multiplication
|
|
|
- * of "abcde" by itself:
|
|
|
- *
|
|
|
- * a b c d e
|
|
|
- * * a b c d e
|
|
|
- * ==================
|
|
|
- * ae be ce de ee
|
|
|
- * ad bd cd dd de
|
|
|
- * ac bc cc cd ce
|
|
|
- * ab bb bc bd be
|
|
|
- * aa ab ac ad ae
|
|
|
- *
|
|
|
- * Note that everything above the main diagonal:
|
|
|
- * ae be ce de = (abcd) * e
|
|
|
- * ad bd cd = (abc) * d
|
|
|
- * ac bc = (ab) * c
|
|
|
- * ab = (a) * b
|
|
|
- *
|
|
|
- * is a copy of everything below the main diagonal:
|
|
|
- * de
|
|
|
- * cd ce
|
|
|
- * bc bd be
|
|
|
- * ab ac ad ae
|
|
|
- *
|
|
|
- * Thus, the sum is 2 * (off the diagonal) + diagonal.
|
|
|
- *
|
|
|
- * This is accumulated beginning with the diagonal (which
|
|
|
- * consist of the squares of the digits of the input), which is then
|
|
|
- * divided by two, the off-diagonal added, and multiplied by two
|
|
|
- * again. The low bit is simply a copy of the low bit of the
|
|
|
- * input, so it doesn't need special care.
|
|
|
- */
|
|
|
-
|
|
|
- // Store the squares, right shifted one bit (i.e., divided by 2)
|
|
|
- int lastProductLowWord = 0;
|
|
|
- for (int j=0, i=0; j < len; j++) {
|
|
|
- std::int64_t piece = (x[j] & LONG_MASK);
|
|
|
- uint64_t product = piece * piece;
|
|
|
- z[i++] = (lastProductLowWord << 31) | (int)(product >> 33);
|
|
|
- z[i++] = (int)(product >> 1);
|
|
|
- lastProductLowWord = (int)product;
|
|
|
- }
|
|
|
-
|
|
|
- // Add in off-diagonal sums
|
|
|
- for (int i=len, offset=1; i > 0; i--, offset+=2) {
|
|
|
- int t = x[i-1];
|
|
|
- t = mulAdd(z, x, offset, i-1, t);
|
|
|
- addOne(z, offset-1, i, t);
|
|
|
- }
|
|
|
-
|
|
|
- // Shift back up and set low bit
|
|
|
- primitiveLeftShift(z, zlen, 1);
|
|
|
- z[zlen-1] |= x[len-1] & 1;
|
|
|
-
|
|
|
- return z;
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- BigInteger squareKaratsuba() {
|
|
|
- int half = (mag.size()+1) / 2;
|
|
|
-
|
|
|
- BigInteger xl = getLower(half);
|
|
|
- BigInteger xh = getUpper(half);
|
|
|
-
|
|
|
- BigInteger xhs = xh.square(); // xhs = xh^2
|
|
|
- BigInteger xls = xl.square(); // xls = xl^2
|
|
|
-
|
|
|
- // xh^2 << 64 + (((xl+xh)^2 - (xh^2 + xl^2)) << 32) + xl^2
|
|
|
- return xhs.shiftLeft(half*32).add(xl.add(xh).square().subtract(xhs.add(xls))).shiftLeft(half*32).add(xls);
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- BigInteger squareToomCook3() {
|
|
|
- int len = mag.size();
|
|
|
-
|
|
|
- // k is the size (in ints) of the lower-order slices.
|
|
|
- int k = (len+2)/3; // Equal to ceil(largest/3)
|
|
|
-
|
|
|
- // r is the size (in ints) of the highest-order slice.
|
|
|
- int r = len - 2*k;
|
|
|
-
|
|
|
- // Obtain slices of the numbers. a2 is the most significant
|
|
|
- // bits of the number, and a0 the least significant.
|
|
|
- BigInteger a0, a1, a2;
|
|
|
- a2 = getToomSlice(k, r, 0, len);
|
|
|
- a1 = getToomSlice(k, r, 1, len);
|
|
|
- a0 = getToomSlice(k, r, 2, len);
|
|
|
- BigInteger v0, v1, v2, vm1, vinf, t1, t2, tm1, da1;
|
|
|
-
|
|
|
- v0 = a0.square();
|
|
|
- da1 = a2.add(a0);
|
|
|
- vm1 = da1.subtract(a1).square();
|
|
|
- da1 = da1.add(a1);
|
|
|
- v1 = da1.square();
|
|
|
- vinf = a2.square();
|
|
|
- v2 = da1.add(a2).shiftLeft(1).subtract(a0).square();
|
|
|
-
|
|
|
- // The algorithm requires two divisions by 2 and one by 3.
|
|
|
- // All divisions are known to be exact, that is, they do not produce
|
|
|
- // remainders, and all results are positive. The divisions by 2 are
|
|
|
- // implemented as right shifts which are relatively efficient, leaving
|
|
|
- // only a division by 3.
|
|
|
- // The division by 3 is done by an optimized algorithm for this case.
|
|
|
- t2 = v2.subtract(vm1).exactDivideBy3();
|
|
|
- tm1 = v1.subtract(vm1).shiftRight(1);
|
|
|
- t1 = v1.subtract(v0);
|
|
|
- t2 = t2.subtract(t1).shiftRight(1);
|
|
|
- t1 = t1.subtract(tm1).subtract(vinf);
|
|
|
- t2 = t2.subtract(vinf.shiftLeft(1));
|
|
|
- tm1 = tm1.subtract(t2);
|
|
|
-
|
|
|
- // Number of bits to shift left.
|
|
|
- int ss = k*32;
|
|
|
-
|
|
|
- return vinf.shiftLeft(ss).add(t2).shiftLeft(ss).add(t1).shiftLeft(ss).add(tm1).shiftLeft(ss).add(v0);
|
|
|
- }
|
|
|
-
|
|
|
- // Division
|
|
|
-
|
|
|
-
|
|
|
- BigInteger divide(BigInteger val) {
|
|
|
- if (val.mag.size() < BURNIKEL_ZIEGLER_THRESHOLD ||
|
|
|
- mag.size() - val.mag.size() < BURNIKEL_ZIEGLER_OFFSET) {
|
|
|
- return divideKnuth(val);
|
|
|
- } else {
|
|
|
- return divideBurnikelZiegler(val);
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- BigInteger divideKnuth(BigInteger val) {
|
|
|
- MutableBigInteger q = MutableBigInteger(),
|
|
|
- a = MutableBigInteger(this->mag),
|
|
|
- b = MutableBigInteger(val.mag);
|
|
|
-
|
|
|
- a.divideKnuth(b, q, false);
|
|
|
- return q.toBigInteger(this->signum * val.signum);
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- std::vector<BigInteger> divideAndRemainder(BigInteger val) {
|
|
|
- if (val.mag.size() < BURNIKEL_ZIEGLER_THRESHOLD ||
|
|
|
- mag.size() - val.mag.size() < BURNIKEL_ZIEGLER_OFFSET) {
|
|
|
- return divideAndRemainderKnuth(val);
|
|
|
- } else {
|
|
|
- return divideAndRemainderBurnikelZiegler(val);
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- std::vector<BigInteger> divideAndRemainderKnuth(BigInteger val) {
|
|
|
- std::vector<BigInteger> result();
|
|
|
- MutableBigInteger q = MutableBigInteger(),
|
|
|
- a = MutableBigInteger(this->mag),
|
|
|
- b = MutableBigInteger(val.mag);
|
|
|
- MutableBigInteger r = a.divideKnuth(b, q);
|
|
|
- result[0] = q.toBigInteger(this->signum == val.signum ? 1 : -1);
|
|
|
- result[1] = r.toBigInteger(this->signum);
|
|
|
- return result;
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- BigInteger remainder(BigInteger val) {
|
|
|
- if (val.mag.size() < BURNIKEL_ZIEGLER_THRESHOLD ||
|
|
|
- mag.size() - val.mag.size() < BURNIKEL_ZIEGLER_OFFSET) {
|
|
|
- return remainderKnuth(val);
|
|
|
- } else {
|
|
|
- return remainderBurnikelZiegler(val);
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- BigInteger remainderKnuth(BigInteger val) {
|
|
|
- MutableBigInteger q = MutableBigInteger(),
|
|
|
- a = MutableBigInteger(this->mag),
|
|
|
- b = MutableBigInteger(val.mag);
|
|
|
-
|
|
|
- return a.divideKnuth(b, q).toBigInteger(this->signum);
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- BigInteger divideBurnikelZiegler(BigInteger val) {
|
|
|
- return divideAndRemainderBurnikelZiegler(val)[0];
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- BigInteger remainderBurnikelZiegler(BigInteger val) {
|
|
|
- return divideAndRemainderBurnikelZiegler(val)[1];
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- std::vector<BigInteger> divideAndRemainderBurnikelZiegler(BigInteger val) {
|
|
|
- MutableBigInteger q = MutableBigInteger();
|
|
|
- MutableBigInteger r = MutableBigInteger(*this).divideAndRemainderBurnikelZiegler(MutableBigInteger(val), q);
|
|
|
- BigInteger qBigInt = q.isZero() ? ZERO : q.toBigInteger(signum*val.signum);
|
|
|
- BigInteger rBigInt = r.isZero() ? ZERO : r.toBigInteger(signum);
|
|
|
- return std::vector<BigInteger> = {qBigInt, rBigInt};
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- BigInteger pow(int exponent) {
|
|
|
- assert(exponent > 0);
|
|
|
- if (signum == 0) {
|
|
|
- return (exponent == 0 ? ONE : this);
|
|
|
- }
|
|
|
-
|
|
|
- BigInteger partToSquare = this->abs();
|
|
|
-
|
|
|
- // Factor out powers of two from the base, as the exponentiation of
|
|
|
- // these can be done by left shifts only.
|
|
|
- // The remaining part can then be exponentiated faster. The
|
|
|
- // powers of two will be multiplied back at the end.
|
|
|
- int powersOfTwo = partToSquare.getLowestSetBit();
|
|
|
- long bitsToShift = (std::int64_t)powersOfTwo * exponent;
|
|
|
- if (bitsToShift > std::numeric_limits<int>::max()) {
|
|
|
- reportOverflow();
|
|
|
- }
|
|
|
-
|
|
|
- int remainingBits;
|
|
|
-
|
|
|
- // Factor the powers of two out quickly by shifting right, if needed.
|
|
|
- if (powersOfTwo > 0) {
|
|
|
- partToSquare = partToSquare.shiftRight(powersOfTwo);
|
|
|
- remainingBits = partToSquare.bitLength();
|
|
|
- if (remainingBits == 1) { // Nothing left but +/- 1?
|
|
|
- if (signum < 0 && (exponent&1) == 1) {
|
|
|
- return NEGATIVE_ONE.shiftLeft(powersOfTwo*exponent);
|
|
|
- } else {
|
|
|
- return ONE.shiftLeft(powersOfTwo*exponent);
|
|
|
- }
|
|
|
- }
|
|
|
- } else {
|
|
|
- remainingBits = partToSquare.bitLength();
|
|
|
- if (remainingBits == 1) { // Nothing left but +/- 1?
|
|
|
- if (signum < 0 && (exponent&1) == 1) {
|
|
|
- return NEGATIVE_ONE;
|
|
|
- } else {
|
|
|
- return ONE;
|
|
|
- }
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
- // This is a quick way to approximate the size of the result,
|
|
|
- // similar to doing log2[n] * exponent. This will give an upper bound
|
|
|
- // of how big the result can be, and which algorithm to use.
|
|
|
- long scaleFactor = (long)remainingBits * exponent;
|
|
|
-
|
|
|
- // Use slightly different algorithms for small and large operands.
|
|
|
- // See if the result will safely fit into a long. (Largest 2^63-1)
|
|
|
- if (partToSquare.mag.size() == 1 && scaleFactor <= 62) {
|
|
|
- // Small number algorithm. Everything fits into a long.
|
|
|
- int newSign = (signum <0 && (exponent&1) == 1 ? -1 : 1);
|
|
|
- long result = 1;
|
|
|
- long baseToPow2 = partToSquare.mag[0] & LONG_MASK;
|
|
|
-
|
|
|
- int workingExponent = exponent;
|
|
|
-
|
|
|
- // Perform exponentiation using repeated squaring trick
|
|
|
- while (workingExponent != 0) {
|
|
|
- if ((workingExponent & 1) == 1) {
|
|
|
- result = result * baseToPow2;
|
|
|
- }
|
|
|
-
|
|
|
- if ((workingExponent >>>= 1) != 0) {
|
|
|
- baseToPow2 = baseToPow2 * baseToPow2;
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
- // Multiply back the powers of two (quickly, by shifting left)
|
|
|
- if (powersOfTwo > 0) {
|
|
|
- if (bitsToShift + scaleFactor <= 62) { // Fits in long?
|
|
|
- return valueOf((result << bitsToShift) * newSign);
|
|
|
- } else {
|
|
|
- return valueOf(result*newSign).shiftLeft((int) bitsToShift);
|
|
|
- }
|
|
|
- }
|
|
|
- else {
|
|
|
- return valueOf(result*newSign);
|
|
|
- }
|
|
|
- } else {
|
|
|
- // Large number algorithm. This is basically identical to
|
|
|
- // the algorithm above, but calls multiply() and square()
|
|
|
- // which may use more efficient algorithms for large numbers.
|
|
|
- BigInteger answer = ONE;
|
|
|
-
|
|
|
- int workingExponent = exponent;
|
|
|
- // Perform exponentiation using repeated squaring trick
|
|
|
- while (workingExponent != 0) {
|
|
|
- if ((workingExponent & 1) == 1) {
|
|
|
- answer = answer.multiply(partToSquare);
|
|
|
- }
|
|
|
-
|
|
|
- if ((workingExponent >>>= 1) != 0) {
|
|
|
- partToSquare = partToSquare.square();
|
|
|
- }
|
|
|
- }
|
|
|
- // Multiply back the (exponentiated) powers of two (quickly,
|
|
|
- // by shifting left)
|
|
|
- if (powersOfTwo > 0) {
|
|
|
- answer = answer.shiftLeft(powersOfTwo*exponent);
|
|
|
- }
|
|
|
-
|
|
|
- if (signum < 0 && (exponent&1) == 1) {
|
|
|
- return answer.negate();
|
|
|
- } else {
|
|
|
- return answer;
|
|
|
- }
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- BigInteger gcd(BigInteger val) {
|
|
|
- if (val.signum == 0)
|
|
|
- return this->abs();
|
|
|
- else if (this.signum == 0)
|
|
|
- return val.abs();
|
|
|
-
|
|
|
- MutableBigInteger a = MutableBigInteger(this);
|
|
|
- MutableBigInteger b = MutableBigInteger(val);
|
|
|
-
|
|
|
- MutableBigInteger result = a.hybridGCD(b);
|
|
|
-
|
|
|
- return result.toBigInteger(1);
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- static int bitLengthForInt(int n) {
|
|
|
- return 32 - Integer.numberOfLeadingZeros(n);
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- static std::vector<int> leftShift(std::vector<int> a, int len, unsigned int n) {
|
|
|
- int nInts = n >> 5;
|
|
|
- int nBits = n&0x1F;
|
|
|
- int bitsInHighWord = bitLengthForInt(a[0]);
|
|
|
-
|
|
|
- // If shift can be done without recopy, do so
|
|
|
- if (n <= (32 - bitsInHighWord)) {
|
|
|
- primitiveLeftShift(a, len, nBits);
|
|
|
- return a;
|
|
|
- } else { // Array must be resized
|
|
|
- if (nBits <= (32 - bitsInHighWord)) {
|
|
|
- std::vector<int> result(nInts + len);
|
|
|
- std::copy(a.begin(),a.begin() + len, result.begin());
|
|
|
- primitiveLeftShift(result, result.size(), nBits);
|
|
|
- return result;
|
|
|
- } else {
|
|
|
- std::vector<int> resul(nInts + len + 1);
|
|
|
- std::copy(a.begin(), a.begin() + len, result.begin());
|
|
|
- primitiveRightShift(result, result.size(), 32 - nBits);
|
|
|
- return result;
|
|
|
- }
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
- // shifts a up to len right n bits assumes no leading zeros, 0<n<32
|
|
|
- static void primitiveRightShift(std::vector<int> a, int len, int n) {
|
|
|
- int n2 = 32 - n;
|
|
|
- for (int i=len-1, c=a[i]; i > 0; i--) {
|
|
|
- unsigned int b = c;
|
|
|
- c = a[i-1];
|
|
|
- a[i] = (c << n2) | (b >> n);
|
|
|
- }
|
|
|
- a[0] = ((unsigned int)a[0]) >> n;
|
|
|
- }
|
|
|
-
|
|
|
- // shifts a up to len left n bits assumes no leading zeros, 0<=n<32
|
|
|
- static void primitiveLeftShift(std::vector<int> a, int len, int n) {
|
|
|
- if (len == 0 || n == 0)
|
|
|
- return;
|
|
|
-
|
|
|
- int n2 = 32 - n;
|
|
|
- for (unsigned int i=0, c=a[i], m=i+len-1; i < m; i++) {
|
|
|
- int b = c;
|
|
|
- c = a[i+1];
|
|
|
- a[i] = (b << n) | (c >> n2);
|
|
|
- }
|
|
|
- a[len-1] <<= n;
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- static int bitLength(std::vector<int> val, int len) {
|
|
|
- if (len == 0)
|
|
|
- return 0;
|
|
|
- return ((len - 1) << 5) + bitLengthForInt(val[0]);
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- BigInteger abs() {
|
|
|
- return (signum >= 0 ? this : this->negate());
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- BigInteger negate() {
|
|
|
- return BigInteger(this->mag, -this->signum);
|
|
|
- }
|
|
|
-
|
|
|
- // Modular Arithmetic Operations
|
|
|
-
|
|
|
-
|
|
|
- BigInteger mod(BigInteger m) {
|
|
|
- assert(m.signum > 0)
|
|
|
-
|
|
|
- BigInteger result = this->remainder(m);
|
|
|
- return (result.signum >= 0 ? result : result.add(m));
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- BigInteger modPow(BigInteger exponent, BigInteger m) {
|
|
|
- assert(m.signum > 0)
|
|
|
-
|
|
|
- // Trivial cases
|
|
|
- if (exponent.signum == 0)
|
|
|
- return (m.equals(ONE) ? ZERO : ONE);
|
|
|
-
|
|
|
- if (this.equals(ONE))
|
|
|
- return (m.equals(ONE) ? ZERO : ONE);
|
|
|
-
|
|
|
- if (this.equals(ZERO) && exponent.signum >= 0)
|
|
|
- return ZERO;
|
|
|
-
|
|
|
- if (this.equals(negConst[1]) && (!exponent.testBit(0)))
|
|
|
- return (m.equals(ONE) ? ZERO : ONE);
|
|
|
-
|
|
|
- bool invertResult;
|
|
|
- if ((invertResult = (exponent.signum < 0)))
|
|
|
- exponent = exponent.negate();
|
|
|
-
|
|
|
- BigInteger base = (this->signum < 0 || this->compareTo(m) >= 0
|
|
|
- ? this->mod(m) : this);
|
|
|
- BigInteger result;
|
|
|
- if (m.testBit(0)) { // odd modulus
|
|
|
- result = base.oddModPow(exponent, m);
|
|
|
- } else {
|
|
|
- /*
|
|
|
- * Even modulus. Tear it into an "odd part" (m1) and power of two
|
|
|
- * (m2), exponentiate mod m1, manually exponentiate mod m2, and
|
|
|
- * use Chinese Remainder Theorem to combine results.
|
|
|
- */
|
|
|
-
|
|
|
- // Tear m apart into odd part (m1) and power of 2 (m2)
|
|
|
- int p = m.getLowestSetBit(); // Max pow of 2 that divides m
|
|
|
-
|
|
|
- BigInteger m1 = m.shiftRight(p); // m/2**p
|
|
|
- BigInteger m2 = ONE.shiftLeft(p); // 2**p
|
|
|
-
|
|
|
- // Calculate new base from m1
|
|
|
- BigInteger base2 = (this->signum < 0 || this->compareTo(m1) >= 0
|
|
|
- ? this->mod(m1) : this);
|
|
|
-
|
|
|
- // Caculate (base ** exponent) mod m1.
|
|
|
- BigInteger a1 = (m1.equals(ONE) ? ZERO :
|
|
|
- base2.oddModPow(exponent, m1));
|
|
|
-
|
|
|
- // Calculate (this ** exponent) mod m2
|
|
|
- BigInteger a2 = base.modPow2(exponent, p);
|
|
|
-
|
|
|
- // Combine results using Chinese Remainder Theorem
|
|
|
- BigInteger y1 = m2.modInverse(m1);
|
|
|
- BigInteger y2 = m1.modInverse(m2);
|
|
|
-
|
|
|
- if (m.mag.size() < MAX_MAG_LENGTH / 2) {
|
|
|
- result = a1.multiply(m2).multiply(y1).add(a2.multiply(m1).multiply(y2)).mod(m);
|
|
|
- } else {
|
|
|
- MutableBigInteger t1;
|
|
|
- new MutableBigInteger(a1.multiply(m2)).multiply(MutableBigInteger(y1), t1);
|
|
|
- MutableBigInteger t2;
|
|
|
- new MutableBigInteger(a2.multiply(m1)).multiply(MutableBigInteger(y2), t2);
|
|
|
- t1.add(t2);
|
|
|
- MutableBigInteger q;
|
|
|
- result = t1.divide(MutableBigInteger(m), q).toBigInteger();
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
- return (invertResult ? result.modInverse(m) : result);
|
|
|
- }
|
|
|
-
|
|
|
- // Montgomery multiplication. These are wrappers for
|
|
|
- // implMontgomeryXX routines which are expected to be replaced by
|
|
|
- // virtual machine intrinsics. We don't use the intrinsics for
|
|
|
- // very large operands: MONTGOMERY_INTRINSIC_THRESHOLD should be
|
|
|
- // larger than any reasonable crypto key.
|
|
|
- static std::vector<int> montgomeryMultiply(std::vector<int> a, std::vector<int> b, std::vector<int> n, int len, long inv,
|
|
|
- std::vector<int> product) {
|
|
|
- implMontgomeryMultiplyChecks(a, b, n, len, product);
|
|
|
- if (len > MONTGOMERY_INTRINSIC_THRESHOLD) {
|
|
|
- // Very long argument: do not use an intrinsic
|
|
|
- product = multiplyToLen(a, len, b, len, product);
|
|
|
- return montReduce(product, n, len, (int)inv);
|
|
|
- } else {
|
|
|
- return implMontgomeryMultiply(a, b, n, len, inv, materialize(product, len));
|
|
|
- }
|
|
|
- }
|
|
|
- static std::vector<int> montgomerySquare(std::vector<int> a, std::vector<int> n, int len, std::int64_t inv,
|
|
|
- std::vector<int> product) {
|
|
|
- implMontgomeryMultiplyChecks(a, a, n, len, product);
|
|
|
- if (len > MONTGOMERY_INTRINSIC_THRESHOLD) {
|
|
|
- // Very long argument: do not use an intrinsic
|
|
|
- product = squareToLen(a, len, product);
|
|
|
- return montReduce(product, n, len, (int)inv);
|
|
|
- } else {
|
|
|
- return implMontgomerySquare(a, n, len, inv, materialize(product, len));
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
- // Range-check everything.
|
|
|
- static void implMontgomeryMultiplyChecks
|
|
|
- (std::vector<int> a, std::vector<int> b, std::vector<int> n, int len, std::vector<int> product) {
|
|
|
- if (len % 2 != 0) {
|
|
|
- throw std::invalid_argument("input array length must be even: " + std::to_string(len));
|
|
|
- }
|
|
|
-
|
|
|
- if (len < 1) {
|
|
|
- throw std::invalid_argument("invalid input length: " + std::to_string(len));
|
|
|
- }
|
|
|
-
|
|
|
- if (len > a.size() ||
|
|
|
- len > b.size() ||
|
|
|
- len > n.size() ||
|
|
|
- (len > product.size())) {
|
|
|
- throw std::invalid_argument("input array length out of bound: " + len);
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
- // Make sure that the int array z (which is expected to contain
|
|
|
- // the result of a Montgomery multiplication) is present and
|
|
|
- // sufficiently large.
|
|
|
- static std::vector<int> materialize(std::vector<int> z, int len) {
|
|
|
- if (z.size() < len)
|
|
|
- z = std::vector<int>(len);
|
|
|
- return z;
|
|
|
- }
|
|
|
-
|
|
|
- // These methods are intended to be be replaced by virtual machine
|
|
|
- // intrinsics.
|
|
|
- static std::vector<int> implMontgomeryMultiply(std::vector<int> a, std::vector<int> b, std::vector<int> n, int len,
|
|
|
- std::int64_t inv, std::vector<int> product) {
|
|
|
- product = multiplyToLen(a, len, b, len, product);
|
|
|
- return montReduce(product, n, len, (int)inv);
|
|
|
- }
|
|
|
- static std::vector<int> implMontgomerySquare(std::vector<int> a, std::vector<int> n, int len,
|
|
|
- std::int64_t inv, std::vector<int> product) {
|
|
|
- product = squareToLen(a, len, product);
|
|
|
- return montReduce(product, n, len, (int)inv);
|
|
|
- }
|
|
|
-
|
|
|
- static std::vector<int> bnExpModThreshTable = {7, 25, 81, 241, 673, 1793,
|
|
|
- std::numeric_limits<int>::max()}; // Sentinel
|
|
|
-
|
|
|
-
|
|
|
- BigInteger oddModPow(BigInteger y, BigInteger z) {
|
|
|
- /*
|
|
|
- * The algorithm is adapted from Colin Plumb's C library.
|
|
|
- *
|
|
|
- * The window algorithm:
|
|
|
- * The idea is to keep a running product of b1 = n^(high-order bits of exp)
|
|
|
- * and then keep appending exponent bits to it. The following patterns
|
|
|
- * apply to a 3-bit window (k = 3):
|
|
|
- * To append 0: square
|
|
|
- * To append 1: square, multiply by n^1
|
|
|
- * To append 10: square, multiply by n^1, square
|
|
|
- * To append 11: square, square, multiply by n^3
|
|
|
- * To append 100: square, multiply by n^1, square, square
|
|
|
- * To append 101: square, square, square, multiply by n^5
|
|
|
- * To append 110: square, square, multiply by n^3, square
|
|
|
- * To append 111: square, square, square, multiply by n^7
|
|
|
- *
|
|
|
- * Since each pattern involves only one multiply, the longer the pattern
|
|
|
- * the better, except that a 0 (no multiplies) can be appended directly.
|
|
|
- * We precompute a table of odd powers of n, up to 2^k, and can then
|
|
|
- * multiply k bits of exponent at a time. Actually, assuming random
|
|
|
- * exponents, there is on average one zero bit between needs to
|
|
|
- * multiply (1/2 of the time there's none, 1/4 of the time there's 1,
|
|
|
- * 1/8 of the time, there's 2, 1/32 of the time, there's 3, etc.), so
|
|
|
- * you have to do one multiply per k+1 bits of exponent.
|
|
|
- *
|
|
|
- * The loop walks down the exponent, squaring the result buffer as
|
|
|
- * it goes. There is a wbits+1 bit lookahead buffer, buf, that is
|
|
|
- * filled with the upcoming exponent bits. (What is read after the
|
|
|
- * end of the exponent is unimportant, but it is filled with zero here.)
|
|
|
- * When the most-significant bit of this buffer becomes set, i.e.
|
|
|
- * (buf & tblmask) != 0, we have to decide what pattern to multiply
|
|
|
- * by, and when to do it. We decide, remember to do it in future
|
|
|
- * after a suitable number of squarings have passed (e.g. a pattern
|
|
|
- * of "100" in the buffer requires that we multiply by n^1 immediately;
|
|
|
- * a pattern of "110" calls for multiplying by n^3 after one more
|
|
|
- * squaring), clear the buffer, and continue.
|
|
|
- *
|
|
|
- * When we start, there is one more optimization: the result buffer
|
|
|
- * is implcitly one, so squaring it or multiplying by it can be
|
|
|
- * optimized away. Further, if we start with a pattern like "100"
|
|
|
- * in the lookahead window, rather than placing n into the buffer
|
|
|
- * and then starting to square it, we have already computed n^2
|
|
|
- * to compute the odd-powers table, so we can place that into
|
|
|
- * the buffer and save a squaring.
|
|
|
- *
|
|
|
- * This means that if you have a k-bit window, to compute n^z,
|
|
|
- * where z is the high k bits of the exponent, 1/2 of the time
|
|
|
- * it requires no squarings. 1/4 of the time, it requires 1
|
|
|
- * squaring, ... 1/2^(k-1) of the time, it reqires k-2 squarings.
|
|
|
- * And the remaining 1/2^(k-1) of the time, the top k bits are a
|
|
|
- * 1 followed by k-1 0 bits, so it again only requires k-2
|
|
|
- * squarings, not k-1. The average of these is 1. Add that
|
|
|
- * to the one squaring we have to do to compute the table,
|
|
|
- * and you'll see that a k-bit window saves k-2 squarings
|
|
|
- * as well as reducing the multiplies. (It actually doesn't
|
|
|
- * hurt in the case k = 1, either.)
|
|
|
- */
|
|
|
- // Special case for exponent of one
|
|
|
- if (y.equals(ONE))
|
|
|
- return *this;
|
|
|
-
|
|
|
- // Special case for base of zero
|
|
|
- if (signum == 0)
|
|
|
- return ZERO;
|
|
|
-
|
|
|
- std::vector<int> base = mag.clone();
|
|
|
- std::vector<int> exp = y.mag;
|
|
|
- std::vector<int> mod = z.mag;
|
|
|
- int modLen = mod.size();
|
|
|
-
|
|
|
- // Make modLen even. It is conventional to use a cryptographic
|
|
|
- // modulus that is 512, 768, 1024, or 2048 bits, so this code
|
|
|
- // will not normally be executed. However, it is necessary for
|
|
|
- // the correct functioning of the HotSpot intrinsics.
|
|
|
- if ((modLen & 1) != 0) {
|
|
|
- std::vector<int> x(modlen + 1);
|
|
|
- std::copy(mod.begin(), mod.begin() + modLen, x.begin());
|
|
|
- mod = std::move(x);
|
|
|
- modLen++;
|
|
|
- }
|
|
|
-
|
|
|
- // Select an appropriate window size
|
|
|
- int wbits = 0;
|
|
|
- int ebits = bitLength(exp, exp.size());
|
|
|
- // if exponent is 65537 (0x10001), use minimum window size
|
|
|
- if ((ebits != 17) || (exp[0] != 65537)) {
|
|
|
- while (ebits > bnExpModThreshTable[wbits]) {
|
|
|
- wbits++;
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
- // Calculate appropriate table size
|
|
|
- int tblmask = 1 << wbits;
|
|
|
-
|
|
|
- // Allocate table for precomputed odd powers of base in Montgomery form
|
|
|
- std::vector<std::vector<int>> table(tblmask);// = new int[tblmask][];
|
|
|
- for (int i=0; i < tblmask; i++)
|
|
|
- table[i] = std::vector<int>(modLen);
|
|
|
-
|
|
|
- // Compute the modular inverse of the least significant 64-bit
|
|
|
- // digit of the modulus
|
|
|
- std::int64_t n0 = (mod[modLen-1] & LONG_MASK) + ((mod[modLen-2] & LONG_MASK) << 32);
|
|
|
- std::int64_t inv = -MutableBigInteger.inverseMod64(n0);
|
|
|
-
|
|
|
- // Convert base to Montgomery form
|
|
|
- std::vector<int> a = leftShift(base, base.size(), modLen << 5);
|
|
|
-
|
|
|
- MutableBigInteger q = MutableBigInteger(),
|
|
|
- a2 = MutableBigInteger(a),
|
|
|
- b2 = MutableBigInteger(mod);
|
|
|
- b2.normalize(); // MutableBigInteger.divide() assumes that its
|
|
|
- // divisor is in normal form.
|
|
|
-
|
|
|
- MutableBigInteger r= a2.divide(b2, q);
|
|
|
- table[0] = r.toIntArray();
|
|
|
-
|
|
|
- // Pad table[0] with leading zeros so its length is at least modLen
|
|
|
- if (table[0].size() < modLen) {
|
|
|
- int offset = modLen - table[0].size();
|
|
|
- std::vector<int> t2(modLen);// = new int[modLen];
|
|
|
- std::copy(table[0].begin(),table[0].end(), t2.begin());
|
|
|
- table[0] = t2;
|
|
|
- }
|
|
|
-
|
|
|
- // Set b to the square of the base
|
|
|
- std::vector<int> b = montgomerySquare(table[0], mod, modLen, inv, null);
|
|
|
-
|
|
|
- // Set t to high half of b
|
|
|
- std::vector<int> t(b.begin(), b.begin() + modLen);
|
|
|
-
|
|
|
- // Fill in the table with odd powers of the base
|
|
|
- for (int i=1; i < tblmask; i++) {
|
|
|
- table[i] = montgomeryMultiply(t, table[i-1], mod, modLen, inv, null);
|
|
|
- }
|
|
|
-
|
|
|
- // Pre load the window that slides over the exponent
|
|
|
- unsigned int bitpos = 1 << ((ebits-1) & (32-1));
|
|
|
-
|
|
|
- unsigned int buf = 0;
|
|
|
- int elen = exp.size();
|
|
|
- int eIndex = 0;
|
|
|
- for (int i = 0; i <= wbits; i++) {
|
|
|
- buf = (buf << 1) | (((exp[eIndex] & bitpos) != 0) ? 1 : 0);
|
|
|
- bitpos >>= 1;
|
|
|
- if (bitpos == 0) {
|
|
|
- eIndex++;
|
|
|
- bitpos = 1 << (32-1);
|
|
|
- elen--;
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
- int multpos = ebits;
|
|
|
-
|
|
|
- // The first iteration, which is hoisted out of the main loop
|
|
|
- ebits--;
|
|
|
- bool isone = true;
|
|
|
-
|
|
|
- multpos = ebits - wbits;
|
|
|
- while ((buf & 1) == 0) {
|
|
|
- buf >>= 1;
|
|
|
- multpos++;
|
|
|
- }
|
|
|
-
|
|
|
- std::vector<int> mult = table[buf >>> 1];
|
|
|
-
|
|
|
- buf = 0;
|
|
|
- if (multpos == ebits)
|
|
|
- isone = false;
|
|
|
-
|
|
|
- // The main loop
|
|
|
- while (true) {
|
|
|
- ebits--;
|
|
|
- // Advance the window
|
|
|
- buf <<= 1;
|
|
|
-
|
|
|
- if (elen != 0) {
|
|
|
- buf |= ((exp[eIndex] & bitpos) != 0) ? 1 : 0;
|
|
|
- bitpos >>= 1;
|
|
|
- if (bitpos == 0) {
|
|
|
- eIndex++;
|
|
|
- bitpos = 1 << (32-1);
|
|
|
- elen--;
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
- // Examine the window for pending multiplies
|
|
|
- if ((buf & tblmask) != 0) {
|
|
|
- multpos = ebits - wbits;
|
|
|
- while ((buf & 1) == 0) {
|
|
|
- buf >>= 1;
|
|
|
- multpos++;
|
|
|
- }
|
|
|
- mult = table[buf >> 1];
|
|
|
- buf = 0;
|
|
|
- }
|
|
|
-
|
|
|
- // Perform multiply
|
|
|
- if (ebits == multpos) {
|
|
|
- if (isone) {
|
|
|
- b = mult;
|
|
|
- isone = false;
|
|
|
- } else {
|
|
|
- t = b;
|
|
|
- a = montgomeryMultiply(t, mult, mod, modLen, inv, a);
|
|
|
- t = a; a = b; b = t;
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
- // Check if done
|
|
|
- if (ebits == 0)
|
|
|
- break;
|
|
|
-
|
|
|
- // Square the input
|
|
|
- if (!isone) {
|
|
|
- t = b;
|
|
|
- a = montgomerySquare(t, mod, modLen, inv, a);
|
|
|
- t = a; a = b; b = t;
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
- // Convert result out of Montgomery form and return
|
|
|
- std::vector<int> t2(2 * modLen);
|
|
|
- std::copy(b.begin(), b.begin() + modLen, t2.begin() + modLen);
|
|
|
-
|
|
|
- b = montReduce(t2, mod, modLen, (int)inv);
|
|
|
-
|
|
|
- t2 = std::vector<int>(b.begin(),b.begin() + modLen);
|
|
|
-
|
|
|
- return BigInteger(1, t2);
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- static std::vector<int> montReduce(std::vector<int> n, std::vector<int> mod, int mlen, int inv) {
|
|
|
- int c = 0;
|
|
|
- int len = mlen;
|
|
|
- int offset = 0;
|
|
|
-
|
|
|
- do {
|
|
|
- int nEnd = n[n.size()-1-offset];
|
|
|
- int carry = mulAdd(n, mod, offset, mlen, inv * nEnd);
|
|
|
- c += addOne(n, offset, mlen, carry);
|
|
|
- offset++;
|
|
|
- } while (--len > 0);
|
|
|
-
|
|
|
- while (c > 0)
|
|
|
- c += subN(n, mod, mlen);
|
|
|
-
|
|
|
- while (intArrayCmpToLen(n, mod, mlen) >= 0)
|
|
|
- subN(n, mod, mlen);
|
|
|
-
|
|
|
- return n;
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- /*
|
|
|
- * Returns -1, 0 or +1 as big-endian unsigned int array arg1 is less than,
|
|
|
- * equal to, or greater than arg2 up to length len.
|
|
|
- */
|
|
|
- static int intArrayCmpToLen(std::vector<int> arg1, std::vector<int> arg2, int len) {
|
|
|
- for (int i=0; i < len; i++) {
|
|
|
- std::int64_t b1 = arg1[i] & LONG_MASK;
|
|
|
- std::int64_t b2 = arg2[i] & LONG_MASK;
|
|
|
- if (b1 < b2)
|
|
|
- return -1;
|
|
|
- if (b1 > b2)
|
|
|
- return 1;
|
|
|
- }
|
|
|
- return 0;
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- static int subN(std::vector<int> a, std::vector<int> b, int len) {
|
|
|
- std::int64_t sum = 0;
|
|
|
-
|
|
|
- while (--len >= 0) {
|
|
|
- sum = (a[len] & LONG_MASK) -
|
|
|
- (b[len] & LONG_MASK) + (sum >> 32);
|
|
|
- a[len] = (int)sum;
|
|
|
- }
|
|
|
-
|
|
|
- return (int)(sum >> 32);
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- static int mulAdd(std::vector<int> out, std::vector<int> in, int offset, int len, int k) {
|
|
|
- implMulAddCheck(out, in, offset, len, k);
|
|
|
- return implMulAdd(out, in, offset, len, k);
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- static void implMulAddCheck(std::vector<int> out, std::vector<int> in, int offset, int len, int k) {
|
|
|
- if (len > in.size()) {
|
|
|
- throw std::invalid_argument("input length is out of bound: " + std::to_string(len) + " > " + std::to_string(in.size()));
|
|
|
- }
|
|
|
- if (offset < 0) {
|
|
|
- throw std::invalid_argument("input offset is invalid: " + std::to_string(offset));
|
|
|
- }
|
|
|
- if (offset > (out.size() - 1)) {
|
|
|
- throw std::invalid_argument("input offset is out of bound: " + std::to_string(offset) + " > " + std::to_string(out.size() - 1));
|
|
|
- }
|
|
|
- if (len > (out.size() - offset)) {
|
|
|
- throw std::invalid_argument("input len is out of bound: " + std::to_string(len) + " > " + std::to_string(out.size() - offset));
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- static int implMulAdd(std::vector<int> out, std::vector<int> in, int offset, int len, int k) {
|
|
|
- std::int64_t kLong = k & LONG_MASK;
|
|
|
- std::int64_t carry = 0;
|
|
|
-
|
|
|
- offset = out.size()-offset - 1;
|
|
|
- for (int j=len-1; j >= 0; j--) {
|
|
|
- std::uint64_t product = (in[j] & LONG_MASK) * kLong +
|
|
|
- (out[offset] & LONG_MASK) + carry;
|
|
|
- out[offset--] = (int)product;
|
|
|
- carry = product >> 32;
|
|
|
- }
|
|
|
- return (int)carry;
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- static int addOne(std::vector<int> a, int offset, int mlen, int carry) {
|
|
|
- offset = a.size()-1-mlen-offset;
|
|
|
- std::uint64_t t = (a[offset] & LONG_MASK) + (carry & LONG_MASK);
|
|
|
-
|
|
|
- a[offset] = (int)t;
|
|
|
- if ((t >> 32) == 0)
|
|
|
- return 0;
|
|
|
- while (--mlen >= 0) {
|
|
|
- if (--offset < 0) { // Carry out of number
|
|
|
- return 1;
|
|
|
- } else {
|
|
|
- a[offset]++;
|
|
|
- if (a[offset] != 0)
|
|
|
- return 0;
|
|
|
- }
|
|
|
- }
|
|
|
- return 1;
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- BigInteger modPow2(BigInteger exponent, int p) {
|
|
|
- /*
|
|
|
- * Perform exponentiation using repeated squaring trick, chopping off
|
|
|
- * high order bits as indicated by modulus.
|
|
|
- */
|
|
|
- BigInteger result = ONE;
|
|
|
- BigInteger baseToPow2 = this->mod2(p);
|
|
|
- int expOffset = 0;
|
|
|
-
|
|
|
- int limit = exponent.bitLength();
|
|
|
-
|
|
|
- if (this.testBit(0))
|
|
|
- limit = (p-1) < limit ? (p-1) : limit;
|
|
|
-
|
|
|
- while (expOffset < limit) {
|
|
|
- if (exponent.testBit(expOffset))
|
|
|
- result = result.multiply(baseToPow2).mod2(p);
|
|
|
- expOffset++;
|
|
|
- if (expOffset < limit)
|
|
|
- baseToPow2 = baseToPow2.square().mod2(p);
|
|
|
- }
|
|
|
-
|
|
|
- return result;
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- BigInteger mod2(unsigned int p) {
|
|
|
- if (bitLength() <= p)
|
|
|
- return this;
|
|
|
-
|
|
|
- // Copy remaining ints of mag
|
|
|
- int numInts = (p + 31) >> 5;
|
|
|
- std::vector<int> mag(numInts,0);
|
|
|
- System.arraycopy(this->mag, (this->mag.size() - numInts), mag, 0, numInts);
|
|
|
-
|
|
|
- // Mask out any excess bits
|
|
|
- int excessBits = (numInts << 5) - p;
|
|
|
- mag[0] &= (1LL << (32-excessBits)) - 1;
|
|
|
-
|
|
|
- return (mag[0] == 0 ? BigInteger(1, mag) : BigInteger(mag, 1));
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- BigInteger modInverse(BigInteger m) {
|
|
|
- assert (m.signum == 1);
|
|
|
-
|
|
|
- if (m.equals(ONE))
|
|
|
- return ZERO;
|
|
|
-
|
|
|
- // Calculate (this mod m)
|
|
|
- BigInteger modVal = *this;
|
|
|
- if (signum < 0 || (this->compareMagnitude(m) >= 0))
|
|
|
- modVal = this->mod(m);
|
|
|
-
|
|
|
- if (modVal.equals(ONE))
|
|
|
- return ONE;
|
|
|
-
|
|
|
- MutableBigInteger a = MutableBigInteger(modVal);
|
|
|
- MutableBigInteger b = MutableBigInteger(m);
|
|
|
-
|
|
|
- MutableBigInteger result = a.mutableModInverse(b);
|
|
|
- return result.toBigInteger(1);
|
|
|
- }
|
|
|
-
|
|
|
- // Shift Operations
|
|
|
-
|
|
|
-
|
|
|
- BigInteger shiftLeft(int n) {
|
|
|
- if (signum == 0)
|
|
|
- return ZERO;
|
|
|
- if (n > 0) {
|
|
|
- return BigInteger(shiftLeft(mag, n), signum);
|
|
|
- } else if (n == 0) {
|
|
|
- return *this;
|
|
|
- } else {
|
|
|
- // Possible int overflow in (-n) is not a trouble,
|
|
|
- // because shiftRightImpl considers its argument unsigned
|
|
|
- return shiftRightImpl(-n);
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- static std::vector<int> shiftLeft(std::vector<int> mag, unsigned int n) {
|
|
|
- int nInts = n >> 5;
|
|
|
- int nBits = n & 0x1f;
|
|
|
- int magLen = mag.size();
|
|
|
- std::vector<int> newMag;
|
|
|
-
|
|
|
- if (nBits == 0) {
|
|
|
- newMag = std::vector<int>(magLen + nInts);
|
|
|
- System.arraycopy(mag.begin(), mag.end(), newMag.begin());
|
|
|
- } else {
|
|
|
- int i = 0;
|
|
|
- int nBits2 = 32 - nBits;
|
|
|
- int highBits = ((unsigned int)mag[0]) >> nBits2;
|
|
|
- if (highBits != 0) {
|
|
|
- newMag = std::vector<int>([magLen + nInts + 1]);
|
|
|
- newMag[i++] = highBits;
|
|
|
- } else {
|
|
|
- newMag = std::vector<int>(magLen + nInts);
|
|
|
- }
|
|
|
- int j=0;
|
|
|
- while (j < magLen-1)
|
|
|
- newMag[i++] = mag[j++] << nBits | ((unsigned int)mag[j]) >> nBits2;
|
|
|
- newMag[i] = mag[j] << nBits;
|
|
|
- }
|
|
|
- return newMag;
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- BigInteger shiftRight(int n) {
|
|
|
- if (signum == 0)
|
|
|
- return ZERO;
|
|
|
- if (n > 0) {
|
|
|
- return shiftRightImpl(n);
|
|
|
- } else if (n == 0) {
|
|
|
- return this;
|
|
|
- } else {
|
|
|
- return BigInteger(shiftLeft(mag, -n), signum);
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- BigInteger shiftRightImpl(unsigned int n) {
|
|
|
- int nInts = n >> 5;
|
|
|
- int nBits = n & 0x1f;
|
|
|
- int magLen = mag.size();
|
|
|
- std::vector<int> newMag;
|
|
|
-
|
|
|
- // Special case: entire contents shifted off the end
|
|
|
- if (nInts >= magLen)
|
|
|
- return (signum >= 0 ? ZERO : negConst[1]);
|
|
|
-
|
|
|
- if (nBits == 0) {
|
|
|
- int newMagLen = magLen - nInts;
|
|
|
-
|
|
|
- newMag = std::vector<int>(mag.begin(), mag.begin() + newMagLen);
|
|
|
- } else {
|
|
|
- int i = 0;
|
|
|
- int highBits = ((unsigned int)mag[0]) >> nBits;
|
|
|
- if (highBits != 0) {
|
|
|
- newMag = std::vector<int>(magLen - nInts);
|
|
|
- newMag[i++] = highBits;
|
|
|
- } else {
|
|
|
- newMag = std::vector<int>(magLen - nInts - 1);
|
|
|
- }
|
|
|
-
|
|
|
- int nBits2 = 32 - nBits;
|
|
|
- int j=0;
|
|
|
- while (j < magLen - nInts - 1)
|
|
|
- newMag[i++] = (mag[j++] << nBits2) | (((unsigned int)mag[j]) >> nBits);
|
|
|
- }
|
|
|
-
|
|
|
- if (signum < 0) {
|
|
|
- // Find out whether any one-bits were shifted off the end.
|
|
|
- bool onesLost = false;
|
|
|
- for (int i=magLen-1, j=magLen-nInts; i >= j && !onesLost; i--)
|
|
|
- onesLost = (mag[i] != 0);
|
|
|
- if (!onesLost && nBits != 0)
|
|
|
- onesLost = (mag[magLen - nInts - 1] << (32 - nBits) != 0);
|
|
|
-
|
|
|
- if (onesLost)
|
|
|
- newMag = javaIncrement(newMag);
|
|
|
- }
|
|
|
-
|
|
|
- return BigInteger(newMag, signum);
|
|
|
- }
|
|
|
-
|
|
|
- std::vector<int> javaIncrement(std::vector<int> val) {
|
|
|
- int lastSum = 0;
|
|
|
- for (int i=val.size()-1; i >= 0 && lastSum == 0; i--)
|
|
|
- lastSum = (val[i] += 1);
|
|
|
- if (lastSum == 0) {
|
|
|
- val = std::vector<int>(val.size()+1);
|
|
|
- val[0] = 1;
|
|
|
- }
|
|
|
- return val;
|
|
|
- }
|
|
|
-
|
|
|
- // Bitwise Operations
|
|
|
-
|
|
|
-
|
|
|
- BigInteger and(BigInteger val) {
|
|
|
- std::vector<int> result(std::max(intLength(), val.intLength()));
|
|
|
- for (int i=0; i < result.size(); i++)
|
|
|
- result[i] = (getInt(result.size()-i-1)
|
|
|
- & val.getInt(result.size()-i-1));
|
|
|
-
|
|
|
- return valueOf(result);
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- BigInteger or(BigInteger val) {
|
|
|
- std::vector<int> result(std::max(intLength(), val.intLength()));
|
|
|
- for (int i=0; i < result.size(); i++)
|
|
|
- result[i] = (getInt(result.size()-i-1)
|
|
|
- | val.getInt(result.size()-i-1));
|
|
|
-
|
|
|
- return valueOf(result);
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- BigInteger xor(BigInteger val) {
|
|
|
- std::vector<int> result(std::max(intLength(), val.intLength()));
|
|
|
- for (int i=0; i < result.size(); i++)
|
|
|
- result[i] = (getInt(result.size()-i-1)
|
|
|
- ^ val.getInt(result.size()-i-1));
|
|
|
-
|
|
|
- return valueOf(result);
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- BigInteger not() {
|
|
|
- std::vector<int> result(intLength());
|
|
|
- for (int i=0; i < result.size(); i++)
|
|
|
- result[i] = ~getInt(result.size()-i-1);
|
|
|
- return valueOf(result);
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- BigInteger andNot(BigInteger val) {
|
|
|
- std::vector<int> result(std::max(intLength(), val.intLength()));
|
|
|
- for (int i=0; i < result.size(); i++)
|
|
|
- result[i] = (getInt(result.size()-i-1)
|
|
|
- & ~val.getInt(result.size()-i-1));
|
|
|
-
|
|
|
- return valueOf(result);
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- // Single Bit Operations
|
|
|
-
|
|
|
-
|
|
|
- bool testBit(unsigned int n) {
|
|
|
- return (getInt(n >> 5) & (1 << (n & 31))) != 0;
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- BigInteger setBit(unsigned int n) {
|
|
|
-
|
|
|
- int intNum = n >> 5;
|
|
|
- std::vector<int> result(std::max(intLength(), intNum+2));
|
|
|
-
|
|
|
- for (int i=0; i < result.size(); i++)
|
|
|
- result[result.size()-i-1] = getInt(i);
|
|
|
-
|
|
|
- result[result.size()-intNum-1] |= (1 << (n & 31));
|
|
|
-
|
|
|
- return valueOf(result);
|
|
|
- }
|
|
|
- BigInteger clearBit(unsigned int n) {
|
|
|
-
|
|
|
- int intNum = n >> 5;
|
|
|
- std::vector<int> result(std::max(intLength(), ((n + 1) >> 5) + 1));
|
|
|
-
|
|
|
- for (int i=0; i < result.size(); i++)
|
|
|
- result[result.size()-i-1] = getInt(i);
|
|
|
-
|
|
|
- result[result.size()-intNum-1] &= ~(1 << (n & 31));
|
|
|
- return valueOf(result);
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- BigInteger flipBit(unsigned int n) {
|
|
|
-
|
|
|
- int intNum = n >> 5;
|
|
|
- std::vector<int> result(std::max(intLength(), intNum+2));
|
|
|
-
|
|
|
- for (int i=0; i < result.size(); i++)
|
|
|
- result[result.size()-i-1] = getInt(i);
|
|
|
-
|
|
|
- result[result.size()-intNum-1] ^= (1 << (n & 31));
|
|
|
-
|
|
|
- return valueOf(result);
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- int getLowestSetBit() {
|
|
|
- if (lsb == -2) { // lowestSetBit not initialized yet
|
|
|
- lsb = 0;
|
|
|
- if (signum == 0) {
|
|
|
- lsb -= 1;
|
|
|
- } else {
|
|
|
- // Search for lowest order nonzero int
|
|
|
- int i,b;
|
|
|
- for (i=0; (b = getInt(i)) == 0; i++)
|
|
|
- ;
|
|
|
- lsb += (i << 5) + __builtin_ctz(b);
|
|
|
- }
|
|
|
- lowestSetBit = lsb + 2;
|
|
|
- }
|
|
|
- return lsb;
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- // Miscellaneous Bit Operations
|
|
|
-
|
|
|
-
|
|
|
- int bitLength() {
|
|
|
- if (n == -1) { // bitLength not initialized yet
|
|
|
- int[] m = mag;
|
|
|
- int len = m.size();
|
|
|
- if (len == 0) {
|
|
|
- n = 0; // offset by one to initialize
|
|
|
- } else {
|
|
|
- // Calculate the bit length of the magnitude
|
|
|
- int magBitLength = ((len - 1) << 5) + bitLengthForInt(mag[0]);
|
|
|
- if (signum < 0) {
|
|
|
- // Check if magnitude is a power of two
|
|
|
- bool pow2 = (__builtin_popcount(mag[0]) == 1);
|
|
|
- for (int i=1; i< len && pow2; i++)
|
|
|
- pow2 = (mag[i] == 0);
|
|
|
-
|
|
|
- n = (pow2 ? magBitLength -1 : magBitLength);
|
|
|
- } else {
|
|
|
- n = magBitLength;
|
|
|
- }
|
|
|
- }
|
|
|
- bitLength = n + 1;
|
|
|
- }
|
|
|
- return n;
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- int bitCount() {
|
|
|
- if (bc == -1) { // bitCount not initialized yet
|
|
|
- bc = 0; // offset by one to initialize
|
|
|
- // Count the bits in the magnitude
|
|
|
- for (int i=0; i < mag.size(); i++)
|
|
|
- bc += __builtin_popcount(mag[i]);
|
|
|
- if (signum < 0) {
|
|
|
- // Count the trailing zeros in the magnitude
|
|
|
- int magTrailingZeroCount = 0, j;
|
|
|
- for (j=mag.size()-1; mag[j] == 0; j--)
|
|
|
- magTrailingZeroCount += 32;
|
|
|
- magTrailingZeroCount += __builtin_ctz(mag[j]);
|
|
|
- bc += magTrailingZeroCount - 1;
|
|
|
- }
|
|
|
- bitCount = bc + 1;
|
|
|
- }
|
|
|
- return bc;
|
|
|
- }
|
|
|
-
|
|
|
- // Primality Testing
|
|
|
-
|
|
|
-
|
|
|
- bool isProbablePrime(int certainty) {
|
|
|
- if (certainty <= 0)
|
|
|
- return true;
|
|
|
- BigInteger w = this->abs();
|
|
|
- if (w.equals(TWO))
|
|
|
- return true;
|
|
|
- if (!w.testBit(0) || w.equals(ONE))
|
|
|
- return false;
|
|
|
-
|
|
|
- return w.primeToCertainty(certainty, null);
|
|
|
- }
|
|
|
-
|
|
|
- // Comparison Operations
|
|
|
-
|
|
|
-
|
|
|
- int compareTo(BigInteger val) {
|
|
|
- if (signum == val.signum) {
|
|
|
- switch (signum) {
|
|
|
- case 1:
|
|
|
- return compareMagnitude(val);
|
|
|
- case -1:
|
|
|
- return val.compareMagnitude(this);
|
|
|
- default:
|
|
|
- return 0;
|
|
|
- }
|
|
|
- }
|
|
|
- return signum > val.signum ? 1 : -1;
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- const int compareMagnitude(BigInteger val) {
|
|
|
- std::vector<int> m1 = mag;
|
|
|
- int len1 = m1.size();
|
|
|
- std::vector<int> m2 = val.mag;
|
|
|
- int len2 = m2.size();
|
|
|
- if (len1 < len2)
|
|
|
- return -1;
|
|
|
- if (len1 > len2)
|
|
|
- return 1;
|
|
|
- for (int i = 0; i < len1; i++) {
|
|
|
- int a = m1[i];
|
|
|
- int b = m2[i];
|
|
|
- if (a != b)
|
|
|
- return ((a & LONG_MASK) < (b & LONG_MASK)) ? -1 : 1;
|
|
|
- }
|
|
|
- return 0;
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- const int compareMagnitude(std::uint64_t val) {
|
|
|
- std::vector<int> m1 = mag;
|
|
|
- int len = m1.size();
|
|
|
- if (len > 2) {
|
|
|
- return 1;
|
|
|
- }
|
|
|
- if (val < 0) {
|
|
|
- val = -val;
|
|
|
- }
|
|
|
- int highWord = (int)(val >> 32);
|
|
|
- if (highWord == 0) {
|
|
|
- if (len < 1)
|
|
|
- return -1;
|
|
|
- if (len > 1)
|
|
|
- return 1;
|
|
|
- int a = m1[0];
|
|
|
- int b = (int)val;
|
|
|
- if (a != b) {
|
|
|
- return ((a & LONG_MASK) < (b & LONG_MASK))? -1 : 1;
|
|
|
- }
|
|
|
- return 0;
|
|
|
- } else {
|
|
|
- if (len < 2)
|
|
|
- return -1;
|
|
|
- int a = m1[0];
|
|
|
- int b = highWord;
|
|
|
- if (a != b) {
|
|
|
- return ((a & LONG_MASK) < (b & LONG_MASK))? -1 : 1;
|
|
|
- }
|
|
|
- a = m1[1];
|
|
|
- b = (int)val;
|
|
|
- if (a != b) {
|
|
|
- return ((a & LONG_MASK) < (b & LONG_MASK))? -1 : 1;
|
|
|
- }
|
|
|
- return 0;
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- /*bool operator==(Object x) {
|
|
|
- // This test is just an optimization, which may or may not help
|
|
|
- if (x == this)
|
|
|
- return true;
|
|
|
-
|
|
|
- if (!(x instanceof BigInteger))
|
|
|
- return false;
|
|
|
-
|
|
|
- BigInteger xInt = (BigInteger) x;
|
|
|
- if (xInt.signum != signum)
|
|
|
- return false;
|
|
|
-
|
|
|
- int[] m = mag;
|
|
|
- int len = m.length;
|
|
|
- int[] xm = xInt.mag;
|
|
|
- if (len != xm.length)
|
|
|
- return false;
|
|
|
-
|
|
|
- for (int i = 0; i < len; i++)
|
|
|
- if (xm[i] != m[i])
|
|
|
- return false;
|
|
|
-
|
|
|
- return true;
|
|
|
- }*/
|
|
|
-
|
|
|
-
|
|
|
- BigInteger min(BigInteger val) {
|
|
|
- return (compareTo(val) < 0 ? *this : val);
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- BigInteger max(BigInteger val) {
|
|
|
- return (compareTo(val) > 0 ? *this : val);
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- // Hash Function
|
|
|
-
|
|
|
-
|
|
|
- int hashCode() {
|
|
|
- int hashCode = 0;
|
|
|
-
|
|
|
- for (int i=0; i < mag.size(); i++)
|
|
|
- hashCode = (int)(31*hashCode + (mag[i] & LONG_MASK));
|
|
|
-
|
|
|
- return hashCode * signum;
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- /*std::string toString(int radix) {
|
|
|
- if (signum == 0)
|
|
|
- return "0";
|
|
|
- if (radix < Character.MIN_RADIX || radix > Character.MAX_RADIX)
|
|
|
- radix = 10;
|
|
|
-
|
|
|
- // If it's small enough, use smallToString.
|
|
|
- if (mag.length <= SCHOENHAGE_BASE_CONVERSION_THRESHOLD)
|
|
|
- return smallToString(radix);
|
|
|
-
|
|
|
- // Otherwise use recursive toString, which requires positive arguments.
|
|
|
- // The results will be concatenated into this StringBuilder
|
|
|
- StringBuilder sb = new StringBuilder();
|
|
|
- if (signum < 0) {
|
|
|
- toString(this.negate(), sb, radix, 0);
|
|
|
- sb.insert(0, '-');
|
|
|
- }
|
|
|
- else
|
|
|
- toString(this, sb, radix, 0);
|
|
|
-
|
|
|
- return sb.toString();
|
|
|
- }*/
|
|
|
-
|
|
|
-
|
|
|
- /*String smallToString(int radix) {
|
|
|
- if (signum == 0) {
|
|
|
- return "0";
|
|
|
- }
|
|
|
-
|
|
|
- // Compute upper bound on number of digit groups and allocate space
|
|
|
- int maxNumDigitGroups = (4*mag.length + 6)/7;
|
|
|
- String digitGroup[] = new String[maxNumDigitGroups];
|
|
|
-
|
|
|
- // Translate number to string, a digit group at a time
|
|
|
- BigInteger tmp = this.abs();
|
|
|
- int numGroups = 0;
|
|
|
- while (tmp.signum != 0) {
|
|
|
- BigInteger d = longRadix[radix];
|
|
|
-
|
|
|
- MutableBigInteger q = MutableBigInteger(),
|
|
|
- a = MutableBigInteger(tmp.mag),
|
|
|
- b = MutableBigInteger(d.mag);
|
|
|
- MutableBigInteger r = a.divide(b, q);
|
|
|
- BigInteger q2 = q.toBigInteger(tmp.signum * d.signum);
|
|
|
- BigInteger r2 = r.toBigInteger(tmp.signum * d.signum);
|
|
|
-
|
|
|
- digitGroup[numGroups++] = Long.toString(r2.longValue(), radix);
|
|
|
- tmp = q2;
|
|
|
- }
|
|
|
-
|
|
|
- // Put sign (if any) and first digit group into result buffer
|
|
|
- StringBuilder buf = new StringBuilder(numGroups*digitsPerLong[radix]+1);
|
|
|
- if (signum < 0) {
|
|
|
- buf.append('-');
|
|
|
- }
|
|
|
- buf.append(digitGroup[numGroups-1]);
|
|
|
-
|
|
|
- // Append remaining digit groups padded with leading zeros
|
|
|
- for (int i=numGroups-2; i >= 0; i--) {
|
|
|
- // Prepend (any) leading zeros for this digit group
|
|
|
- int numLeadingZeros = digitsPerLong[radix]-digitGroup[i].length();
|
|
|
- if (numLeadingZeros != 0) {
|
|
|
- buf.append(zeros[numLeadingZeros]);
|
|
|
- }
|
|
|
- buf.append(digitGroup[i]);
|
|
|
- }
|
|
|
- return buf.toString();
|
|
|
- }*/
|
|
|
-
|
|
|
-
|
|
|
- /*static void toString(BigInteger u, StringBuilder sb, int radix, int digits) {
|
|
|
- if (u.mag.length <= SCHOENHAGE_BASE_CONVERSION_THRESHOLD) {
|
|
|
- String s = u.smallToString(radix);
|
|
|
-
|
|
|
- // Pad with internal zeros if necessary.
|
|
|
- // Don't pad if we're at the beginning of the string.
|
|
|
- if ((s.length() < digits) && (sb.length() > 0)) {
|
|
|
- for (int i=s.length(); i < digits; i++) { // May be a faster way to
|
|
|
- sb.append('0'); // do this?
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
- sb.append(s);
|
|
|
- return;
|
|
|
- }
|
|
|
-
|
|
|
- int b, n;
|
|
|
- b = u.bitLength();
|
|
|
-
|
|
|
- // Calculate a value for n in the equation radix^(2^n) = u
|
|
|
- // and subtract 1 from that value. This is used to find the
|
|
|
- // cache index that contains the best value to divide u.
|
|
|
- n = (int) Math.round(Math.log(b * LOG_TWO / logCache[radix]) / LOG_TWO - 1.0);
|
|
|
- BigInteger v = getRadixConversionCache(radix, n);
|
|
|
- BigInteger[] results;
|
|
|
- results = u.divideAndRemainder(v);
|
|
|
-
|
|
|
- int expectedDigits = 1 << n;
|
|
|
-
|
|
|
- // Now recursively build the two halves of each number.
|
|
|
- toString(results[0], sb, radix, digits-expectedDigits);
|
|
|
- toString(results[1], sb, radix, expectedDigits);
|
|
|
- }*/
|
|
|
-
|
|
|
-
|
|
|
- /*static BigInteger getRadixConversionCache(int radix, int exponent) {
|
|
|
- BigInteger[] cacheLine = powerCache[radix]; // volatile read
|
|
|
- if (exponent < cacheLine.length) {
|
|
|
- return cacheLine[exponent];
|
|
|
- }
|
|
|
-
|
|
|
- int oldLength = cacheLine.length;
|
|
|
- cacheLine = Arrays.copyOf(cacheLine, exponent + 1);
|
|
|
- for (int i = oldLength; i <= exponent; i++) {
|
|
|
- cacheLine[i] = cacheLine[i - 1].pow(2);
|
|
|
- }
|
|
|
-
|
|
|
- BigInteger[][] pc = powerCache; // volatile read again
|
|
|
- if (exponent >= pc[radix].length) {
|
|
|
- pc = pc.clone();
|
|
|
- pc[radix] = cacheLine;
|
|
|
- powerCache = pc; // volatile write, publish
|
|
|
- }
|
|
|
- return cacheLine[exponent];
|
|
|
- }*/
|
|
|
-
|
|
|
- /* zero[i] is a string of i consecutive zeros. */
|
|
|
- static std::string* zeros = new String[64];
|
|
|
- static {
|
|
|
- zeros[63] =
|
|
|
- "000000000000000000000000000000000000000000000000000000000000000";
|
|
|
- for (int i=0; i < 63; i++)
|
|
|
- zeros[i] = std::string(zeros[63].begin(),zeros[63].begin() + i));
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- String toString() {
|
|
|
- return toString(10);
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- std::vector<char> toByteArray() {
|
|
|
- int byteLen = bitLength()/8 + 1;
|
|
|
- std::vector<char> byteArray(byteLen);
|
|
|
-
|
|
|
- for (int i=byteLen-1, bytesCopied=4, nextInt=0, intIndex=0; i >= 0; i--) {
|
|
|
- if (bytesCopied == 4) {
|
|
|
- nextInt = getInt(intIndex++);
|
|
|
- bytesCopied = 1;
|
|
|
- } else {
|
|
|
- nextInt >> 8;
|
|
|
- bytesCopied++;
|
|
|
- }
|
|
|
- byteArray[i] = (byte)nextInt;
|
|
|
- }
|
|
|
- return byteArray;
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- int intValue() {
|
|
|
- int result = 0;
|
|
|
- result = getInt(0);
|
|
|
- return result;
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- std::int64_t longValue() {
|
|
|
- std::int64_t result = 0;
|
|
|
-
|
|
|
- for (int i=1; i >= 0; i--)
|
|
|
- result = (result << 32) + (getInt(i) & LONG_MASK);
|
|
|
- return result;
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- /*float floatValue() {
|
|
|
- if (signum == 0) {
|
|
|
- return 0.0f;
|
|
|
- }
|
|
|
-
|
|
|
- int exponent = ((mag.size() - 1) << 5) + bitLengthForInt(mag[0]) - 1;
|
|
|
-
|
|
|
- // exponent == floor(log2(abs(this)))
|
|
|
- if (exponent < 64 - 1) {
|
|
|
- return longValue();
|
|
|
- } else if (exponent > 127) {
|
|
|
- return signum > 0 ? (1.0f / 0.0f) : (-1.0f / 0.0f);
|
|
|
- }
|
|
|
-
|
|
|
- int shift = exponent - 24;
|
|
|
-
|
|
|
- int twiceSignifFloor;
|
|
|
- // twiceSignifFloor will be == abs().shiftRight(shift).intValue()
|
|
|
- // We do the shift into an int directly to improve performance.
|
|
|
-
|
|
|
- int nBits = shift & 0x1f;
|
|
|
- int nBits2 = 32 - nBits;
|
|
|
-
|
|
|
- if (nBits == 0) {
|
|
|
- twiceSignifFloor = mag[0];
|
|
|
- } else {
|
|
|
- twiceSignifFloor = ((unsigned)mag[0]) >> nBits;
|
|
|
- if (twiceSignifFloor == 0) {
|
|
|
- twiceSignifFloor = (mag[0] << nBits2) | (((unsigned)mag[1]) >> nBits);
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
- int signifFloor = twiceSignifFloor >> 1;
|
|
|
- signifFloor &= FloatConsts.SIGNIF_BIT_MASK; // remove the implied bit
|
|
|
- bool increment = (twiceSignifFloor & 1) != 0
|
|
|
- && ((signifFloor & 1) != 0 || abs().getLowestSetBit() < shift);
|
|
|
- int signifRounded = increment ? signifFloor + 1 : signifFloor;
|
|
|
- int bits = ((exponent + FloatConsts.EXP_BIAS))
|
|
|
- << (FloatConsts.SIGNIFICAND_WIDTH - 1);
|
|
|
- bits += signifRounded;
|
|
|
- bits |= signum & FloatConsts.SIGN_BIT_MASK;
|
|
|
- return Float.intBitsToFloat(bits);
|
|
|
- }*/
|
|
|
-
|
|
|
-
|
|
|
- /*double doubleValue() {
|
|
|
- if (signum == 0) {
|
|
|
- return 0.0;
|
|
|
- }
|
|
|
-
|
|
|
- int exponent = ((mag.length - 1) << 5) + bitLengthForInt(mag[0]) - 1;
|
|
|
-
|
|
|
- // exponent == floor(log2(abs(this))Double)
|
|
|
- if (exponent < Long.SIZE - 1) {
|
|
|
- return longValue();
|
|
|
- } else if (exponent > Double.MAX_EXPONENT) {
|
|
|
- return signum > 0 ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY;
|
|
|
- }
|
|
|
-
|
|
|
- int shift = exponent - DoubleConsts.SIGNIFICAND_WIDTH;
|
|
|
-
|
|
|
- long twiceSignifFloor;
|
|
|
- // twiceSignifFloor will be == abs().shiftRight(shift).longValue()
|
|
|
- // We do the shift into a long directly to improve performance.
|
|
|
-
|
|
|
- int nBits = shift & 0x1f;
|
|
|
- int nBits2 = 32 - nBits;
|
|
|
-
|
|
|
- int highBits;
|
|
|
- int lowBits;
|
|
|
- if (nBits == 0) {
|
|
|
- highBits = mag[0];
|
|
|
- lowBits = mag[1];
|
|
|
- } else {
|
|
|
- highBits = mag[0] >>> nBits;
|
|
|
- lowBits = (mag[0] << nBits2) | (mag[1] >>> nBits);
|
|
|
- if (highBits == 0) {
|
|
|
- highBits = lowBits;
|
|
|
- lowBits = (mag[1] << nBits2) | (mag[2] >>> nBits);
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
- twiceSignifFloor = ((highBits & LONG_MASK) << 32)
|
|
|
- | (lowBits & LONG_MASK);
|
|
|
-
|
|
|
- long signifFloor = twiceSignifFloor >> 1;
|
|
|
- signifFloor &= DoubleConsts.SIGNIF_BIT_MASK; // remove the implied bit
|
|
|
- bool increment = (twiceSignifFloor & 1) != 0
|
|
|
- && ((signifFloor & 1) != 0 || abs().getLowestSetBit() < shift);
|
|
|
- long signifRounded = increment ? signifFloor + 1 : signifFloor;
|
|
|
- long bits = (long) ((exponent + DoubleConsts.EXP_BIAS))
|
|
|
- << (DoubleConsts.SIGNIFICAND_WIDTH - 1);
|
|
|
- bits += signifRounded;
|
|
|
- bits |= signum & DoubleConsts.SIGN_BIT_MASK;
|
|
|
- return Double.longBitsToDouble(bits);
|
|
|
- }*/
|
|
|
-
|
|
|
-
|
|
|
- static std::vector<int> stripLeadingZeroInts(std::vector<int> val) {
|
|
|
- int vlen = val.size();
|
|
|
- int keep;
|
|
|
-
|
|
|
- // Find first nonzero byte
|
|
|
- for (keep = 0; keep < vlen && val[keep] == 0; keep++)
|
|
|
- ;
|
|
|
- return std::vector<int>(val.begin() + keep,val.begin() + keep + vlen);
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- static std::vector<int> trustedStripLeadingZeroInts(std::vector<int> val) {
|
|
|
- int vlen = val.size();
|
|
|
- int keep;
|
|
|
-
|
|
|
- // Find first nonzero byte
|
|
|
- for (keep = 0; keep < vlen && val[keep] == 0; keep++)
|
|
|
- ;
|
|
|
- return keep == 0 ? val : std::vector<int>(val.begin() + keep,val.begin() + keep + vlen);
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- static std::vector<int> stripLeadingZeroBytes(std::vector<char> a) {
|
|
|
- unsigned int byteLength = a.size();
|
|
|
- int keep;
|
|
|
-
|
|
|
- // Find first nonzero byte
|
|
|
- for (keep = 0; keep < byteLength && a[keep] == 0; keep++)
|
|
|
- ;
|
|
|
-
|
|
|
- // Allocate new array and copy relevant part of input array
|
|
|
- int intLength = ((byteLength - keep) + 3u) >> 2;
|
|
|
- std::vector<int> result(intLength);
|
|
|
- int b = byteLength - 1;
|
|
|
- for (int i = intLength-1; i >= 0; i--) {
|
|
|
- result[i] = a[b--] & 0xff;
|
|
|
- int bytesRemaining = b - keep + 1;
|
|
|
- int bytesToTransfer = std::min(3, bytesRemaining);
|
|
|
- for (int j=8; j <= (bytesToTransfer << 3); j += 8)
|
|
|
- result[i] |= ((a[b--] & 0xff) << j);
|
|
|
- }
|
|
|
- return result;
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- static std::vector<int> makePositive(std::vector<char> a) {
|
|
|
- int keep, k;
|
|
|
- unsigned int byteLength = a.size();
|
|
|
-
|
|
|
- // Find first non-sign (0xff) byte of input
|
|
|
- for (keep=0; keep < byteLength && a[keep] == -1; keep++);
|
|
|
-
|
|
|
-
|
|
|
- /* Allocate output array. If all non-sign bytes are 0x00, we must
|
|
|
- * allocate space for one extra output byte. */
|
|
|
- for (k=keep; k < byteLength && a[k] == 0; k++);
|
|
|
-
|
|
|
- int extraByte = (k == byteLength) ? 1 : 0;
|
|
|
- int intLength = ((byteLength - keep + extraByte) + 3u) >> 2;
|
|
|
- std::vector<int> result(intLength);
|
|
|
-
|
|
|
- /* Copy one's complement of input into output, leaving extra
|
|
|
- * byte (if it exists) == 0x00 */
|
|
|
- int b = byteLength - 1;
|
|
|
- for (int i = intLength-1; i >= 0; i--) {
|
|
|
- result[i] = a[b--] & 0xff;
|
|
|
- int numBytesToTransfer = Math.min(3, b-keep+1);
|
|
|
- if (numBytesToTransfer < 0)
|
|
|
- numBytesToTransfer = 0;
|
|
|
- for (int j=8; j <= 8*numBytesToTransfer; j += 8)
|
|
|
- result[i] |= ((a[b--] & 0xff) << j);
|
|
|
-
|
|
|
- // Mask indicates which bits must be complemented
|
|
|
- int mask = -1 >>> (8*(3-numBytesToTransfer));
|
|
|
- result[i] = ~result[i] & mask;
|
|
|
- }
|
|
|
-
|
|
|
- // Add one to one's complement to generate two's complement
|
|
|
- for (int i=result.size()-1; i >= 0; i--) {
|
|
|
- result[i] = (int)((result[i] & LONG_MASK) + 1);
|
|
|
- if (result[i] != 0)
|
|
|
- break;
|
|
|
- }
|
|
|
-
|
|
|
- return result;
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- static std::vector<int> makePositive(std::vector<int> a) {
|
|
|
- int keep, j;
|
|
|
-
|
|
|
- // Find first non-sign (0xffffffff) int of input
|
|
|
- for (keep=0; keep < a.size() && a[keep] == -1; keep++)
|
|
|
- ;
|
|
|
-
|
|
|
- /* Allocate output array. If all non-sign ints are 0x00, we must
|
|
|
- * allocate space for one extra output int. */
|
|
|
- for (j=keep; j < a.size() && a[j] == 0; j++)
|
|
|
- ;
|
|
|
- int extraInt = (j == a.size() ? 1 : 0);
|
|
|
- std::vector<int> result(a.size() - keep + extraInt);
|
|
|
-
|
|
|
- /* Copy one's complement of input into output, leaving extra
|
|
|
- * int (if it exists) == 0x00 */
|
|
|
- for (int i = keep; i < a.size(); i++)
|
|
|
- result[i - keep + extraInt] = ~a[i];
|
|
|
-
|
|
|
- // Add one to one's complement to generate two's complement
|
|
|
- for (int i = result.size()-1; ++result[i] == 0; i--);
|
|
|
-
|
|
|
- return result;
|
|
|
- }
|
|
|
-
|
|
|
- /*
|
|
|
- * The following two arrays are used for fast String conversions. Both
|
|
|
- * are indexed by radix. The first is the number of digits of the given
|
|
|
- * radix that can fit in a Java long without "going negative", i.e., the
|
|
|
- * highest integer n such that radix**n < 2**63. The second is the
|
|
|
- * "long radix" that tears each number into "long digits", each of which
|
|
|
- * consists of the number of digits in the corresponding element in
|
|
|
- * digitsPerLong (longRadix[i] = i**digitPerLong[i]). Both arrays have
|
|
|
- * nonsense values in their 0 and 1 elements, as radixes 0 and 1 are not
|
|
|
- * used.
|
|
|
- */
|
|
|
- static int digitsPerLong[] = {0, 0,
|
|
|
- 62, 39, 31, 27, 24, 22, 20, 19, 18, 18, 17, 17, 16, 16, 15, 15, 15, 14,
|
|
|
- 14, 14, 14, 13, 13, 13, 13, 13, 13, 12, 12, 12, 12, 12, 12, 12, 12};
|
|
|
-
|
|
|
- static BigInteger longRadix[] = {valueOf(0x4000000000000000L), valueOf(0x4000000000000000L),
|
|
|
- valueOf(0x4000000000000000L), valueOf(0x383d9170b85ff80bL),
|
|
|
- valueOf(0x4000000000000000L), valueOf(0x6765c793fa10079dL),
|
|
|
- valueOf(0x41c21cb8e1000000L), valueOf(0x3642798750226111L),
|
|
|
- valueOf(0x1000000000000000L), valueOf(0x12bf307ae81ffd59L),
|
|
|
- valueOf( 0xde0b6b3a7640000L), valueOf(0x4d28cb56c33fa539L),
|
|
|
- valueOf(0x1eca170c00000000L), valueOf(0x780c7372621bd74dL),
|
|
|
- valueOf(0x1e39a5057d810000L), valueOf(0x5b27ac993df97701L),
|
|
|
- valueOf(0x1000000000000000L), valueOf(0x27b95e997e21d9f1L),
|
|
|
- valueOf(0x5da0e1e53c5c8000L), valueOf( 0xb16a458ef403f19L),
|
|
|
- valueOf(0x16bcc41e90000000L), valueOf(0x2d04b7fdd9c0ef49L),
|
|
|
- valueOf(0x5658597bcaa24000L), valueOf( 0x6feb266931a75b7L),
|
|
|
- valueOf( 0xc29e98000000000L), valueOf(0x14adf4b7320334b9L),
|
|
|
- valueOf(0x226ed36478bfa000L), valueOf(0x383d9170b85ff80bL),
|
|
|
- valueOf(0x5a3c23e39c000000L), valueOf( 0x4e900abb53e6b71L),
|
|
|
- valueOf( 0x7600ec618141000L), valueOf( 0xaee5720ee830681L),
|
|
|
- valueOf(0x1000000000000000L), valueOf(0x172588ad4f5f0981L),
|
|
|
- valueOf(0x211e44f7d02c1000L), valueOf(0x2ee56725f06e5c71L),
|
|
|
- valueOf(0x41c21cb8e1000000L)};
|
|
|
-
|
|
|
- /*
|
|
|
- * These two arrays are the integer analogue of above.
|
|
|
- */
|
|
|
- static int digitsPerInt[] = {0, 0, 30, 19, 15, 13, 11,
|
|
|
- 11, 10, 9, 9, 8, 8, 8, 8, 7, 7, 7, 7, 7, 7, 7, 6, 6, 6, 6,
|
|
|
- 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 5};
|
|
|
-
|
|
|
- static int intRadix[] = {0, 0,
|
|
|
- 0x40000000, 0x4546b3db, 0x40000000, 0x48c27395, 0x159fd800,
|
|
|
- 0x75db9c97, 0x40000000, 0x17179149, 0x3b9aca00, 0xcc6db61,
|
|
|
- 0x19a10000, 0x309f1021, 0x57f6c100, 0xa2f1b6f, 0x10000000,
|
|
|
- 0x18754571, 0x247dbc80, 0x3547667b, 0x4c4b4000, 0x6b5a6e1d,
|
|
|
- 0x6c20a40, 0x8d2d931, 0xb640000, 0xe8d4a51, 0x1269ae40,
|
|
|
- 0x17179149, 0x1cb91000, 0x23744899, 0x2b73a840, 0x34e63b41,
|
|
|
- 0x40000000, 0x4cfa3cc1, 0x5c13d840, 0x6d91b519, 0x39aa400
|
|
|
- };
|
|
|
-
|
|
|
-
|
|
|
-
|
|
|
-
|
|
|
- int intLength() {
|
|
|
- return ((unsigned)bitLength() >> 5) + 1;
|
|
|
- }
|
|
|
-
|
|
|
- /* Returns sign bit */
|
|
|
- int signBit() {
|
|
|
- return signum < 0 ? 1 : 0;
|
|
|
- }
|
|
|
-
|
|
|
- /* Returns an int of sign bits */
|
|
|
- int signInt() {
|
|
|
- return signum < 0 ? -1 : 0;
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- int getInt(int n) {
|
|
|
- if (n < 0)
|
|
|
- return 0;
|
|
|
- if (n >= mag.size())
|
|
|
- return signInt();
|
|
|
-
|
|
|
- int magInt = mag[mag.size() - n - 1];
|
|
|
-
|
|
|
- return (signum >= 0 ? magInt :
|
|
|
- (n <= firstNonzeroIntNum() ? -magInt : ~magInt));
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- int firstNonzeroIntNum() {
|
|
|
- int fn = firstNonzeroIntNum - 2;
|
|
|
- if (fn == -2) { // firstNonzeroIntNum not initialized yet
|
|
|
- fn = 0;
|
|
|
-
|
|
|
- // Search for the first nonzero int
|
|
|
- int i;
|
|
|
- int mlen = mag.size();
|
|
|
- for (i = mlen - 1; i >= 0 && mag[i] == 0; i--);
|
|
|
- fn = mlen - i - 1;
|
|
|
- firstNonzeroIntNum = fn + 2; // offset by two to initialize
|
|
|
- }
|
|
|
- return fn;
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- static const long serialVersionUID = -8287574255936472291L;
|
|
|
-
|
|
|
-
|
|
|
-
|
|
|
-
|
|
|
- std::vector<char> magSerializedForm() {
|
|
|
- int len = mag.size();
|
|
|
-
|
|
|
- int bitLen = (len == 0 ? 0 : ((len - 1) << 5) + bitLengthForInt(mag[0]));
|
|
|
- int byteLen = (bitLen + 7) >> 3;
|
|
|
- std::vector<char> result(byteLen);
|
|
|
-
|
|
|
- for (int i = byteLen - 1, bytesCopied = 4, intIndex = len - 1, nextInt = 0;
|
|
|
- i >= 0; i--) {
|
|
|
- if (bytesCopied == 4) {
|
|
|
- nextInt = mag[intIndex--];
|
|
|
- bytesCopied = 1;
|
|
|
- } else {
|
|
|
- nextInt >>= 8;
|
|
|
- bytesCopied++;
|
|
|
- }
|
|
|
- result[i] = (byte)nextInt;
|
|
|
- }
|
|
|
- return result;
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- std::int64_t longValueExact() {
|
|
|
- if (mag.size() <= 2 && bitLength() <= 63)
|
|
|
- return longValue();
|
|
|
- else
|
|
|
- throw std::logic_error("BigInteger out of long range");
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- int intValueExact() {
|
|
|
- if (mag.size() <= 1 && bitLength() <= 31)
|
|
|
- return intValue();
|
|
|
- else
|
|
|
- throw std::logic_error("BigInteger out of int range");
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- short shortValueExact() {
|
|
|
- if (mag.size() <= 1 && bitLength() <= 31) {
|
|
|
- int value = intValue();
|
|
|
- if (value >= Short.MIN_VALUE && value <= std::numeric_limits<short>::max())
|
|
|
- return shortValue();
|
|
|
- }
|
|
|
- throw std::logic_error("BigInteger out of short range");
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- char byteValueExact() {
|
|
|
- if (mag.size() <= 1 && bitLength() <= 31) {
|
|
|
- int value = intValue();
|
|
|
- if (value >= Byte.MIN_VALUE && value <= std::numeric_limits<char>::max())
|
|
|
- return byteValue();
|
|
|
- }
|
|
|
- throw std::logic_error("BigInteger out of byte range");
|
|
|
- }
|
|
|
-};
|
|
|
-int main() {
|
|
|
- BigInteger a(6);
|
|
|
- a = a.pow(100);
|
|
|
-}
|
mawinkle