数论算法模板，最近不更新。
/**-------------------------模意义下运算-------------------------*/
///龟速乘取余，O(logb)
const int mod = 1e9 + 7;
ll qmul(ll a, ll b) {
ll ans = 0;
while (b) {
if (b & 1) ans = (ans + a) % mod;
b >>= 1;
a = (a << 1) % mod;
}
return ans;
}
///快速幂取余，O(logk)，底数越大速度慢的越快
const int mod = 1e9 + 7;
ll qpow(ll a, ll k) {
ll ans = 1;
while (k) {
if (k & 1) ans = (ans * a) % mod;
k >>= 1;
a = (a * a) % mod;
}
return ans;
}
///矩阵快速幂取余，O(n^3logk)
const int mod = 1e9 + 7;
struct Matrix {
int n, m;//不能const，快速幂要复制
vector<vector<ll>> mat;
explicit Matrix(int _n) :n(_n), m(_n), mat(_n + 1, vector<ll>(_n + 1)) {
for (int i = 1;i <= n;i++)
for (int j = 1;j <= m;j++)
mat[i][j] = i == j;
}//初始化n阶单位矩阵，explicit防止误用类型转换
Matrix(int _n, int _m) :n(_n), m(_m), mat(_n + 1, vector<ll>(_m + 1)) {}//初始化nxm零矩阵
friend Matrix operator*(const Matrix &A, const Matrix &B) {
Matrix ans(A.n, B.m);
if (A.m != B.n) return ans;
for (int i = 1;i <= A.n;i++)
for (int j = 1;j <= B.m;j++)
for (int k = 1;k <= A.m;k++) //a.m == b.n
ans.mat[i][j] = (ans.mat[i][j] + A.mat[i][k] * B.mat[k][j]) % mod;
return ans;
}//矩阵乘法取余
friend Matrix operator^(Matrix A, ll k) {
Matrix ans(A.n);//A必须是方阵
while (k) {
if (k & 1) ans = ans * A;
k >>= 1;
A = A * A;
}
return ans;
}//快速幂取余
void output() const {
for (int i = 1; i <= n; i++) {
for (int j = 1; j <= m; j++)
cout << mat[i][j] << ' ';
cout << '\n';
}
cout << '\n';
}//输出检测
};
/**-------------------------素数判断-------------------------*/
///试除法，O(n^1/2)，若是合数，那么在[1,sqrt(n)]中至少有一个质因子
bool isPrime_naive(int n) {
if (n == 2) return 1;
if (n == 1) return 0;
for (int i = 2;i * i <= n;i++) if (!(n % i)) return 0;
return 1;
}
///kn+i法，k=30，O(n^2/15)，kn+i中的i = 1,7,11,13,17,19,23,29时，kn+i的值才可能成为待判数字的质因子
bool isPrime_kni(int n) {
if (n == 2 || n == 3 || n == 5) return 1;
if (n == 1 || !(n % 2) || !(n % 3) || !(n % 5)) return 0;
ll a[8] = { 4,2,4,2,4,6,2,4 }, p = 0;//下一个数字的增量
for (int i = 7;i * i <= n;i += a[p++], p %= 8) if (!(n % i)) return 0;
return 1;
}
///预处理法，对于多组数据，如果n是合数就一定有在[1,n^1/2]的质因子，先预处理出[1,n^1/2]的所有素数,然后对n测试
/**-------------------------素数筛-------------------------*/
///埃氏筛，O(nloglogn)，通过标记素数的倍数筛掉合数
const int N = 1e6 + 7;
bool vis[N];
int prime[N], cnt;
void eratosthenes_screen(int n) {
for (int i = 2;i <= n;i++) {
if (vis[i]) continue;
prime[cnt++] = i;
for (int j = 2;j * i <= n;j++) vis[i * j] = 1;
}
}
///欧拉筛(线性筛)，O(n)，每个合数只会被最小质因子筛掉
const int N = 1e6 + 7;
bool vis[N];
int prime[N], cnt;
void euler_screen(int n) {
for (int i = 2;i <= n;i++) {
if (!vis[i]) prime[cnt++] = i;
for (int j = 0;j < cnt && i * prime[j] <= n;j++) {
vis[i * prime[j]] = 1;
if (!(i % prime[j])) break;//如果到了i的最小质因子就不用继续，因为接下去的数x一定能被(i,x)之间的数筛掉
}
}
}
/**-------------------------反素数-------------------------*/
/**-------------------------唯一分解定理-------------------------*/
/**-------------------------正整数集的定理-------------------------*/
/**-------------------------剩余类集结构-------------------------*/
/**-------------------------约数-------------------------*/
/**-------------------------最大公约数-------------------------*/
///辗转相除法(欧几里得算法)，O（logn)，大整数取模耗时大。b!=0，gcd(n,m) = gcd(b,a mod b)。
int gcd(int a, int b) {
return b ? gcd(b, a % b) : a;
}
///更相减损术，gcd(n,m) = gcd(n,n-m)。
///Stein算法，O(logn),避免大整数取模的耗时。更相减损术的应用，gcd(ka,kb) = k*gcd(a,b)。
//递归版
int stein_recursion(int a, int b) {
if (a < b) a ^= b, b ^= a, a ^= b;//交换使a为较大数
if (!b) return a;//两数相等时，gcd=a
if (!(a & 1) && !(b & 1)) return stein_recursion(a >> 1, b >> 1) << 1;//最后要把因子2乘上
else if (a & 1 && !(b & 1)) return stein_recursion(a, b >> 1);
else if (!(a & 1) && b & 1) return stein_recursion(a >> 1, b);
else return stein_recursion(a - b, b);
}
//迭代版
int stein_iteration(int a, int b) {
int k = 1;
while (!(a & 1) && !(b & 1)) k <<= 1, a >>= 1, b >>= 1;//k记录因子2乘积
while (!(a & 1)) a >>= 1;
while (!(b & 1)) b >>= 1;
if (a < b) a ^= b, b ^= a, a ^= b;//使a较大
while (a != b) {
a -= b;
if (a < b) a ^= b, b ^= a, a ^= b;
}
return k * a;
}
/**-------------------------最小公倍数-------------------------*/
///最大公倍数，利用了gcd(a,b)*lcm(a,b)=a*b的性质求解
int lcm(int a, int b) {
return a / gcd(a, b) * b;//先除后乘避免溢出
}
/**-------------------------欧拉函数-------------------------*/
///分解质因数求欧拉函数，O(n^1/2),利用公式phi[n] = n*(1-1/n的素因子)连乘。搭配线性筛直接枚举素数能达到O((n^1/2)/logn)
int euler_one(int n) {
int ans = n;
for (int i = 2;i * i <= n;i++) {
if (!(n % i)) {
ans = ans / i * (i - 1);
while (!(n % i)) n /= i;
}
}
if (n > 1)ans = ans / n * (n - 1);//大于n^1/2的素因子至多只有一个且次数为1，如此特判(包括n本身就是素数)即可
return ans;
}
/**-------------------------扩展欧几里得-------------------------*/
///扩展欧几里得算法，在gcd基础上增加了扩展，可以得到ax+by = gcd(a,b)的特解。
int exgcd(int a, int b, int &x, int &y) {
if (!b) {
x = 1, y = 0;
return a;
}
int d = exgcd(b, a % b, x, y);
int t = x;
x = y;
y = t - (a / b) * y;
return d;
}
/**-------------------------乘法逆元-------------------------*/
///扩展欧几里得算法求逆元，b，m互质，解丢番图方程ax+my=1，得到逆元
int exgcd_inverse(int a, int m) {
int x, y;
exgcd(a, m, x, y);
return (m + x % m) % m;
}
///线性递推求乘法逆元
///求阶乘逆元
/**-------------------------省略-------------------------*/
/**-------------------------线性筛法-------------------------*/
///线性筛求欧拉函数，O(n)，利用筛数时的最小质因子来递推欧拉函数
const int N = 1e6 + 7;
bool vis[N];
int prime[N], cnt, phi[N];
void get_euler(int n) {
phi[1] = 1;//莫忘
for (int i = 2;i <= n;i++) {
if (!vis[i]) {
prime[cnt++] = i;
phi[i] = i - 1;//素数
}
for (int j = 0;j < cnt && i * prime[j] <= n;j++) {
vis[i * prime[j]] = 1;
if (!(i % prime[j])) {//prime[j]是i最小质因子，则i*prime[j]有至少两个最小质因子
phi[i * prime[j]] = phi[i] * prime[j];
break;
}
phi[i * prime[j]] = phi[i] * (prime[j] - 1);//prime[j]非i最小质因子，则自身就是唯一最小质因子
}
}
}
///线性筛求莫比乌斯函数
/**-------------------------杜教筛-------------------------*/
