[专题-数论]数论基础

数论大礼包：
[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，