From d24ab714b958b9fece4631076e240739ad0dd23f Mon Sep 17 00:00:00 2001
From: tomsmeding <tom.smeding@gmail.com>
Date: Mon, 3 Oct 2016 22:15:25 +0200
Subject: Progress

---
 bigint.cpp | 106 +++++++++++++++++++++++++++++++++++++++++++++++++++----------
 1 file changed, 89 insertions(+), 17 deletions(-)

(limited to 'bigint.cpp')

diff --git a/bigint.cpp b/bigint.cpp
index 5182d07..644ffb3 100644
--- a/bigint.cpp
+++ b/bigint.cpp
@@ -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);
-- 
cgit v1.2.3-70-g09d2