数论的模板

//欧拉筛

#include<iostream>
#include<algorithm>
#include<vector>
#include<stdlib.h>
#define ll long long 
using namespace std;

const int maxn = 1000010;
bool visited[maxn];
vector<int>prime;

//欧拉筛模板
void isprime()
{
	for (int i = 2; i < maxn; i++) {
		if (visited[i] == false) {
			prime.push_back(i);
			visited[i] = true;
		}
		for (int j = 0; j < prime.size()&&i*prime[j]<maxn; j++) {
			visited[i*prime[j]] = true;
			if (i%prime[j] == 0)break;
		}
	}
}



//gcd模板,非递归
ll gcd(ll a, ll b)
{
	if (a == 0)return 1;
	if (a < 0)gcd(-a, b);
	while (b)
	{
		ll k = b;
		b = a%b;
		a = k;
	}
	return a;
}
//ex_gcd模板
//ax+by=gcd(a,b)
ll ex_gcd(ll a, ll b,ll &x, ll &y)
{
	if (b == 0) {
		x = 1;
		y = 0;
		return a;
	}
	ll temp, d;
	d = ex_gcd(b, a%b, x, y);
	temp = x;
	x = y;
	y = temp - (a / b)*y;
	return d;
}
//用扩欧得到一组解x0,y0之后，a*x0+b*y0=gcd(a,b),两边同除gcd(a,b)再乘上n，得到方程的解
//若gcd(a,b)=1,则通解为x=x0+b*t,y=y0-a*t
//求解最小正整数解时，t=b/gdc(a,b),x=(x%t+t)%t


//求解逆元
ll inv(ll a, ll m)
{
	ll d, x, y;
	d = ex_gcd(a, m, x, y);
	if (d == 1)return (x%m + m) % m;
	else return -1;
}

//分解质因数模板
vector<int>P;
void div(ll n)
{
	for (int i = 2; i*i <= n; i++) {
		if (n%i == 0) {
			P.push_back(i);
			while (n%i == 0)n /= i;
		}
	}
	if (n > 1)P.push_back(n);
}

//中国剩余定理模板
ll china_mod(ll mod[], ll a[],int n)
{
	ll lcm=1, i, ans, Mi, x, y;
	for (int i = 0; i < n; i++)
		lcm *= mod[i];
	ans = 0;
	for (int i = 0; i < n; i++) {
		Mi = lcm / mod[i];
		ex_gcd(Mi, mod[i], x, y);
		ans += Mi*x*a[i]%lcm;
	}
	ans %= lcm;
	while (ans>0)
	{
		ans += lcm;
	}
	return ans;
}


//ll a*ll b%c的算法
ll mul_mod(ll a, ll b, ll c)
{
	a %= c;
	b %= c;
	ll ret = 0;
	while (b)
	{
		if (b & 1) {
			ret += a;
			ret %= c;
		}
		a <<= 1;
		if (a >= c)a %= c;
		b >>= 1;
	}
	return ret;
}

//快速幂模板
ll pow_mod(ll a, ll n, ll m)
{
	ll res = 1;
	while (n)
	{
		if (n & 1)res =res*a%m;
		a = a*a%m;
		n >>=1;
	}
	return res;
}
//快速幂模板2
ll pow_mod2(ll x, ll n, ll mod)
{
	if (n == 1)return x%mod;
	x %= mod;
	ll temp = x;
	ll ret = 1;
	while (n)
	{
		if (n & 1)ret = mul_mod(ret, temp, mod);
		temp = mul_mod(temp, temp, mod);
		n >>= 1;
	}
	return ret;
}

//欧拉函数筛
int phi[maxn];
vector<int>is_prime;
void init()
{
	phi[1] = 1;
	for (int i = 2; i < maxn; i++) {
		if (phi[i] == 0) {
			is_prime.push_back(i);
			phi[i] = i - 1;
		}
		for (int j = 0; j < is_prime.size(); j++) {
			if (i%is_prime[j] == 0) {
				phi[i*is_prime[j]] = phi[i] * prime[j];
				break;
			}
			else phi[i*prime[j]] = phi[i] * (prime[j] - 1);
		}
	}
}


