From 27472e604eeb74b3b60313d38d537c7e0e83151b Mon Sep 17 00:00:00 2001
From: tomsmeding <tom.smeding@gmail.com>
Date: Sat, 8 Oct 2016 15:07:23 +0200
Subject: Fix sign issues in divmod

---
 bigint.cpp    | 62 +++++++++++++++++++++++++++++++++++++++++++++++++----------
 bigint.h      |  8 ++++++--
 biginttest.py | 13 +++++++------
 3 files changed, 65 insertions(+), 18 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;
diff --git a/bigint.h b/bigint.h
index 02f3da3..1e8515a 100644
--- a/bigint.h
+++ b/bigint.h
@@ -26,7 +26,7 @@ private:
 	void normalise();
 	void checkconsistent();
 
-	std::pair<Bigint,Bigint> divmod(const Bigint&,int depth,int maxdepth) const;
+	std::pair<Bigint,Bigint> divmod(const Bigint&,int depth,int maxdepth) const; //ignores all signs
 
 public:
 	Bigint();
@@ -63,7 +63,11 @@ public:
 	Bigint operator*(slongdigit_t) const;
 	Bigint operator<<(int) const;
 	Bigint operator>>(int) const;
-	std::pair<Bigint,Bigint> divmod(const Bigint&) const; //rounds towards zero; returns {quotient,remainder}
+
+	//Uses *mathematical* division-with-remainder.
+	//If we look at {q,r} = x.divmod(y), then x = q*y + r, and 0<=r<|q|. This is *not* compatible with many other
+	//programming languages (and, for that matter, the C++ built-in behaviour).
+	std::pair<Bigint,Bigint> divmod(const Bigint&) const;
 
 	bool operator==(const Bigint&) const;
 	bool operator!=(const Bigint&) const;
diff --git a/biginttest.py b/biginttest.py
index f125953..71c5417 100755
--- a/biginttest.py
+++ b/biginttest.py
@@ -11,7 +11,7 @@ def check(desc,x,y):
 
 def gendata():
 	for _ in range(ntimes):
-		yield random.randint(0,maxn), random.randint(1,maxn)
+		yield random.randint(-maxn,maxn), random.randint(-maxn,maxn)
 
 def proctest():
 	proc=subprocess.Popen(["./main"],stdin=subprocess.PIPE,stdout=subprocess.PIPE,stderr=sys.stderr)
@@ -21,11 +21,12 @@ def proctest():
 		proc.stdin.write("mod {} {}\n".format(a,b).encode("ascii"))
 		proc.stdin.flush()
 
-		ans=int(proc.stdout.readline())
-		check("{}/{}".format(a,b),ans,a//b)
-
-		ans=int(proc.stdout.readline())
-		check("{}%{}".format(a,b),ans,a%b)
+		q=int(proc.stdout.readline())
+		r=int(proc.stdout.readline())
+		if r<0 or r>=abs(b) or a!=q*b+r:
+			print("Error: {} divmod {}".format(a,b))
+			print("Diff: {}".format(a-q*b-r))
+			sys.exit(1)
 
 	proc.kill()
 
-- 
cgit v1.2.3-70-g09d2