diff options
Diffstat (limited to 'bigint.cpp')
| -rw-r--r-- | bigint.cpp | 110 |
1 files changed, 12 insertions, 98 deletions
@@ -63,7 +63,6 @@ void Bigint::add(Bigint &a,const Bigint &b){ longdigit_t bdig=i<(int)b.digits.size()?b.digits[i]:0; longdigit_t sum=a.digits[i]+bdig+carry; a.digits[i]=sum; - // carry=sum>=((longdigit_t)1<<digit_bits); carry=sum>>digit_bits; } if(carry)a.digits.push_back(1); @@ -88,12 +87,10 @@ void Bigint::subtract(Bigint &a,const Bigint &b){ digit_t adig=a.digits[i]; digit_t bdig=i<(int)b.digits.size()?b.digits[i]:0; digit_t res=adig-(bdig+carry); - // cerr<<"carry="<<carry<<" res="<<res<<" adig="<<adig<<" bdig="<<bdig<<endl; carry=(bdig||carry)&&res>=adig; a.digits[i]=res; } if(carry){ - // cerr<<"2s complement"<<endl; //we do a fake 2s complement, sort of carry=0; for(int i=0;i<sz;i++){ @@ -120,14 +117,12 @@ Bigint Bigint::product(const Bigint &a,const Bigint &b){ longdigit_t newd=pr+res.digits[i+j]; //this always fits, I checked res.digits[i+j]=(digit_t)newd; carry=newd>>digit_bits; - // cerr<<"carry="<<carry<<endl; } for(int j=bsz;carry;j++){ assert(i+j<(int)res.digits.size()); longdigit_t newd=res.digits[i+j]+carry; res.digits[i+j]=newd; carry=newd>>digit_bits; - // cerr<<"(2) carry="<<carry<<endl; } } res.sign=a.sign*b.sign; @@ -322,63 +317,24 @@ Bigint Bigint::operator>>(int sh) const { return Bigint(*this)>>=sh; } -int depthrecord; - pair<Bigint,Bigint> Bigint::divmod(const Bigint &div) const { // cerr<<"divmod("<<hex<<*this<<','<<div<<dec<<')'<<endl; int bitcdiff=bitcount()-div.bitcount(); if(bitcdiff<0)bitcdiff=0; - pair<Bigint,Bigint> res=divmod(div,1,bitcdiff/29+10); - //cerr<<hex<<*this<<' '<<div<<' '; - //cerr<<dec<<depthrecord<<endl; - //cerr<<" -> "<<hex<<res.first<<' '<<res.second<<dec<<endl; - return res; -} - -inline void record(int depth){ - depthrecord=depth; + return divmod(div,1,bitcdiff/29+10); } pair<Bigint,Bigint> Bigint::divmod(const Bigint &div,int depth,int maxdepth) const { if(depth>maxdepth)assert(false); // cerr<<"divmod("<<hex<<*this<<','<<div<<dec<<") depth="<<depth<<" maxdepth="<<maxdepth<<endl; if(div.digits.size()==0)throw domain_error("Bigint divide by zero"); - if(digits.size()==0){record(depth); return make_pair(Bigint::zero,Bigint::zero);} + if(digits.size()==0)return make_pair(Bigint::zero,Bigint::zero); int cmp=compareAbs(div); - if(cmp==0){record(depth); return make_pair(Bigint(sign*div.sign),Bigint::zero);} - if(cmp<0){record(depth); return make_pair(Bigint::zero,*this);} + if(cmp==0)return make_pair(Bigint(sign*div.sign),Bigint::zero); + if(cmp<0)return make_pair(Bigint::zero,*this); //now *this is greater in magnitude than the divisor -#if 0 - int twoexp=bitcount()-div.bitcount(); - if(twoexp<0)twoexp=0; - Bigint quotient(sign*div.sign); - quotient<<=twoexp; - Bigint guess(div); //guess == quotient * div - guess<<=twoexp; - //cerr<<"guess="<<hex<<guess<<endl; - //twoexp now becomes irrelevant -#else -#if 0 - Bigint guess,quotient; - if(digits.size()==div.digits.size()){ - quotient=digits.back()/div.digits.back(); - quotient.sign=sign*div.sign; - guess=quotient; - guess*=div; - } else { - assert(digits.size()>div.digits.size()); - assert(digits.size()>=2&&div.digits.size()>=1); - longdigit_t factor=((longdigit_t)1<<digit_bits)*digits.back()+digits[digits.size()-2]; - factor/=div.digits.back()+1; - quotient=factor; - quotient<<=(digits.size()-1-div.digits.size())*digit_bits; - quotient.sign=sign*div.sign; - guess=quotient*div; - } - cerr<<"guess="<<hex<<guess<<" quotient="<<quotient<<dec<<endl; -#else int thisbtc=bitcount(),divbtc=div.bitcount(); assert(divbtc<=thisbtc); Bigint quotient,guess; @@ -386,80 +342,47 @@ pair<Bigint,Bigint> Bigint::divmod(const Bigint &div,int depth,int maxdepth) con //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){record(depth); return make_pair(*this,Bigint::zero);} - record(depth); return make_pair( + if(divnum==1)return make_pair(*this,Bigint::zero); + return make_pair( Bigint(sign*div.sign*(slongdigit_t)(thisnum/divnum)), Bigint(sign*div.sign*(slongdigit_t)(thisnum%divnum))); } else if(divbtc>=digit_bits){ //the large case //take 2 digits of *this and 1 digit of div; quotient gives a good guess int spill=__builtin_clz(digits.back()); - //cerr<<"spill="<<spill<<endl; 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); - //cerr<<"thishead2="<<hex<<thishead2<<dec<<endl; longdigit_t divhead=((longdigit_t)div.digits.back()<<digit_bits)|div.digits[div.digits.size()-2]; divhead>>=digit_bits-__builtin_clz(div.digits.back()); - //cerr<<"divhead="<<hex<<divhead<<dec<<endl; longdigit_t factor=thishead2/(divhead+1); //+1 to make sure the quotient guess is <= the actual quotient - // cerr<<"factor="<<factor<<" thishead2="<<thishead2<<" divhead="<<divhead<<endl; 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 quotient.sign=sign*div.sign; guess=quotient*div; - // cerr<<"quotient = "<<hex<<quotient<<dec<<endl; - // cerr<<"guess = "<<hex<<guess<<dec<<endl; } 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]; - //cerr<<"factor="<<factor<<endl; quotient=factor; quotient<<=thisbtc-2*digit_bits; quotient.sign=sign*div.sign; - //cerr<<"quotient="<<hex<<quotient<<dec<<endl; guess=quotient*div; - //cerr<<"guess= "<<hex<<guess<<dec<<endl; } - //cerr<<"guess= "<<hex<<guess<<" quotient="<<quotient<<dec<<endl; -#endif -#endif + + //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 guess2(guess); - while(true){ - guess2<<=1; - int cmp=guess2.compareAbs(*this); - if(cmp>0)break; - guess<<=1; - quotient<<=1; - if(cmp==0){record(depth); return make_pair(quotient,Bigint::zero);} - }*/ Bigint rest(*this); 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); } - if(cmp==0){record(depth); return make_pair(quotient,Bigint::zero);} - //then cmp>0, so our guess is too large - do { - if(quotient.digits[0]&1){ - guess-=div; - // quotient-=Bigint::one; // not necessary, since we shift the bit out anyway - } - guess>>=1; - quotient>>=1; - cmp=guess.compareAbs(*this); - } while(cmp>0); - if(cmp==0){record(depth); return make_pair(quotient,Bigint::zero);} - Bigint rest(*this); - rest-=guess; - pair<Bigint,Bigint> dm=rest.divmod(div,depth+1,maxdepth); - dm.first+=quotient; - return dm; } bool Bigint::operator==(const Bigint &o) const {return compare(o)==0;} @@ -619,7 +542,6 @@ istream& operator>>(istream &is,Bigint &b){ if(!is)break; b*=ten; b+=c-'0'; - // cerr<<"b="<<b<<endl; } if(!acted)is.setstate(ios_base::failbit); b.normalise(); @@ -641,14 +563,7 @@ std::ostream& operator<<(std::ostream &os,Bigint b){ } return os; } -#if 0 - assert(b.digits.size()<=2); - if(b.sign==-1)os<<'-'; - if(b.digits.size()==2)os<<((uint64_t)b.digits[1]<<32)+b.digits[0]; - else if(b.digits.size()==1)os<<b.digits[0]; - else os<<'0'; - return os; -#else + if(b==0)return os<<'0'; Bigint div((int64_t)1000000000000000000LL); vector<Bigint::longdigit_t> outbuf; @@ -667,5 +582,4 @@ std::ostream& operator<<(std::ostream &os,Bigint b){ (i==(int)outbuf.size()-1?os:os<<setfill('0')<<setw(18))<<outbuf[i]; } return os; -#endif } |