path: root/bigint.cpp

#include <iomanip>
#include <sstream>
#include <stdexcept>
#include <cctype>
#include <cassert>
#include "bigint.h"
#include "numalgo.h"

using namespace std;

Bigint Bigint::mone(-1);
Bigint Bigint::zero(0);
Bigint Bigint::one(1);
Bigint Bigint::two(2);

Bigint::Bigint()
		:sign(1){}

Bigint::Bigint(Bigint &&o)
		:digits(move(o.digits)),sign(o.sign){
	o.sign=1;
	checkconsistent();
	o.checkconsistent();
}

Bigint::Bigint(const string &repr){
	stringstream(repr)>>*this;
	checkconsistent();
}

Bigint::Bigint(slongdigit_t v)
		:digits(1,abs(v)),sign(v>=0?1:-1){
	static_assert(sizeof(longdigit_t)==2*sizeof(digit_t),
		"longdigit_t should be twice as large as digit_t");
	v=abs(v);
	if(v>digits[0])digits.push_back(v>>digit_bits);
	else if(v==0){
		digits.clear();
		sign=1;
	}
	checkconsistent();
}

Bigint::Bigint(longdigit_t v)
		:digits(1,(digit_t)v),sign(1){
	if(v>digits[0])digits.push_back(v>>digit_bits);
	else if(v==0)digits.clear();
	checkconsistent();
}

Bigint::Bigint(sdigit_t v)
		:digits(1,abs(v)),sign(v>=0?1:-1){
	if(v==0){
		digits.clear();
		sign=1;
	}
	checkconsistent();
}

Bigint::Bigint(digit_t v)
		:digits(1,v),sign(1){
	if(v==0)digits.clear();
	checkconsistent();
}

//ignores sign of arguments
void Bigint::add(Bigint &a,const Bigint &b){
	if(a.digits.size()<b.digits.size())a.digits.resize(b.digits.size());
	int sz=a.digits.size();
	int carry=0;
	for(int i=0;i<sz;i++){
		longdigit_t bdig=i<(int)b.digits.size()?b.digits[i]:0;
		longdigit_t sum=a.digits[i]+bdig+carry;
		a.digits[i]=sum;
		carry=sum>>digit_bits;
	}
	if(carry)a.digits.push_back(1);
	a.normalise();
	a.checkconsistent();
}

//ignores sign of arguments
void Bigint::subtract(Bigint &a,const Bigint &b){
	if(a.digits.size()<b.digits.size()){
		a.digits.resize(b.digits.size()); //adds zeros
	}
	assert(a.digits.size()>=b.digits.size());
	if(a.digits.size()==0){ //then a==b==0
		a.checkconsistent();
		return;
	}
	assert(a.digits.size()>0);
	int sz=a.digits.size();
	int carry=0;
	for(int i=0;i<sz;i++){
		if(i>=(int)b.digits.size()&&!carry)break;
		digit_t adig=a.digits[i];
		digit_t bdig=i<(int)b.digits.size()?b.digits[i]:0;
		digit_t res=adig-(bdig+carry);
		carry=(bdig||carry)&&res>=adig;
		a.digits[i]=res;
	}
	if(carry){ //Apparently, a<b
		//we do a fake 2s complement, sort of
		carry=0;
		for(int i=0;i<sz;i++){
			a.digits[i]=~a.digits[i];
			a.digits[i]+=(i==0)+carry;
			carry=a.digits[i]<=(digit_t)carry;
		}
		a.sign=-a.sign;
	}
	a.shrink();
	a.normalise();
	a.checkconsistent();
}

