数学模板 - 数论基础

int PowerMod(int a, int b, int k) ///a ^ b % k
{
int tmp = a, ret = 1;
while(b)
{
if(b & 1)
ret = ret * tmp % k;
tmp = tmp * tmp % k;
b >>= 1;
}
return ret;
}

int gcd(int a, int b) ///最大公约数
{
return b == 0 ? a : gcd(b, a % b);
}
int lcm(int a, int b) ///最小公倍数
{
if (a * b == 0) return 0;
return a / gcd(a, b) * b;
}

int Ext_gcd(int a,int b,int &x,int &y)///求解方程 ax + by = gcd(a,b) 的整数解(其中a,b都为非零整数)
{
int d;
if(b == 0)
{
x = 1;  y = 0;
return a;
}
d = Ext_gcd(b, a%b, y, x);
y -= a / b * x;
return d;///返回gcd(a,b)
}

int inv(int a, int n)///a在模n乘法下的逆元,没有则返回-1(a * b % n == 1,已知a,n,求b就是乘法逆元)
{
int x, y;
int t = Ext_gcd(a, n, x, y);///扩展欧几里得算法
if(t != 1)
return -1;
return (x % n + n) % n;
}

int Chinese_Remainder(int a[],int w[],int len)///中国剩余定理,a[]存放余数,w[]存放两两互质的数
{
int d,x,y,m,n,ret;
ret = 0;
n = 1;
for(int i = 0; i < len; i++)
n *= w[i];
for(int i = 0; i < len; i++)
{
m = n / w[i];
d = Ext_gcd(w[i],m,x,y);///扩展欧几里得算法
ret = (ret + y * m * a[i]) % n;
}
return (n + ret % n) % n;
}

bool merge(int a1, int n1, int a2, int n2, int& a3, int& n3)///将两个方程合并为一个
{
int d = gcd(n1, n2);
int c = a2 - a1;
if(c % d)
return false;
c = (c % n2 + n2) % n2;
c /= d;
n1 /= d;
n2 /= d;
c *= inv(n1,n2);///乘法逆元
c %= n2;
c *= n1 * d;
c += a1;
n3 = n1 * n2 * d;
a3 = (c % n3 + n3) % n3;
return true;
}

int China_Reminder2(int len, int* a, int* n)///求模线性方程组x=ai(mod ni),ni可以不互质
{
int a1 = a[0], n1 = n[0];
int a2, n2;
for(int i = 1; i < len; i++)
{
int aa, nn;
a2 = a[i], n2 = n[i];
if(!merge(a1,n1,a2,n2,aa,nn))
return -1;///不存在在结果
a1 = aa;
n1 = nn;
}
Mod = n1;
return (a1 % n1 + n1) % n1;
}

void Euler()
{
int i, j;
for (i = 1; i <= Max; i++)
num[i] = i;
for (i = 2; i <= Max; i += 2)
num[i] /= 2;
for (i = 3; i <= Max; i += 2)
{
if(num[i] == i)
{
for (j = i; j <= Max; j += i)
num[j] = num[j] / i * (i - 1);
}
}
}

int Euler(int x)
{
int res = x;
for(int i = 2; i < (int)sqrt(x * 1.0) + 1; i++)
{
if(x % i == 0)
{
res = res / i * (i - 1);
while (x % i == 0)/// 保证i一定是素数
x /= i;
}
}
if(x > 1)
res = res / x * (x - 1);
return res;
}

void isprime(int Max)
{
int i, j;
memset(num, 0, sizeof(num)); /// 初始赋值所有数都是素数.
num[1] = 1;
num[0] = 1; ///0和1不是素数
for (i = 2; i < Max; i++)
{
if (!num[i])
{
if (1LL * i * i > 1LL * Max)
continue;
for(j = i * i; j < Max; j += i)
num[j] = 1;
}
}
}

bool vis[Max]; /// 判断是否为素数
int PrimeList[Max]; /// 素数列表
void Prime_Linear(void)
{
int num = 0;
memset(vis, true, sizeof(vis));
vis[1] = 0;
vis[0] = 0;
for (int i = 4; i < Max; i += 2)
vis[i] = 0;
for (i = 3; i < Max; i += 2)
{
if (vis[i])///如果是素数则加入素数列表.
PrimeList[num++] = i;
for (j = 0; j < num; j++)
{
if (i * PrimeList[j] > Max)
break;
vis[i * PrimeList[j]] = 0;
if (i % PrimeList[j] == 0)
break;
}
}
}

int PrimeList[maxn];
bool vis[maxn]; /// vis[i] = 1表示i + L这个数是素数.
int num;
void SegmentPrime(int L, int U)///求区间[L, U]中的素数.
{
int i, j;
int SU = sqrt(1.0 * U);
int d = U - L + 1;
for (i = 0; i < d; i++)/// 初始全是素数
vis[i] = 1;
for (i = (L % 2 != 0); i < d; i += 2)
vis[i] = 0;
for (i = 3; i <= SU; i += 2)
{
if (i > L && !vis[i - L])
continue;
j = (L / i) * i;
if (j < L)
j += i;
if (j == i)
j += i;
j = j - L;
for (; j < d; j += i)
vis[j] = 0;
}
if (L <= 1)
vis[1 - L] = 0;
if (L <= 2)
vis[2 - L] = 1;
num = 0;
for (i = 0; i < d; i++)
if (vis[i])
PrimeList[num++] = i + L;
}

#include<stdlib.h>
bool Witness(ll a,ll n)
{
ll t, d = 1, x;
int i = ceil(log(n - 1.0) / log(2.0)) - 1;
for(; i >= 0; i--)
{
x = d;
d = (d * d) % n;
if(d == 1 && x != 1 && x != n-1)
return true;
if( ((n-1) & (1<<i)) > 0)
d = (d * a) % n;
}
return d == 1 ? false : true;
}

bool Miller_Rabin(ll n)///随机素数判断
{
if(n == 2)
return true;
if(n == 1 || ((n & 1) == 0))
return false;
for(int i = 0; i < 50; i++)
{
ll a = rand() * (n-2) / RAND_MAX + 1;
if(Witness(a, n))
return false;
}
return true;
}

void PrimeFactor(void)///factor存储最小质因数
{
int i, j;
memset(factor, 0, sizeof(factor));
factor[1] = 1;
factor[0] = 1;
int PrimeNum = 0;
for (i = 2; i < Max; i++)
{
if (!factor[i])
PrimeList[PrimeNum++] = i;
for (j = 0; j < PrimeNum && i * PrimeList[j] < Max && (PrimeList[j] <= factor[i] || factor[i] == 0); j++)
{
factor[i * PrimeList[j]] = PrimeList[j];
if (i % PrimeList[j] == 0)
break;
}
}
/*
for(int i = 1; i <= Max; i++)
printf("%d %d\n",i, factor[i]);
*/
}