diff options
author | tomsmeding <tom.smeding@gmail.com> | 2016-10-03 22:15:25 +0200 |
---|---|---|
committer | tomsmeding <tom.smeding@gmail.com> | 2016-10-03 22:15:25 +0200 |
commit | d24ab714b958b9fece4631076e240739ad0dd23f (patch) | |
tree | 7d4d35dc7a7b2a03863b9a5bbde9382e14b7d1aa /bigint.cpp | |
parent | 2bf5effe95641667a1ed51c04eff7760f6a42ef4 (diff) |
Progress
Diffstat (limited to 'bigint.cpp')
-rw-r--r-- | bigint.cpp | 106 |
1 files changed, 89 insertions, 17 deletions
@@ -7,6 +7,10 @@ using namespace std; +Bigint Bigint::zero(0); +Bigint Bigint::one(1); +Bigint Bigint::mone(-1); + Bigint::Bigint() :sign(1){} @@ -163,6 +167,11 @@ Bigint& Bigint::operator*=(const Bigint &o){ return *this; } +//TODO: optimise these functions +Bigint& Bigint::operator+=(slongdigit_t n){return *this+=Bigint(n);} +Bigint& Bigint::operator-=(slongdigit_t n){return *this-=Bigint(n);} +Bigint& Bigint::operator*=(slongdigit_t n){return *this*=Bigint(n);} + Bigint& Bigint::operator<<=(int sh){ if(sh==0)return *this; if(digits.size()==0)return *this; @@ -233,12 +242,26 @@ Bigint Bigint::operator*(const Bigint &o) const { return product(*this,o); } +//TODO: optimise these functions +Bigint Bigint::operator+(slongdigit_t n) const {return *this+Bigint(n);} +Bigint Bigint::operator-(slongdigit_t n) const {return *this-Bigint(n);} +Bigint Bigint::operator*(slongdigit_t n) const {return *this*Bigint(n);} + +Bigint Bigint::operator<<(int sh) const { + return Bigint(*this)<<=sh; +} + +Bigint Bigint::operator>>(int sh) const { + return Bigint(*this)>>=sh; +} + int depthrecord; pair<Bigint,Bigint> Bigint::divmod(const Bigint &div) const { pair<Bigint,Bigint> res=divmod(div,1); //cerr<<hex<<*this<<' '<<div<<' '; //cerr<<dec<<depthrecord<<endl; + //cerr<<" -> "<<hex<<res.first<<' '<<res.second<<dec<<endl; return res; } @@ -247,13 +270,14 @@ inline void record(int depth){ } pair<Bigint,Bigint> Bigint::divmod(const Bigint &div,int depth) const { - //cerr<<"divmod("<<hex<<*this<<','<<hex<<div<<')'<<endl; + if(depth>100)assert(false); + //cerr<<"divmod("<<hex<<*this<<','<<hex<<div<<dec<<") depth="<<depth<<endl; if(div.digits.size()==0)throw domain_error("Bigint divide by zero"); - if(digits.size()==0){record(depth); return make_pair(Bigint(0),Bigint(0));} + if(digits.size()==0){record(depth); return make_pair(Bigint::zero,Bigint::zero);} int cmp=compareAbs(div); - if(cmp==0){record(depth); return make_pair(Bigint(sign*div.sign),Bigint(0));} - if(cmp<0){record(depth); return make_pair(Bigint(0),*this);} + 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);} //now *this is greater in magnitude than the divisor #if 0 @@ -266,9 +290,11 @@ pair<Bigint,Bigint> Bigint::divmod(const Bigint &div,int depth) const { //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 { @@ -277,42 +303,81 @@ pair<Bigint,Bigint> Bigint::divmod(const Bigint &div,int depth) const { longdigit_t factor=((longdigit_t)1<<digit_bits)*digits.back()+digits[digits.size()-2]; factor/=div.digits.back()+1; quotient=factor; - quotient<<=(digits.size()-div.digits.size()-1)*digit_bits; - guess=quotient; - guess*=div; + quotient<<=(digits.size()-1-div.digits.size())*digit_bits; + quotient.sign=sign*div.sign; + guess=quotient*div; } - // cerr<<"guess="<<hex<<guess<<" quotient="<<quotient<<endl; + cerr<<"guess="<<hex<<guess<<" quotient="<<quotient<<dec<<endl; +#else + 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){record(depth); return make_pair(*this,Bigint::zero);} + record(depth); 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 + quotient=factor; + quotient<<=thisbtc-digit_bits-divbtc; //shift amount may be negative if thisbtc and divbtc < digit_bits apart + quotient.sign=sign*div.sign; + guess=quotient*div; + } 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))|(digits[digits.size()-2]<<spill); + if(spill>0)thishead2|=digits[digits.size()-3]>>(digit_bits-spill); + longdigit_t factor=thishead2/(div.digits.back()+1); //+1 to make sure the quotient guess is <= the actual quotient + quotient=factor; + quotient<<=thisbtc-2*digit_bits; + quotient.sign=sign*div.sign; + guess=quotient*div; + } + //cerr<<"guess= "<<hex<<guess<<" quotient="<<quotient<<dec<<endl; +#endif #endif cmp=guess.compareAbs(*this); if(cmp<0){ - Bigint guess2(guess); + /*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(0));} - } + 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); dm.first+=quotient; return dm; } - if(cmp==0){record(depth); return make_pair(quotient,Bigint(0));} + if(cmp==0){record(depth); return make_pair(quotient,Bigint::zero);} //then cmp>0, so our guess is too large - Bigint one(1); do { if(quotient.digits[0]&1){ guess-=div; - // quotient-=one; // not necessary, since we shift the bit out anyway + // 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(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); @@ -372,7 +437,7 @@ int Bigint::compareAbs(slongdigit_t v) const { int Bigint::bitcount() const { if(digits.size()==0)return 0; - return (digits.size()-1)*8+ilog2(digits.back())+1; + return (digits.size()-1)*digit_bits+ilog2(digits.back())+1; } Bigint::slongdigit_t Bigint::lowdigits() const { @@ -382,6 +447,13 @@ Bigint::slongdigit_t Bigint::lowdigits() const { return ((slongdigit_t)1<<digit_bits)*(digits[1]&mask)+digits[0]; } +bool Bigint::even() const { + return digits.size()==0||(digits[0]&1)==0; +} +bool Bigint::odd() const { + return !even(); +} + vector<char> Bigint::serialise() const { vector<char> v(1+digits.size()*sizeof(digit_t)); v[0]=sign; @@ -446,7 +518,7 @@ istream& operator>>(istream &is,Bigint &b){ return is; } b*=ten; - b+=Bigint(c-'0'); + b+=c-'0'; // cerr<<"b="<<b<<endl; } if(!acted)is.setstate(ios_base::failbit); |