//This is simple O(n^2) multiplication, with no optimisation for a==b. Also, no Karatsuba here.
//This works and is simple (KISS), but doesn't perform that well. In fact, it seems to be the
//bottleneck of the entire bigint implementation (through its use divmod, which is used A LOT
//for modulo in RSA).
Bigint Bigint::product(const Bigint &a,const Bigint &b){
	int asz=a.digits.size(),bsz=b.digits.size();
	if(asz==0||bsz==0)return Bigint();
	Bigint res;
	res.digits.resize(asz+bsz);
	for(int i=0;i<asz;i++){
		digit_t carry=0;
		for(int j=0;j<bsz;j++){
			longdigit_t pr=(longdigit_t)a.digits[i]*b.digits[j]+carry;
			longdigit_t newd=pr+res.digits[i+j]; //this always fits, I checked
			res.digits[i+j]=(digit_t)newd;
			carry=newd>>digit_bits;
		}
		for(int j=bsz;carry;j++){
			assert(i+j<(int)res.digits.size());
			longdigit_t newd=res.digits[i+j]+carry;
			res.digits[i+j]=newd;
			carry=newd>>digit_bits;
		}
	}
	res.sign=a.sign*b.sign;
	res.shrink();
	res.normalise();
	res.checkconsistent();
	return res;
}

void Bigint::shrink(){
	while(digits.size()&&digits.back()==0)digits.pop_back();
}

void Bigint::normalise(){
	if(digits.size()==0&&sign==-1)sign=1;
}

void Bigint::checkconsistent(){
	assert(digits.size()==0||digits.back()!=0);
	assert(digits.size()!=0||sign==1);
}

Bigint& Bigint::operator=(Bigint &&o){
	sign=o.sign;
	digits=move(o.digits);
	o.sign=1;
	normalise();
	checkconsistent();
	o.checkconsistent();
	return *this;
}

Bigint& Bigint::operator=(slongdigit_t v){
	digits.clear();
	if(v==0){
		sign=1;
		checkconsistent();
		return *this;
	}
	sign=v>=0?1:-1;
	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;
}

Bigint& Bigint::operator=(longdigit_t v){
	digits.clear();
	if(v!=0){
		digits.push_back((digit_t)v);
		if(v>digits[0])digits.push_back(v>>digit_bits);
	}
	checkconsistent();
	return *this;
}

Bigint& Bigint::operator=(sdigit_t v){
	digits.clear();
	if(v==0)sign=1;
	else {
		sign=v>=0?1:-1;
		v*=sign;
		digits.push_back(v);
	}
	checkconsistent();
	return *this;
}

Bigint& Bigint::operator=(digit_t v){
	digits.clear();
	sign=1;
	if(v!=0)digits.push_back(v);
	checkconsistent();
	return *this;
}

Bigint& Bigint::operator+=(const Bigint &o){
	if(&o==this){ //*this + *this = *this<<1
		operator<<=(1);
		return *this;
	}
	if(sign==1){
		if(o.sign==1)add(*this,o);
		else subtract(*this,o);
	} else {
		if(o.sign==1)subtract(*this,o);
		else add(*this,o);
	}
	checkconsistent();
	return *this;
}

Bigint& Bigint::operator-=(const Bigint &o){
	if(&o==this){ // *this - *this = 0
		sign=1;
		digits.clear();
		return *this;
	}
	if(sign==1){
		if(o.sign==1)subtract(*this,o);
		else add(*this,o);
	} else {
		if(o.sign==1)add(*this,o);
		else subtract(*this,o);
	}
	checkconsistent();
	return *this;
}

Bigint& Bigint::operator*=(const Bigint &o){
	*this=product(*this,o);
	checkconsistent();
	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;
	if(sh<0)return *this>>=-sh; //we support negative shifting
	if(sh/digit_bits>0){ //first shift by a multiple of our digit size, since that's easy
		digits.insert(digits.begin(),sh/digit_bits,0);
		sh%=digit_bits;
		if(sh==0){
			checkconsistent();
			return *this;
		}
	}
	digits.push_back(0); //afterwards, shift by the remaining amount
	for(int i=digits.size()-2;i>=0;i--){
		digits[i+1]|=digits[i]>>(digit_bits-sh);
		digits[i]<<=sh;
	}
	shrink();
	normalise();
	checkconsistent();
	return *this;
}

