#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(const string &repr){
	stringstream(repr)>>*this;
}

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

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

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){
		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){
		//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();
}

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=(slongdigit_t v){
	digits.clear();
	if(v==0){
		sign=1;
		checkconsistent();
		return *this;
	}
	sign=v>=0?1:-1;
	v*=sign;
	digits[0]=v;
	if(v>digits[0])digits.push_back(v>>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
		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
		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;
	if(sh/digit_bits>0){
		digits.insert(digits.begin(),sh/digit_bits,0);
		sh%=digit_bits;
		if(sh==0){
			checkconsistent();
			return *this;
		}
	}
	digits.push_back(0);
	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;
	if(sh/digit_bits>0){
		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;
	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 {
	// cerr<<"divmod("<<hex<<*this<<','<<div<<dec<<')'<<endl;
	int bitcdiff=bitcount()-div.bitcount();
	if(bitcdiff<0)bitcdiff=0;
	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;
	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){ //the large case
		//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);
	}
}

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

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))&255;
		}
	}
	return s;
}

void Bigint::deserialiseMantissa(const string &s){
	assert(s.size()%4==0);
	sign=1;
	int sz=s.size()/4;
	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]|=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;
}

Bigint::digit_t Bigint::_digit(int idx) const {
	return digits.at(idx);
}

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'){
			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;
}