/**************以下为常用算法************/
//把n的约数的莫比乌斯函数值用map的形式返回O(sqrt(n))
map<int,int > moebius (int n)
{
map<int,int > res;
vector<int > pri;
for (int i=2;i*i<=n;i++)
if (n%i==0){
pri.push_back(i);
while (n%i==0) n/=i;
}
if (n!=1) pri.push_back(n);
int mm=pri.size();
for (int i=0;i<(1<<mm);i++)
{
int mu=1,d=1;
for (int j=0;j<mm;j++)
if (i>>j&1){
mu*=-1;
d*=pri[j];
}
res[d]=mu;
}
return res;
}
//-------------------------------
int euler(int n)
{
int res=n;
for (int i=2;i*i<=n;i++)
if (n%i==0)
{
res=res/i*(i-1);
for (;n%i==0;n/=i);
}
if (n!=1) res=res/n*(n-1);
return res;
}
//-------------------------------
bool v[N];//筛法的辅助数组
int num[N],o,e[N],mu[N];//num存素数,e存欧拉函数,mu存莫比乌斯函数
void get_euler_and_prime()
{
memset(v,0,sizeof(v));e[1]=1; o=0;
for (int i=2;i<N;i++)
{ if (!v[i]) {num[o++]=i; e[i]=i-1; mu[i]=-1;}
for (int j=0;j<o;j++)
{
if (i*num[j]>=N) break;
v[i*num[j]]=1;
if (i%num[j]==0) {e[i*num[j]]=e[i]*num[j];break;
mu[i*num[j]]=0;
}
else {e[i*num[j]]=e[i]*e[num[j]];
mu[i*num[j]]=-mu[i];
}
}
}
}
//-------------------------------
ll mul_(ll a,ll b,ll mod)
{
ll exp=a%mod,res=0;
while (b)
{
if (b&1){ res+=exp;
if (res>=mod) res-=mod;
}
exp<<=1;
if (exp>=mod) exp-=mod;
b>>=1;
}
return res;
}
//-------------------------------
ll pow_(ll a,ll n,ll mod)
{ ll res=1;
while (n)
{ if (n&1) res=mul_(res,a,mod);
n>>=1;
a=mul_(a,a,mod);
}
return res;
}
//-------------------------------
ll pow2_(ll a,ll n,ll mod)
{
if (n==0) return 1;
if (n%2==1) return mul_(pow2_(a,n/2,mod),pow_(a,n/2+1,mod)+1,mod);
return (pow2_(a,n-1,mod)+pow_(a,n,mod)+1)%mod;
}
//-------------------------------
void exgcd(ll a,ll b,ll &d,ll &x,ll&y)//ax+by=d;且返回的x,y满足:abs(x)+abs(y)最小
{
if (!b) {d=a;x=1;y=0;}
else {exgcd(b,a%b,d,y,x); y-=a/b*x;}
}
//-------------------------------
ll minx_(ll a,ll b,ll c)//ax≡c(mod b),求最小非负x
{
ll d,x,y;
exgcd(a,b,d,x,y);
if (c%d!=0) return -1;
a/=d;
b/=d;
c/=d;
return (x*c%b+b)%b;
}
//-------------------------------
ll inv(ll a,ll n)
{
ll d,x,y;
exgcd(a,n,d,x,y);
return d==1?(x+n)%n:-1;
}
//-------------------------------
ll fac[P]; 预处理的n!%p的表 O(p)
ll mod_fact(ll n,ll p,ll &e) //要求p为素数,O(log2(n)/log2(p))
{
e=0;if (n==0) return 1;
ll res=mod_fact(n/p,p,e);
e+=n/p;
if (n/p%2!=0) return res*(p-fac[n%p])%p;
return res*fac[n%p]%p;
}
//-------------------------------
ll mod_comb(ll n,ll k,ll p) //要求p为素数,O(log2(n)/log2(p))
{
if (n<0||k<0||n<k) return 0;
ll e1,e2,e3,a1,a2,a3;
a1=mod_fact(n,p,e1);
a2=mod_fact(k,p,e2);
a3=mod_fact(n-k,p,e3);
if (e1>e2+e3) return 0;
return a1*inv(a2*a3%p,p)%p;
}
//-------------------------------
map<ll,ll> ma;
ll log_mod(ll a,ll b,ll n)//a^x≡b(mod n)
{
ll m=1,v,e=1,i;
m=ll(sqrt(n+0.5));
v=inv(pow_(a,m,n),n);
ma.clear();ma[1]=0;
for (int i=1;i<m;i++)
{
e=mul_(e,a,n);
if (!ma.count(e)) ma[e]=i;
}
for (int i=0;i<m;i++)
{
if (ma.count(b)) return i*m+ma[b];
b=mul_(b,v,n);
}
return -1;
}
/***************以下为解方程组的******************/
pair<ll,ll> china(int n,int *a,int *m)
{
ll M=1,d,y,x,re=0;
for (int i=0;i<n;i++) M*=m[i];
for (int i=0;i<n;i++)
{
ll w=M/m[i];
exgcd(m[i],w,d,x,y);
re=(re+y*w*a[i])%M;
}
return make_pair((re+M)%M,M);
}
//-------------------------------
pair<ll,ll> equs(int *aa,int *bb,int *cc,int nn)
{
ll re=0,m=1,t,a,b,d;
for (int i=0;i<nn;i++)
{
a=aa[i]*m,b=bb[i]-aa[i]*re,d=gcd(cc[i],a);
if (b%d!=0) return make_pair(ll(0),ll(-1));
t=b/d*inv(a/d,cc[i]/d)%(cc[i]/d);
re=re+m*t;
m*=cc[i]/d;
}
return make_pair((re%m+m)%m,m);
}
/**************以下为两个随机算法************/
//需要主程序开头加srand(time(NULL));以及ctime头文件
bool check(ll a,ll n,ll r,ll s)
{
ll ans,last,i;
ans=pow_(a,r,n);
last=ans;
for(i=1;i<=s;i++)
{
ans=mul_(ans,ans,n);
if(ans==1&&last!=1&&last!=n-1)return true;
last=ans;
}
if(ans!=1)return true;
return false;
}
bool Miller_Rabin(ll n)//Miller_Rabin算法,判断n是否为素数
{
if(n<2)return false;
if(n==2)return true;
if(!(n&1))return false;
ll r,s,a;
r=n-1;s=0;
while(!(r&1)){r=r>>1;s++;}
for(int i=0;i<20;i++)
{
a=rand()%(n-1)+1;
if(check(a,n,r,s))
return false;
}
return true;
}
//-------------------------------
ll gcd(ll a,ll b){return b==0?a:gcd(b,a%b);}
ll Pollard_rho(ll n,ll c)
{
ll i=1,j,k=2,x,y,d,p;
x=rand()%n;
y=x;
while(true)
{
i++;
x=(mul_(x,x,n)+c)%n;
if(y==x)return n;
if(y>x)p=y-x;else p=x-y;
d=gcd(p,n);
if(d!=1&&d!=n)return d;
if(i==k)
{
y=x;
k+=k;
}
}
}
void find(ll n)//找出n的所有质因子
{
if(Miller_Rabin(n))
{
Q.push(n);return;
}
ll p=n;
while(p>=n)p=Pollard_rho(p,rand()%(n-1)+1);
find(p);
find(n/p);//如果只要求找到一个质因子,可以注释掉
}