//求一个数的欧拉函数
ll ol(ll x)
{
	ll i, res = x;
	for (int i = 2; i*i <= x; i++) {
		if (x%i == 0) {
			res = res - res / i;
			while (x%i == 0)x /= i;
		}
	}
	if (x > 1)res - res / x;
	return res;
}


//二次探测原理判断是否是素数
// Miller_Rabin()算法素数判定
//是素数返回true.(可能是伪素数，但概率极小)
//合数返回false;
bool check(ll a, ll n, ll d, ll r)
{
	ll x = pow_mod(a, d, n);
	ll last = x;
	for (int i = 1; i <= r; i++) {
		x = x*x%n;      //为防止溢出可以使用x=mul_mod(x,x,n);
		if (x == 1 && last != 1 && last != n - 1)
			return true;      //合数
		last = x;
	}
	if (x != 1)return true;
	return false;
}

bool Miller_Rabin(ll n)
{
	if (n == 1)return false;
	if (n == 2)return true;
	if ((n & 1) == 0)return false;
	ll d = n - 1;
	ll r = 0;
	while ((d & 1) == 0) {
		d >>= 1; r++;
	}
	srand(time_t(0));
	for(int i=0;i<50;i++)//测试次数
	{
		ll a = rand() % (n - 1) + 1;
		if (check(a, n, d, r))
			return false;
	}
	return true;
}


//************************************************
//pollard_rho 算法进行质因数分解
//************************************************
long long factor[100];//质因数分解结果（刚返回时是无序的）
int tol;//质因数的个数。数组小标从0开始


long long Pollard_rho(long long x, long long c)
{
	long long i = 1, k = 2;
	long long x0 = rand() % x;
	long long y = x0;
	while (1)
	{
		i++;
		x0 = (mul_mod(x0, x0, x) + c) % x;
		long long d = gcd(y - x0, x);
		if (d != 1 && d != x) return d;
		if (y == x0) return x;
		if (i == k) { y = x0; k += k; }
	}
}
//对n进行素因子分解
void findfac(long long n)
{
	if (Miller_Rabin(n))//素数
	{
		factor[tol++] = n;
		return;
	}
	long long p = n;
	while (p >= n)p = Pollard_rho(p, rand() % (n - 1) + 1);
	findfac(p);
	findfac(n / p);
}


//二次剩余的模板：x*x%p=a,求解x,称x在mod p下的二次剩余,若有解，a称为mod p下的完全平方数
//若a^2-n不是Mod p下的完全平方数，则x=(a+aqrt(a^2-n))^(p+1)/2为mod p 的二次剩余
struct T
{
	ll x, y;
	T() {};
	T(ll _x, ll _y)
	{
		x = _x, y = _y;
	}
};
ll w;
T mul(T a, T b, ll p)
{
	T ans;
	ans.x = (a.x*b.x%p + a.y*b.y%p*w%p) % p;
	ans.y = (a.x*b.y%p + a.y*b.x%p) % p;
	return ans;
}
T T_mod_pow(T a, ll b, ll p)
{
	T ans = T(1, 0);
	while (b)
	{
		if (b & 1)ans = mul(ans, a, p);
		a = mul(a, a, p);
		b >>= 1;
	}
	return ans;
}

ll Solve(ll a, ll p)//求解x^2=a(mod p) 
{
	a %= p;
	if (p == 2)return a;
	if (pow_mod(a, (p - 1) / 2, p) == p - 1)return -1;
	int b, t;
	while (1)
	{
		b = rand() % p;
		t = b*b - a;
		w = (t%p + p) % p;
		if (pow_mod(w, (p - 1) / 2, p) == p - 1)break;
	}
	T ans = T(b, 1);
	ans = T_mod_pow(ans, (p + 1) / 2, p);
	return ans.x;
}