#include <cstdlib>
#include <cassert>
#include "numalgo.h"

using namespace std;

int64_t gcd(int64_t a,int64_t b){
	while(true){
		if(a==0)return b;
		if(b==0)return a;
		if(abs(a)>abs(b))a%=b;
		else b%=a;
	}
}

Bigint gcd(Bigint a,Bigint b){
	while(true){
		if(a==0)return b;
		if(b==0)return a;
		if(a.compareAbs(b)==1)a=a.divmod(b).second;
		else b=b.divmod(a).second;
	}
}

Bigint egcd(const Bigint &a,const Bigint &b,Bigint &x,Bigint &y){
	Bigint x2(0),y2(1),r(a),r2(b);
	x=1; y=0;
	while(r2!=0){
		pair<Bigint,Bigint> dm=r.divmod(r2);
		Bigint xn=x-dm.first*x2;
		Bigint yn=y-dm.first*y2;
		x=x2; x2=xn;
		y=y2; y2=yn;
		r=r2; r2=dm.second;
	}
	return r;
}

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--){
		res*=res;
		if(bits[i])res*=b;
		res=res.divmod(m).second;
	}
	return res;
}

int jacobiSymbol(Bigint a,Bigint n){
	//https://en.wikipedia.org/wiki/Jacobi_symbol#Calculating_the_Jacobi_symbol
	assert(n>0&&n.odd());
	int res=1;
	while(true){
		a=a.divmod(n).second;
		if(a<0)a+=n;
		int nmod8=n.divmod(Bigint(8)).second.lowdigits();
		int two_over_n=nmod8==1||nmod8==7?1:-1;
		if(a==0)return 0;
		while(a.even()){
			res*=two_over_n;
			a>>=1;
		}
		if(a==1)return res;
		if(gcd(a,n)!=1)return 0;
		if(a.divmod(Bigint(4)).second.lowdigits()==3&&nmod8%4==3)res*=-1;
		swap(a,n);
	}
}

Bigint isqrt(const Bigint &n){
	assert(n>=0);
	if(n<=1)return n;
	const int maxiter=20; //empirically, this should happen around 3.5 million bits in n. (niter ~= -1.87+1.45ln(bits))
	Bigint x(n);
	x>>=n.bitcount()/2;
	int iter;
	for(iter=0;iter<maxiter;iter++){
		Bigint xnext(x*x);  // xnext = x - (x*x-n)/(2*x)  [Newton's method]
		xnext-=n;
		xnext>>=1;
		xnext=xnext.divmod(x).first;
		xnext.negate();
		xnext+=x;
		if(xnext==x)break;
		x=xnext;
	}
	// cerr<<iter<<endl;
	assert(iter<maxiter);
	switch((x*x).compare(n)){
		case -1:{
			Bigint x1(x);
			x1+=Bigint::one;
			assert(x1*x1>n);
			return x;
		}
		case 0: return x;
		case 1:{
			x-=Bigint::one;
			assert(x*x<=n);
			return x;
		}
		default: assert(false);
	}
}

int ilog2(uint64_t i){
	assert(i);
	int l=0;
	while(i>>=1)l++;
	return l;
}

Bigint bigrandom(Rng &rng,const Bigint &bound){
	const int blocksize=32;
	int btc=bound.bitcount();
	int nblocks=btc/blocksize,rest=btc%blocksize;
	while(true){
		Bigint res;
		for(int i=0;i<nblocks;i++){
			if(i!=0)res<<=blocksize;
			res+=rng.get_uniform((uint32_t)(((uint64_t)1<<blocksize)-1)); //make sure we don't shift out of our int
		}
		if(rest){
			res<<=rest;
			res+=rng.get_uniform((uint32_t)1<<rest);
		}
		if(res<=bound)return res;
	}
}