Bigint& Bigint::operator>>=(int sh){
	if(sh==0)return *this;
	if(digits.size()==0)return *this;
	if(sh<0)return *this<<=-sh; //we support negative shifting
	if(sh/digit_bits>0){ //first shift by a multiple of our digit size, since that's easy
		if(sh/digit_bits>=(int)digits.size()){
			digits.clear();
			sign=1;
			checkconsistent();
			return *this;
		}
		digits.erase(digits.begin(),digits.begin()+sh/digit_bits);
		sh%=digit_bits;
		if(sh==0){
			checkconsistent();
			return *this;
		}
	}
	digits[0]>>=sh; //afterwards, shift by the remaining amount
	int sz=digits.size();
	for(int i=1;i<sz;i++){
		digits[i-1]|=digits[i]<<(digit_bits-sh);
		digits[i]>>=sh;
	}
	shrink();
	normalise();
	checkconsistent();
	return *this;
}

Bigint& Bigint::negate(){
	sign=-sign;
	return *this;
}

Bigint Bigint::operator+(const Bigint &o) const {
	return Bigint(*this)+=o;
}

Bigint Bigint::operator-(const Bigint &o) const {
	return Bigint(*this)-=o;
}

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;
}

pair<Bigint,Bigint> Bigint::divmod(const Bigint &div) const {
	int bitcdiff=bitcount()-div.bitcount();
	if(bitcdiff<0)bitcdiff=0;
	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.
	 * (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!
//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<=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");

	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;}
bool Bigint::operator!=(const Bigint &o) const {return compare(o)!=0;}
bool Bigint::operator<(const Bigint &o) const {return compare(o)<0;}
bool Bigint::operator>(const Bigint &o) const {return compare(o)>0;}
bool Bigint::operator<=(const Bigint &o) const {return compare(o)<=0;}
bool Bigint::operator>=(const Bigint &o) const {return compare(o)>=0;}

bool Bigint::operator==(slongdigit_t v) const {return compare(v)==0;}
bool Bigint::operator!=(slongdigit_t v) const {return compare(v)!=0;}
bool Bigint::operator<(slongdigit_t v) const {return compare(v)<0;}
bool Bigint::operator>(slongdigit_t v) const {return compare(v)>0;}
bool Bigint::operator<=(slongdigit_t v) const {return compare(v)<=0;}
bool Bigint::operator>=(slongdigit_t v) const {return compare(v)>=0;}

int Bigint::compare(const Bigint &o) const {
	if(sign>o.sign)return 1;
	if(sign<o.sign)return -1;
	return sign*compareAbs(o);
}

int Bigint::compare(slongdigit_t v) const {
	if(sign==-1&&v>=0)return -1;
	if(sign==1&&v<0)return 1;
	return sign*compareAbs(v);
}

int Bigint::compareAbs(const Bigint &o) const {
	int sz=digits.size(),osz=o.digits.size();
	if(sz>osz)return 1;
	if(sz<osz)return -1;
	for(int i=sz-1;i>=0;i--){
		if(digits[i]>o.digits[i])return 1;
		if(digits[i]<o.digits[i])return -1;
	}
	return 0;
}

int Bigint::compareAbs(slongdigit_t v) const {
	v=abs(v);
	if(digits.size()>2)return 1;
	if(digits.size()==0)return v==0?0:-1;
	if(digits.size()==2){
		if(digits[1]>(digit_t)(v>>digit_bits))return 1;
		if(digits[1]<(digit_t)(v>>digit_bits))return -1;
	}
	if(digits[0]<(digit_t)v)return -1;
	if(digits[0]>(digit_t)v)return 1;
	return 0;
}

int Bigint::bitcount() const {
	if(digits.size()==0)return 0;
	return (digits.size()-1)*digit_bits+ilog2(digits.back())+1;
}

Bigint::slongdigit_t Bigint::lowdigits() const {
	if(digits.size()==0)return 0;
	if(digits.size()==1)return digits[0];
	longdigit_t mask=~((longdigit_t)1<<(digit_bits-1));
	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();
}

//Produces a string with the bytes of the mantissa in little-endian order.
string Bigint::serialiseMantissa() const {
	string s;
	s.resize(digits.size()*sizeof(digit_t));
	int sz=digits.size();
	for(int i=0;i<sz;i++){
		for(int j=0;j<(int)sizeof(digit_t);j++){
			s[i*sizeof(digit_t)+j]=(digits[i]>>(8*j))&0xff;
		}
	}
	return s;
}

//Inverse of serialiseMantissa
void Bigint::deserialiseMantissa(const string &s){
	if(s.size()%sizeof(digit_t)!=0)throw invalid_argument("Not a serialised Bigint");
	sign=1;
	int sz=s.size()/sizeof(digit_t);
	digits.resize(sz);
	for(int i=0;i<sz;i++){
		digits[i]=0;
		for(int j=0;j<(int)sizeof(digit_t);j++){
			digits[i]|=(uint8_t)s[i*sizeof(digit_t)+j]<<(8*j);
		}
	}
	shrink();
	normalise();
	checkconsistent();
}

vector<bool> Bigint::bits() const {
	if(digits.size()==0)return {};
	vector<bool> v(digit_bits*(digits.size()-1)+ilog2(digits.back())+1);
	int sz=digits.size();
	for(int i=0;i<sz;i++){
		digit_t dig=digits[i];
		for(int j=0;dig;j++){
			if(dig&1)v[digit_bits*i+j]=true;
			dig>>=1;
		}
	}
	return v;
}

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()=='-'){
		negative=true;
		is.get();
	}
	b.sign=1;
	bool acted=false;
	if(is.peek()=='0'){
		is.get();
		if(is.peek()=='x'){ //hex value
			is.get();
			acted=false;
			while(true){
				char c=is.peek();
				if(!isdigit(c)&&(c<'a'||c>'f')&&(c<'A'||c>'F'))break;
				acted=true;
				is.get();
				if(!is)break;
				int n;
				if(c<='9')n=c-'0';
				else if(c<='F')n=c-'A'+10;
				else n=c-'a'+10;
				b<<=4;
				b+=n;
			}
			if(!acted)is.setstate(ios_base::failbit);
			else if(negative)b.sign=-1;
			b.normalise();
			b.checkconsistent();
			return is;
		} else acted=true;
	}
	Bigint ten(10);
	while(true){
		char c=is.peek();
		if(!isdigit(c))break;
		acted=true;
		is.get();
		if(!is)break;
		b*=ten;
		b+=c-'0';
	}
	if(!acted)is.setstate(ios_base::failbit);
	else if(negative)b.sign=-1;
	b.normalise();
	b.checkconsistent();
	return is;
}

