Fix sign issues in divmod

1 files changed, 52 insertions, 10 deletions

diff --git a/bigint.cpp b/bigint.cpp
index a9763dd..639d263 100644
--- a/bigint.cpp
+++ b/bigint.cpp
@@ -321,9 +321,42 @@ 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;
-	return divmod(div,1,bitcdiff/29+10);
+	pair<Bigint,Bigint> p=divmod(div,1,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.
+	 * (1)  4 =  1* 3 + 1    6 =  2* 3 + 0
+	 * (2)  4 = -1*-3 + 1    6 = -2*-3 + 0
+	 * (3) -4 = -2* 3 + 2   -6 = -2* 3 + 0
+	 * (4) -4 =  2*-3 + 2   -6 =  2*-3 + 0
+	 */
+	if(sign==1){
+		if(div.sign==1){ // (1)
+			//nothing to do
+		} else { // (2)
+			p.first.sign=-1;
+		}
+	} else {
+		if(div.sign==1){ // (3)
+			p.first.sign=-1;
+			if(p.second!=0){
+				p.first-=1;
+				p.second=div-p.second;
+			}
+		} else { // (4)
+			if(p.second!=0){
+				p.first+=1;
+				p.second.sign=-1;
+				p.second-=div;
+			}
+		}
+	}
+	p.first.normalise();
+	p.second.normalise();
+	return p;
 }
 
+//ignores all signs, and always returns positive numbers!
 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;
@@ -331,8 +364,12 @@ pair<Bigint,Bigint> Bigint::divmod(const Bigint &div,int depth,int maxdepth) con
 	if(digits.size()==0)return make_pair(Bigint::zero,Bigint::zero);
 
 	int cmp=compareAbs(div);
-	if(cmp==0)return make_pair(Bigint(sign*div.sign),Bigint::zero);
-	if(cmp<0)return make_pair(Bigint::zero,*this);
+	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();
@@ -344,10 +381,10 @@ pair<Bigint,Bigint> Bigint::divmod(const Bigint &div,int depth,int maxdepth) con
 		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(sign*div.sign*(slongdigit_t)(thisnum/divnum)),
-			Bigint(sign*div.sign*(slongdigit_t)(thisnum%divnum)));
+			Bigint((slongdigit_t)(thisnum/divnum)),
+			Bigint((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
+		//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);
@@ -357,8 +394,8 @@ pair<Bigint,Bigint> Bigint::divmod(const Bigint &div,int depth,int maxdepth) con
 		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;
+		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());
@@ -367,8 +404,8 @@ pair<Bigint,Bigint> Bigint::divmod(const Bigint &div,int depth,int maxdepth) con
 		longdigit_t factor=thishead2/div.digits[0];
 		quotient=factor;
 		quotient<<=thisbtc-2*digit_bits;
-		quotient.sign=sign*div.sign;
 		guess=quotient*div;
+		guess.sign=1;
 	}
 
 	//Now actually subtract out our guess
@@ -376,6 +413,7 @@ pair<Bigint,Bigint> Bigint::divmod(const Bigint &div,int depth,int maxdepth) con
 	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;
@@ -504,10 +542,12 @@ istream& operator>>(istream &is,Bigint &b){
 	while(isspace(is.peek()))is.get();
 	if(!is)return is;
 	b.digits.resize(0);
+	bool negative=false;
 	if(is.peek()=='-'){
-		b.sign=-1;
+		negative=true;
 		is.get();
-	} else b.sign=1;
+	}
+	b.sign=1;
 	bool acted=false;
 	if(is.peek()=='0'){
 		is.get();
@@ -528,6 +568,7 @@ istream& operator>>(istream &is,Bigint &b){
 				b+=n;
 			}
 			if(!acted)is.setstate(ios_base::failbit);
+			else if(negative)b.sign=-1;
 			b.normalise();
 			b.checkconsistent();
 			return is;
@@ -544,6 +585,7 @@ istream& operator>>(istream &is,Bigint &b){
 		b+=c-'0';
 	}
 	if(!acted)is.setstate(ios_base::failbit);
+	else if(negative)b.sign=-1;
 	b.normalise();
 	b.checkconsistent();
 	return is;