数论大礼包：
[NWPU][2018暑假集训]day10&11 素数+GCD
[NWPU][2018暑假集训]day10&11 逆元+快速幂
[NWPU][2018暑假集训]day10&11 同余+欧拉函数
[NWPU][2018暑假集训]day10&11 勾股数+佩尔方程
模版
求乘法逆元
扩展欧几里德求逆元 O(logn)
原理即ax+by=c两边膜b，得ax=c(mod b)，利用扩欧求解即可
int exgcd(int a, int b, int &x, int &y){
if (b==0){x=1; y=0; return a;}
int gcd=exgcd(b, a%b, y, x);
y-=(a/b)*x;
return gcd;
}
int inv(int a, int p){
int x, y, gcd=exgcd(a, p, x, y);
if (gcd==1) return (x%p+p)%p;
return -1;
}
求逆元表 O(n)
对于不能求逆元的情况(A/B)%mod = (A%(B*mod))/B%mod
long long inv[mod];
void init(void){
inv[1]=1;
for (int i=2; i<mod; i++)
inv[i]=(long long)(mod-mod/i)*inv[mod%i]%mod;
}
筛法求素数
埃氏筛 O(nlog(logn))
const int maxn=1e7, maxp=7e5;
int primes[maxp+5], psize;
bool isprime[maxn+5];
void initPrime(void){
memset(isprime, true, sizeof(isprime));
for (int i=2; i<=maxn; i++) if (isprime[i]){
for (int j=i*2; j<=maxn; j+=i)
isprime[j]=false;
primes[psize++]=i;
}
}
欧拉筛 O(n)
const int maxn=1e5+20;
int primes[maxn/10], psize;
bool isprime[maxn];
void initPrimes(void){
memset(isprime, true, sizeof(isprime));
isprime[0]=isprime[1]=false;
for (int i=2; i<=maxn; i++){
if(isprime[i]) primes[psize++]=i;
for (int j=0; j<psize && i*primes[j]<=maxn; j++){
isprime[primes[j]*i]=false;
if (i%primes[j]==0) break;
}
}
}
质因数分解（唯一分解定理）
利用素数表 O(sqrt(n)/logn) + O(n)
注意maxn可为需分解的最大数的开方
const int maxn=1e5+20;
int factors[100][2], fsize, primes[maxn/10], psize;
void getFactors(long long n){
fsize=0;
for (int i=0; i<psize && primes[i]<=n/primes[i]; i++){
if (n%primes[i]==0){
factors[fsize][0]=primes[i];
factors[fsize][1]=0;
while (n%primes[i]==0) factors[fsize][1]++, n/=primes[i];
fsize++;
}
}
if (n>1){
factors[fsize][0]=n;
factors[fsize++][1]=1;
}
}
不用素数表 O(sqrt(n))
int factors[100][2], fsize;
void getFactors(long long n){
fsize=0;
for (int i=2; i<=n/i; i++){
if (n%i==0){
factors[fsize][0]=i;
factors[fsize][1]=0;
while (n%i==0) factors[fsize][1]++, n/=i
fsize++;
}
}
if (n>1){
factors[fsize][0]=n;
factors[fsize++][1]=1;
}
}
最大公因数&最小公倍数
// gcd(a, b, c)==gcd(gcd(a, b), c)
// lcm(a, b, c)==lcm(lcm(a, b), c)
long long gcd(long long a, long long b){
return (b==0)?a:gcd(b, a%b);
}
long long lcm(long long a, long long b){
return a/gcd(a, b)*b;
}
long long exgcd(long long a, long long b, long long &x, long long &y){
if (b==0){x=1; y=0; return a;}
long long gcd=exgcd(b, a%b, y, x);
y-=(a/b)*x;
return gcd;
}
分块打表
例：求1~n的调和级数(LightOJ-1234)
const int maxn=1e6;
double h[maxn+5];
void init(void){
h[0]=h[1]=0;
for (int i=1, ptr=1; i<=1e8; i++){
if (i%100==0) {ptr++; h[ptr]=h[ptr-1];}
h[ptr]+=1/(double)i;
}
}
double calc(int n){
double sum=0;
for (int i=(n/100)*100; i<=n; i++)
if (i!=0) sum+=1/(double)i;
return sum+h[n/100];
}
快速幂
这里贴快速幂的原因是某些数太大，直接快速幂容易long long溢出
long long quickMult(long long a, long long n, long long mod){
long long ans=0, tmp=a;
for (int i=0; (1<<i)<=n; i++){
if ((1<<i)&n) ans=(ans+tmp)%mod;
tmp=(long long)(tmp+tmp)%mod;
}return ans;
}
long long quickPow(long long a, long long n, long long mod){
int ans=1, tmp=a;
for (int i=0; (1<<i)<=n; i++){
if ((1<<i)&n) ans=quickMult(ans, tmp)%mod;
tmp=quickMult(tmp, tmp)%mod;
}return ans;
}
矩阵快速幂
const int maxn=20;
const long long mod=1000000007;
struct Matrix{
int r, c;
long long mat[maxn][maxn];
Matrix(int r, int c):r(r), c(c) {}
void clear(void){memset(mat, 0, sizeof(mat));}
};
Matrix operator + (Matrix a, Matrix b){
Matrix s(a.r, a.c);
for(int i = 0; i < a.r; i++)
for(int j = 0; j < a.c; j++)
s.mat[i][j]=(a.mat[i][j]+b.mat[i][j])%mod;
return s;
}
Matrix operator * (Matrix a, Matrix b){
Matrix s(a.r, b.c);
for(int i = 0; i < a.r; i++)
for(int k = 0; k < a.c; k++)
for(int j = 0; j < b.c; j++)
s.mat[i][j]=(s.mat[i][j]+a.mat[i][k]*b.mat[k][j])%mod;
return s;
}
Matrix pow(Matrix a, long long n){
Matrix ret(a.r, a.c), tmp(a);
for(int i = 0; i < a.r; i++)
ret.mat[i][i]=1;
while(n){
if(n&1) ret=ret*tmp;
tmp=tmp*tmp;
n>>=1;
}return ret;
}
同余方程
解一元线性同余方程
// 求ax=b (mod m)最小解
long long solve(long long a, long long b, long long m){
long long x, y, gcd=exgcd(a, m, x, y);
if(b % d == 0) {
x = x * (b /gcd);
x = (x%(m/gcd) + (m/gcd)) % (m/gcd);
return x;
}
return -1;
}
迭代法解同余方程组
为什么不用中国剩余定理？迭代法不要求互素
long long a[maxn], m[maxn];
bool solve(long long &m0, long long &a0, int n){
m0 = 1; a0 = 0;
for(int i = 0; i < n; i++) {
long long t, s, t0;
long long d = exgcd(m0, m[i], t, s);
if((a[i] - a0) % d != 0) return false;
t *= (a[i] - a0) / d;
t0 = (t % (m[i] / d) + (m[i] / d)) % (m[i] / d);
a0 += m0 * t0;
m0 *= (m[i] / d);
a0 %= m0;
}
return true;
}
欧拉函数
费马小定理
p为素数，则a^p=a (mod p)
若a, p互素，则a^(p-1)=1 (mod p)
欧拉定理
任意正整数n，a^(phi(n)+1)=a (mod n)
若a, n互素，a^phi(n)=1 (mod n)
欧拉函数、莫比乌斯函数是积性函数( f(nm)=f(n)*f(m) )
应用
求逆元：a,n互素时，inv(a, p)==a^(phi(n)-1)
降幂：当b>phi(n)时（即指数超大），a^b=a^( b%phi(n)+phi(n) ) (mod n)注意a, n不需要互素
分解质因数方法 O(sqrt(n))
long long eular(long long n) {
getFactors(n);
long long ret = n;
for(int i=0; i<cnt; i++)
ret=ret/factor[i][0]*(factor[i][0]-1);
return ret;
}
筛法
const int maxn=1e6;
int phi[maxn+5];
long long sum[maxn+5];
void initPhi(void){
memset(phi, 0, sizeof(phi));
phi[1]=1;
for (int i=2; i<=maxn; i++) if (!phi[i])
for (int j=i; j<=maxn; j+=i){
if (!phi[j]) phi[j]=j;
phi[j]=phi[j]/i*(i-1);
}
}
非线性丢番图方程
毕达哥拉斯三元组
x^2+y^2=z^2
当gcd(x, y, z)=1时，称此三元组是本原的
存在互素且奇偶性不同的正整数n, m
x=m^2-n^2, y=2mn, z=m^2+n^2
费马大定理
x^n+y^n=z^n在n>2时无非零整数解
佩尔方程
x^2-Dy^2=1
最小整数解，D<30时暴力枚举
递推法求通解
对于 x2 - Dy2 = M，其中 M = ±1, ±2, ±4，若方程存在基本解 (x1, y1)，则有
xn = C(xn-1 - xn-2) + xn-3，
yn = C(yn-1 - yn-2) + yn-3，
当 M 为 1, 2, 4, -1, -2, -4 时，
C 分别为 2x1+1, 2x12-1, x1+1, 4x12+3, 2x12+3, x12+3，