std::ostream& operator<<(std::ostream &os,Bigint b){
	if(b<0){
		os<<'-';
		b.negate();
	}
	if(os.flags()&ios_base::hex){
		os<<"0x";
		if(b.digits.size()==0)return os<<'0';
		os<<b.digits.back();
		for(int i=b.digits.size()-2;i>=0;i--){
			os<<setw(Bigint::digit_bits/4)<<setfill('0')<<b.digits[i];
		}
		return os;
	}

	if(b==0)return os<<'0';
	Bigint div((int64_t)1000000000000000000LL);
	vector<Bigint::longdigit_t> outbuf;
	while(b!=0){
		pair<Bigint,Bigint> dm=b.divmod(div);
		b=dm.first;
		Bigint::longdigit_t val=0;
		assert(dm.second.digits.size()<=2);
		if(dm.second.digits.size()>=2)
			val+=((Bigint::longdigit_t)1<<Bigint::digit_bits)*dm.second.digits[1];
		if(dm.second.digits.size()>=1)
			val+=dm.second.digits[0];
		outbuf.push_back(val);
	}
	for(int i=outbuf.size()-1;i>=0;i--){
		(i==(int)outbuf.size()-1?os:os<<setfill('0')<<setw(18))<<outbuf[i];
	}
	return os;
}