divmod: recursive -> iterative

3 files changed, 95 insertions, 68 deletions

diff --git a/bigint.cpp b/bigint.cpp
index 7c9fec7..e57752f 100644
--- a/bigint.cpp
+++ b/bigint.cpp
@@ -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;}
diff --git a/bigint.h b/bigint.h
index cff7bd4..c66dcf6 100644
--- a/bigint.h
+++ b/bigint.h
@@ -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--){