diff options
| -rw-r--r-- | bigint.cpp | 160 | ||||
| -rw-r--r-- | bigint.h | 2 | ||||
| -rw-r--r-- | numalgo.cpp | 1 |
3 files changed, 95 insertions, 68 deletions
@@ -177,9 +177,9 @@ Bigint& Bigint::operator=(slongdigit_t v){ return *this; } sign=v>=0?1:-1; - v*=sign; - digits[0]=v; - if(v>digits[0])digits.push_back(v>>digit_bits); + longdigit_t uv=sign*v; //unsigned version of v + digits.push_back(uv); + if(uv>digits[0])digits.push_back(uv>>digit_bits); checkconsistent(); return *this; } @@ -344,7 +344,7 @@ Bigint Bigint::operator>>(int sh) const { pair<Bigint,Bigint> Bigint::divmod(const Bigint &div) const { int bitcdiff=bitcount()-div.bitcount(); if(bitcdiff<0)bitcdiff=0; - pair<Bigint,Bigint> p=divmod(div,1,bitcdiff/29+10); //ignores all signs + pair<Bigint,Bigint> p=divmod(*this,div,bitcdiff/29+10); //ignores all signs /* To let the result come out correctly, we apply case analysis to the signs of the arguments. * As a guiding example, these two cases can be examined. * The code was tested with many large random numbers (also negative), using a Python script. @@ -381,71 +381,97 @@ pair<Bigint,Bigint> Bigint::divmod(const Bigint &div) const { //ignores all signs, and always returns positive numbers! //This function is opaque and way more complicated than it should be. Sorry for that. -//Strategy: find a quotient (and guess=quotient*div) such that guess<=*this, but *this-guess -//is small. Then subtract guess from *this, and repeat. -pair<Bigint,Bigint> Bigint::divmod(const Bigint &div,int depth,int maxdepth) const { - if(depth>maxdepth)assert(false); //Sanity check +//Strategy: find a quotient (and guess=quotient*div) such that guess<=a, but a-guess +//is small. Then subtract guess from a, and repeat. +pair<Bigint,Bigint> Bigint::divmod(Bigint a,const Bigint &div,const int maxiter){ if(div.digits.size()==0)throw domain_error("Bigint divide by zero"); - if(digits.size()==0)return make_pair(Bigint::zero,Bigint::zero); - - int cmp=compareAbs(div); - if(cmp==0)return make_pair(Bigint::one,Bigint::zero); - if(cmp<0){ - pair<Bigint,Bigint> p(make_pair(Bigint::zero,*this)); - p.second.sign=1; - return p; - } - //now *this is greater in magnitude than the divisor - - int thisbtc=bitcount(),divbtc=div.bitcount(); - assert(divbtc<=thisbtc); - Bigint quotient,guess; - if(thisbtc<=2*digit_bits){ - //simple integral division - longdigit_t thisnum=(digits.size()==2?((longdigit_t)1<<digit_bits)*digits[1]:0)+digits[0]; - longdigit_t divnum=(div.digits.size()==2?((longdigit_t)1<<digit_bits)*div.digits[1]:0)+div.digits[0]; - if(divnum==1)return make_pair(*this,Bigint::zero); - return make_pair( - Bigint((slongdigit_t)(thisnum/divnum)), - Bigint((slongdigit_t)(thisnum%divnum))); - } else if(divbtc>=digit_bits){ //both *this and div are large - //take 2*digit_bits of *this and 1*digit_bits of div; quotient gives a good guess - int spill=__builtin_clz(digits.back()); - longdigit_t thishead2=((longdigit_t)digits.back()<<(spill+digit_bits))|((longdigit_t)digits[digits.size()-2]<<spill); - if(spill>0)thishead2|=digits[digits.size()-3]>>(digit_bits-spill); - longdigit_t divhead=((longdigit_t)div.digits.back()<<digit_bits)|div.digits[div.digits.size()-2]; - divhead>>=digit_bits-__builtin_clz(div.digits.back()); - longdigit_t factor=thishead2/(divhead+1); //+1 to make sure the quotient guess is <= the actual quotient - quotient=factor; - quotient<<=thisbtc-digit_bits-divbtc; //shift amount may be negative if thisbtc and divbtc are less than digit_bits apart - if(quotient==0)quotient=1; //prevents against (HUGE+1)/HUGE where HUGE==HUGE - guess=quotient*div; - guess.sign=1; - } else { //divbtc<digit_bits, but *this is large - //take 2 digits of *this and all of div - int spill=__builtin_clz(digits.back()); - longdigit_t thishead2=((longdigit_t)digits.back()<<(spill+digit_bits))|((longdigit_t)digits[digits.size()-2]<<spill); - if(spill>0)thishead2|=digits[digits.size()-3]>>(digit_bits-spill); - longdigit_t factor=thishead2/div.digits[0]; - quotient=factor; - quotient<<=thisbtc-2*digit_bits; - guess=quotient*div; - guess.sign=1; - } - - //Now actually subtract out our guess - cmp=guess.compareAbs(*this); - assert(cmp<=0); //we specifically made our guessed quotient to be <= the actual quotient - if(cmp<0){ - Bigint rest(*this); - rest.sign=1; - rest-=guess; //also correct for *this and guess negative - pair<Bigint,Bigint> dm=rest.divmod(div,depth+1,maxdepth); - dm.first+=quotient; - return dm; - } else { //then cmp==0 - return make_pair(quotient,Bigint::zero); + + pair<Bigint,Bigint> result; + + a.sign=1; + + for(int iter=0;iter<maxiter;iter++){ //the maxiter is there to make sure we don't loop infinitely + if(a.digits.size()==0){ + result.second=0; + break; + } + + int cmp=a.compareAbs(div); + if(cmp==0){ + result.first+=1; + result.second=0; + break; + } + if(cmp<0){ + result.second=a; + break; + } + //now a is greater in magnitude than the divisor + + int abtc=a.bitcount(),divbtc=div.bitcount(); + assert(divbtc<=abtc); + Bigint quotient,guess; + if(abtc<=2*digit_bits){ + //simple integral division + longdigit_t anum=(a.digits.size()==2?((longdigit_t)1<<digit_bits)*a.digits[1]:0)+a.digits[0]; + longdigit_t divnum=(div.digits.size()==2?((longdigit_t)1<<digit_bits)*div.digits[1]:0)+div.digits[0]; + if(divnum==1){ + result.first+=a; + result.second=0; + break; + } + result.first+=(slongdigit_t)(anum/divnum); + result.second=(slongdigit_t)(anum%divnum); + break; + } else if(divbtc>=digit_bits){ //both a and div are large + //We're going to take 2*digit_bits of a and 1*digit_bits of div + int spill=__builtin_clz(a.digits.back()); //number of zero bits on top of a + longdigit_t ahead2= //top 64 bits of a (not yet) + ((longdigit_t)a.digits.back()<<(spill+digit_bits))| + ((longdigit_t)a.digits[a.digits.size()-2]<<spill); + if(spill>0){ //if there was some spill, we need to pull in a bit of the third digit + ahead2|=a.digits[a.digits.size()-3]>>(digit_bits-spill); + } + //now ahead is top 64 bits of a + + longdigit_t divhead= //top 2 digits of divisor + ((longdigit_t)div.digits.back()<<digit_bits)| + div.digits[div.digits.size()-2]; + divhead>>=digit_bits-__builtin_clz(div.digits.back()); //shift out some bits such that divhead is top 32 bits of div + + longdigit_t factor=ahead2/(divhead+1); //+1 to make sure the quotient guess is <= the actual quotient + quotient=factor; + quotient<<=abtc-digit_bits-divbtc; //shift amount may be negative if abtc and divbtc are less than digit_bits apart + if(quotient==0)quotient=1; //prevents against (HUGE+1)/HUGE where HUGE==HUGE + guess=quotient*div; + guess.sign=1; //for if div is negative + } else { //divbtc<digit_bits, but a is large + //We're going to take 2 digits of a and all of div (partly analogous to previous case) + int spill=__builtin_clz(a.digits.back()); + longdigit_t ahead2= + ((longdigit_t)a.digits.back()<<(spill+digit_bits))| + ((longdigit_t)a.digits[a.digits.size()-2]<<spill); + if(spill>0){ + ahead2|=a.digits[a.digits.size()-3]>>(digit_bits-spill); + } + + longdigit_t factor=ahead2/div.digits[0]; + quotient=factor; + quotient<<=abtc-2*digit_bits; + guess=quotient*div; + guess.sign=1; + } + + //Now actually subtract out our guess + a-=guess; + result.first+=quotient; + if(a==0){ + result.second=0; + break; + } } + + return result; } bool Bigint::operator==(const Bigint &o) const {return compare(o)==0;} @@ -32,7 +32,7 @@ private: //basically a debug check, but it can't hurt to have it since it takes few cpu cycles. void checkconsistent(); - std::pair<Bigint,Bigint> divmod(const Bigint&,int depth,int maxdepth) const; //ignores all signs + static std::pair<Bigint,Bigint> divmod(Bigint a,const Bigint &div,int maxdepth); //ignores all signs public: Bigint(); diff --git a/numalgo.cpp b/numalgo.cpp index ee77afc..a8f2321 100644 --- a/numalgo.cpp +++ b/numalgo.cpp @@ -40,6 +40,7 @@ Bigint expmod(const Bigint &b,const Bigint &e,const Bigint &m){ assert(e>=0); assert(m>=1); if(m==1)return Bigint::zero; + if(e==0)return b; Bigint res(1); vector<bool> bits(e.bits()); for(int i=bits.size()-1;i>=0;i--){ |