12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112111311141115111611171118111911201121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157115811591160116111621163116411651166116711681169117011711172117311741175117611771178117911801181118211831184118511861187118811891190119111921193119411951196119711981199120012011202120312041205120612071208120912101211121212131214121512161217121812191220122112221223122412251226122712281229123012311232123312341235123612371238123912401241124212431244124512461247124812491250125112521253125412551256125712581259126012611262126312641265126612671268126912701271127212731274127512761277127812791280128112821283128412851286128712881289129012911292129312941295129612971298129913001301130213031304130513061307130813091310131113121313131413151316131713181319132013211322132313241325132613271328132913301331133213331334133513361337133813391340134113421343134413451346134713481349135013511352135313541355135613571358135913601361136213631364136513661367136813691370137113721373137413751376137713781379138013811382138313841385138613871388138913901391139213931394139513961397139813991400140114021403140414051406140714081409141014111412141314141415141614171418141914201421142214231424142514261427142814291430143114321433143414351436143714381439144014411442144314441445144614471448144914501451145214531454145514561457145814591460146114621463146414651466146714681469147014711472147314741475147614771478147914801481148214831484148514861487148814891490149114921493149414951496149714981499150015011502150315041505150615071508150915101511151215131514151515161517151815191520152115221523152415251526152715281529153015311532153315341535153615371538153915401541154215431544154515461547154815491550155115521553155415551556155715581559156015611562156315641565156615671568156915701571157215731574157515761577157815791580158115821583158415851586158715881589159015911592159315941595159615971598159916001601160216031604160516061607160816091610161116121613161416151616161716181619162016211622162316241625162616271628162916301631163216331634163516361637163816391640164116421643164416451646164716481649165016511652165316541655165616571658165916601661166216631664166516661667166816691670167116721673167416751676167716781679168016811682168316841685168616871688168916901691169216931694169516961697169816991700170117021703170417051706170717081709171017111712171317141715171617171718171917201721172217231724172517261727172817291730173117321733173417351736173717381739174017411742174317441745174617471748174917501751175217531754175517561757175817591760176117621763176417651766176717681769177017711772177317741775177617771778177917801781178217831784178517861787178817891790179117921793179417951796179717981799180018011802180318041805180618071808180918101811181218131814181518161817181818191820182118221823182418251826182718281829183018311832183318341835183618371838183918401841184218431844184518461847184818491850185118521853185418551856185718581859186018611862186318641865186618671868186918701871187218731874187518761877187818791880188118821883188418851886188718881889189018911892189318941895189618971898189919001901190219031904190519061907190819091910191119121913191419151916191719181919192019211922192319241925192619271928192919301931193219331934193519361937193819391940194119421943194419451946194719481949195019511952195319541955195619571958195919601961196219631964196519661967196819691970197119721973197419751976197719781979198019811982198319841985198619871988198919901991199219931994199519961997199819992000200120022003200420052006200720082009201020112012201320142015201620172018201920202021202220232024202520262027202820292030203120322033203420352036203720382039204020412042204320442045204620472048204920502051205220532054205520562057205820592060206120622063206420652066206720682069207020712072207320742075207620772078207920802081208220832084208520862087208820892090209120922093209420952096209720982099210021012102210321042105210621072108210921102111211221132114211521162117211821192120212121222123212421252126212721282129213021312132213321342135213621372138213921402141214221432144214521462147214821492150215121522153215421552156215721582159216021612162216321642165216621672168216921702171217221732174217521762177217821792180218121822183218421852186218721882189219021912192219321942195219621972198219922002201220222032204220522062207220822092210221122122213221422152216221722182219222022212222222322242225222622272228222922302231223222332234223522362237223822392240224122422243224422452246224722482249225022512252225322542255225622572258225922602261226222632264226522662267226822692270227122722273227422752276227722782279228022812282228322842285228622872288228922902291229222932294229522962297229822992300230123022303230423052306230723082309231023112312231323142315231623172318231923202321232223232324232523262327232823292330233123322333233423352336233723382339234023412342234323442345234623472348234923502351235223532354235523562357235823592360236123622363236423652366236723682369237023712372237323742375237623772378237923802381238223832384238523862387238823892390239123922393239423952396239723982399240024012402240324042405240624072408240924102411241224132414241524162417241824192420242124222423242424252426242724282429243024312432243324342435243624372438243924402441244224432444244524462447244824492450245124522453245424552456245724582459246024612462246324642465246624672468246924702471247224732474247524762477247824792480248124822483248424852486248724882489249024912492249324942495249624972498249925002501250225032504250525062507250825092510251125122513251425152516251725182519252025212522252325242525252625272528252925302531253225332534253525362537253825392540254125422543254425452546254725482549255025512552255325542555255625572558255925602561256225632564256525662567256825692570257125722573257425752576257725782579258025812582258325842585258625872588258925902591259225932594259525962597259825992600260126022603260426052606260726082609261026112612261326142615261626172618261926202621262226232624262526262627262826292630263126322633263426352636263726382639264026412642264326442645264626472648264926502651265226532654265526562657265826592660266126622663266426652666266726682669267026712672267326742675267626772678267926802681268226832684268526862687268826892690269126922693269426952696269726982699270027012702270327042705270627072708270927102711271227132714271527162717271827192720272127222723272427252726272727282729273027312732273327342735273627372738273927402741274227432744274527462747274827492750275127522753275427552756275727582759276027612762276327642765276627672768276927702771277227732774277527762777277827792780278127822783278427852786278727882789279027912792279327942795279627972798279928002801280228032804280528062807280828092810281128122813281428152816281728182819282028212822282328242825282628272828282928302831283228332834283528362837283828392840284128422843284428452846284728482849285028512852285328542855285628572858285928602861286228632864286528662867286828692870287128722873287428752876287728782879288028812882288328842885288628872888288928902891289228932894289528962897289828992900290129022903290429052906290729082909291029112912291329142915291629172918291929202921292229232924292529262927292829292930293129322933293429352936293729382939294029412942294329442945294629472948294929502951295229532954295529562957295829592960296129622963296429652966296729682969297029712972297329742975297629772978297929802981298229832984298529862987298829892990299129922993299429952996299729982999300030013002300330043005300630073008300930103011301230133014301530163017301830193020302130223023302430253026302730283029303030313032303330343035303630373038303930403041304230433044304530463047304830493050305130523053305430553056305730583059306030613062306330643065306630673068306930703071307230733074307530763077307830793080308130823083308430853086308730883089309030913092309330943095309630973098309931003101310231033104310531063107310831093110311131123113311431153116311731183119312031213122312331243125312631273128312931303131313231333134313531363137313831393140314131423143314431453146314731483149315031513152315331543155315631573158315931603161316231633164316531663167316831693170317131723173317431753176317731783179318031813182318331843185318631873188318931903191319231933194319531963197319831993200 |
- #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();
- }
- }
-
-
- BigInteger(std::vector<char> val, int sign, int len) {
- int cursor = 0, numDigits;
-
- while (cursor < len && Character.digit(val[cursor], 10) == 0) {
- cursor++;
- }
- if (cursor == len) {
- signum = 0;
- mag = ZERO.mag;
- return;
- }
- numDigits = len - cursor;
- signum = sign;
-
- 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);
-
- int firstGroupLen = numDigits % digitsPerInt[10];
- if (firstGroupLen == 0)
- firstGroupLen = digitsPerInt[10];
- magnitude[numWords - 1] = parseInt(val, cursor, cursor += firstGroupLen);
-
- 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');
- }
-
-
-
- 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;
- }
-
-
-
-
- static void destructiveMulAdd(std::vector<int>& x, int y, int z) {
-
- 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;
- }
-
- 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);
- std::vector<char> randomBits(numBytes);
- std::uniform_int_distribution<char> dis(-128,127);
-
- 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;
- }
-
-
-
- static const int SMALL_PRIME_THRESHOLD = 95;
-
- static const int DEFAULT_PRIME_CERTAINTY = 100;
-
-
-
- bool primeToCertainty(int certainty, std::mt19937_64& random) {
- int rounds = 0;
- int n = (std::min(certainty, Integer.MAX_VALUE-1)+1)/2;
-
-
-
- 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);
-
- int d = 5;
- while (jacobiSymbol(d, this) != -1) {
-
- d = (d < 0) ? std::abs(d)+2 : -(d+2);
- }
-
- BigInteger u = lucasLehmerSequence(d, thisPlusOne, *this);
-
- return u.mod(*this).equals(ZERO);
- }
-
- static int jacobiSymbol(int p, BigInteger n) {
- if (p == 0)
- return 0;
-
- int j = 1;
- int u = n.mag[n.mag.size()-1];
-
- if (p < 0) {
- p = -p;
- int n8 = u & 7;
- if ((n8 == 3) || (n8 == 7))
- j = -j;
- }
-
- while ((p & 3) == 0)
- p >>= 2;
- if ((p & 1) == 0) {
- p >>= 1;
- if (((u ^ (u>>1)) & 2) != 0)
- j = -j;
- }
- if (p == 1)
- return j;
-
- if ((p & u & 2) != 0)
- j = -j;
-
- u = n.mod(BigInteger::valueOf(p)).intValue();
-
- while (u != 0) {
- while ((u & 3) == 0)
- u >>= 2;
- if ((u & 1) == 0) {
- u >>= 1;
- if (((p ^ (p>>1)) & 2) != 0)
- j = -j;
- }
- if (u == 1)
- return j;
-
- assert (u < p);
- int t = u; u = p; p = t;
- if ((u & p & 2) != 0)
- j = -j;
-
- 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) {
-
- BigInteger thisMinusOne = this->subtract(ONE);
- BigInteger m = thisMinusOne;
- int a = m.getLowestSetBit();
- m = m.shiftRight(a);
-
- for (int i=0; i < iterations; i++) {
-
- 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 BigInteger valueOf(std::int64_t val) {
-
- 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));
- }
-
-
- 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);
- }
-
- 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);
-
-
- 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;
- }
- }
-
- bool carry = (sum >>> 32 != 0);
- while (xIndex > 0 && carry)
- carry = ((result[--xIndex] = x[xIndex] + 1) == 0);
-
- while (xIndex > 0)
- result[--xIndex] = x[xIndex];
-
- 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> add(std::vector<int> x, std::vector<int> y) {
-
- if (x.size() < y.size()) {
- std::swap(x,y);
-
- }
- int xIndex = x.size();
- int yIndex = y.size();
- std::vector<int> result(xIndex);
- long sum = 0;
- if (yIndex == 1) {
- sum = (x[--xIndex] & LONG_MASK) + (y[0] & LONG_MASK) ;
- result[xIndex] = (int)sum;
- } else {
-
- while (yIndex > 0) {
- sum = (x[--xIndex] & LONG_MASK) +
- (y[--yIndex] & LONG_MASK) + (sum >>> 32);
- result[xIndex] = (int)sum;
- }
- }
-
- bool carry = (sum >>> 32 != 0);
- while (xIndex > 0 && carry)
- carry = ((result[--xIndex] = x[xIndex] + 1) == 0);
-
- while (xIndex > 0)
- result[--xIndex] = x[xIndex];
-
- 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;
-
- bool borrow = (difference >> 32 != 0);
- if (borrow) {
- result[0] = highWord - 1;
- } else {
- result[0] = highWord;
- }
- return result;
- } else {
- 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;
- }
-
- bool borrow = (difference >> 32 != 0);
- while (bigIndex > 0 && borrow)
- borrow = ((result[--bigIndex] = big[bigIndex] - 1) == -1);
-
- 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;
-
- while (littleIndex > 0) {
- difference = (big[--bigIndex] & LONG_MASK) -
- (little[--littleIndex] & LONG_MASK) +
- (difference >> 32);
- result[bigIndex] = (int)difference;
- }
-
- bool borrow = (difference >> 32 != 0);
- while (bigIndex > 0 && borrow)
- borrow = ((result[--bigIndex] = big[bigIndex] - 1) == -1);
-
- 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;
- uint64_t dl = v & LONG_MASK;
- 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();
-
- int half = (std::max(xlen, ylen)+1) / 2;
-
-
- 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);
- BigInteger p2 = xl.multiply(yl);
-
- BigInteger p3 = xh.add(xl).multiply(yh.add(yl));
-
- 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);
-
- int k = (largest+2)/3;
-
- int r = largest - 2*k;
-
-
- 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);
-
-
-
-
-
-
- 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);
-
- 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;
- }
-
-
- 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) {
- borrow = 1LL;
- } else {
- borrow = 0LL;
- }
-
-
-
- q = (w * 0xAAAAAAABL) & LONG_MASK;
- result[i] = (int) q;
-
-
- 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);
- }
-
-
- 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);
-
- 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) {
-
-
- 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;
- }
-
- 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);
- }
-
- 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();
- BigInteger xls = xl.square();
-
- 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();
-
- int k = (len+2)/3;
-
- int r = len - 2*k;
-
-
- 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();
-
-
-
-
-
-
- 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);
-
- int ss = k*32;
- return vinf.shiftLeft(ss).add(t2).shiftLeft(ss).add(t1).shiftLeft(ss).add(tm1).shiftLeft(ss).add(v0);
- }
-
-
- 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();
-
-
-
-
- int powersOfTwo = partToSquare.getLowestSetBit();
- long bitsToShift = (std::int64_t)powersOfTwo * exponent;
- if (bitsToShift > std::numeric_limits<int>::max()) {
- reportOverflow();
- }
- int remainingBits;
-
- if (powersOfTwo > 0) {
- partToSquare = partToSquare.shiftRight(powersOfTwo);
- remainingBits = partToSquare.bitLength();
- if (remainingBits == 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) {
- if (signum < 0 && (exponent&1) == 1) {
- return NEGATIVE_ONE;
- } else {
- return ONE;
- }
- }
- }
-
-
-
- long scaleFactor = (long)remainingBits * exponent;
-
-
- if (partToSquare.mag.size() == 1 && scaleFactor <= 62) {
-
- int newSign = (signum <0 && (exponent&1) == 1 ? -1 : 1);
- long result = 1;
- long baseToPow2 = partToSquare.mag[0] & LONG_MASK;
- int workingExponent = exponent;
-
- while (workingExponent != 0) {
- if ((workingExponent & 1) == 1) {
- result = result * baseToPow2;
- }
- if ((workingExponent >>>= 1) != 0) {
- baseToPow2 = baseToPow2 * baseToPow2;
- }
- }
-
- if (powersOfTwo > 0) {
- if (bitsToShift + scaleFactor <= 62) {
- return valueOf((result << bitsToShift) * newSign);
- } else {
- return valueOf(result*newSign).shiftLeft((int) bitsToShift);
- }
- }
- else {
- return valueOf(result*newSign);
- }
- } else {
-
-
-
- BigInteger answer = ONE;
- int workingExponent = exponent;
-
- while (workingExponent != 0) {
- if ((workingExponent & 1) == 1) {
- answer = answer.multiply(partToSquare);
- }
- if ((workingExponent >>>= 1) != 0) {
- partToSquare = partToSquare.square();
- }
- }
-
-
- 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 (n <= (32 - bitsInHighWord)) {
- primitiveLeftShift(a, len, nBits);
- return a;
- } else {
- 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;
- }
- }
- }
-
- 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;
- }
-
- 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);
- }
-
-
- 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)
-
- 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)) {
- result = base.oddModPow(exponent, m);
- } else {
-
-
- int p = m.getLowestSetBit();
- BigInteger m1 = m.shiftRight(p);
- BigInteger m2 = ONE.shiftLeft(p);
-
- BigInteger base2 = (this->signum < 0 || this->compareTo(m1) >= 0
- ? this->mod(m1) : this);
-
- BigInteger a1 = (m1.equals(ONE) ? ZERO :
- base2.oddModPow(exponent, m1));
-
- BigInteger a2 = base.modPow2(exponent, p);
-
- 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);
- }
-
-
-
-
-
- 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) {
-
- 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) {
-
- product = squareToLen(a, len, product);
- return montReduce(product, n, len, (int)inv);
- } else {
- return implMontgomerySquare(a, n, len, inv, materialize(product, len));
- }
- }
-
- 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);
- }
- }
-
-
-
- static std::vector<int> materialize(std::vector<int> z, int len) {
- if (z.size() < len)
- z = std::vector<int>(len);
- return z;
- }
-
-
- 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()};
-
- BigInteger oddModPow(BigInteger y, BigInteger z) {
-
-
- if (y.equals(ONE))
- return *this;
-
- 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();
-
-
-
-
- if ((modLen & 1) != 0) {
- std::vector<int> x(modlen + 1);
- std::copy(mod.begin(), mod.begin() + modLen, x.begin());
- mod = std::move(x);
- modLen++;
- }
-
- int wbits = 0;
- int ebits = bitLength(exp, exp.size());
-
- if ((ebits != 17) || (exp[0] != 65537)) {
- while (ebits > bnExpModThreshTable[wbits]) {
- wbits++;
- }
- }
-
- int tblmask = 1 << wbits;
-
- std::vector<std::vector<int>> table(tblmask);
- for (int i=0; i < tblmask; i++)
- table[i] = std::vector<int>(modLen);
-
-
- std::int64_t n0 = (mod[modLen-1] & LONG_MASK) + ((mod[modLen-2] & LONG_MASK) << 32);
- std::int64_t inv = -MutableBigInteger.inverseMod64(n0);
-
- std::vector<int> a = leftShift(base, base.size(), modLen << 5);
- MutableBigInteger q = MutableBigInteger(),
- a2 = MutableBigInteger(a),
- b2 = MutableBigInteger(mod);
- b2.normalize();
-
- MutableBigInteger r= a2.divide(b2, q);
- table[0] = r.toIntArray();
-
- if (table[0].size() < modLen) {
- int offset = modLen - table[0].size();
- std::vector<int> t2(modLen);
- std::copy(table[0].begin(),table[0].end(), t2.begin());
- table[0] = t2;
- }
-
- std::vector<int> b = montgomerySquare(table[0], mod, modLen, inv, null);
-
- std::vector<int> t(b.begin(), b.begin() + modLen);
-
- for (int i=1; i < tblmask; i++) {
- table[i] = montgomeryMultiply(t, table[i-1], mod, modLen, inv, null);
- }
-
- 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;
-
- 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;
-
- while (true) {
- ebits--;
-
- buf <<= 1;
- if (elen != 0) {
- buf |= ((exp[eIndex] & bitpos) != 0) ? 1 : 0;
- bitpos >>= 1;
- if (bitpos == 0) {
- eIndex++;
- bitpos = 1 << (32-1);
- elen--;
- }
- }
-
- if ((buf & tblmask) != 0) {
- multpos = ebits - wbits;
- while ((buf & 1) == 0) {
- buf >>= 1;
- multpos++;
- }
- mult = table[buf >> 1];
- buf = 0;
- }
-
- 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;
- }
- }
-
- if (ebits == 0)
- break;
-
- if (!isone) {
- t = b;
- a = montgomerySquare(t, mod, modLen, inv, a);
- t = a; a = b; b = t;
- }
- }
-
- 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;
- }
-
- 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) {
- return 1;
- } else {
- a[offset]++;
- if (a[offset] != 0)
- return 0;
- }
- }
- return 1;
- }
-
- BigInteger modPow2(BigInteger exponent, int p) {
-
- 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;
-
- int numInts = (p + 31) >> 5;
- std::vector<int> mag(numInts,0);
- System.arraycopy(this->mag, (this->mag.size() - numInts), mag, 0, numInts);
-
- 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;
-
- 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);
- }
-
-
- 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 {
-
-
- 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;
-
- 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) {
-
- 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;
- }
-
-
- 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);
- }
-
-
- 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) {
- lsb = 0;
- if (signum == 0) {
- lsb -= 1;
- } else {
-
- int i,b;
- for (i=0; (b = getInt(i)) == 0; i++)
- ;
- lsb += (i << 5) + __builtin_ctz(b);
- }
- lowestSetBit = lsb + 2;
- }
- return lsb;
- }
-
-
- int bitLength() {
- if (n == -1) {
- int[] m = mag;
- int len = m.size();
- if (len == 0) {
- n = 0;
- } else {
-
- int magBitLength = ((len - 1) << 5) + bitLengthForInt(mag[0]);
- if (signum < 0) {
-
- 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) {
- bc = 0;
-
- for (int i=0; i < mag.size(); i++)
- bc += __builtin_popcount(mag[i]);
- if (signum < 0) {
-
- 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;
- }
-
-
- 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);
- }
-
-
- 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;
- }
- }
-
-
-
- BigInteger min(BigInteger val) {
- return (compareTo(val) < 0 ? *this : val);
- }
-
- BigInteger max(BigInteger val) {
- return (compareTo(val) > 0 ? *this : val);
- }
-
-
- int hashCode() {
- int hashCode = 0;
- for (int i=0; i < mag.size(); i++)
- hashCode = (int)(31*hashCode + (mag[i] & LONG_MASK));
- return hashCode * signum;
- }
-
-
-
-
-
-
-
-
-
- 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;
- }
-
-
-
-
-
- static std::vector<int> stripLeadingZeroInts(std::vector<int> val) {
- int vlen = val.size();
- int keep;
-
- 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;
-
- 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;
-
- for (keep = 0; keep < byteLength && a[keep] == 0; keep++)
- ;
-
- 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();
-
- for (keep=0; keep < byteLength && a[keep] == -1; keep++);
-
- 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);
-
- 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);
-
- int mask = -1 >>> (8*(3-numBytesToTransfer));
- result[i] = ~result[i] & mask;
- }
-
- 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;
-
- for (keep=0; keep < a.size() && a[keep] == -1; keep++)
- ;
-
- 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);
-
- for (int i = keep; i < a.size(); i++)
- result[i - keep + extraInt] = ~a[i];
-
- for (int i = result.size()-1; ++result[i] == 0; i--);
- return result;
- }
-
- 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)};
-
- 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;
- }
-
- int signBit() {
- return signum < 0 ? 1 : 0;
- }
-
- 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) {
- fn = 0;
-
- int i;
- int mlen = mag.size();
- for (i = mlen - 1; i >= 0 && mag[i] == 0; i--);
- fn = mlen - i - 1;
- firstNonzeroIntNum = fn + 2;
- }
- 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);
- }
|