数论代码补充

1.素数的检验（试除法）

bool is_prime(int x)
{
    for(int i = 2; i<=sqrt(n); i++)
    {
        if(x%i == 0)
            return false;
    }
    return true;
}

2.素数的筛选

Eratothenes筛法

int vis[maxn];
void primes(int n)
{
    memset(vis, 0, sizeof vis);
    for(int i = 2; i<=n; i++)
    {
        if(vis[i])
            continue;
        cout<<"i是质数"<<endl;
        for(int j = i; j<=n/i; j++)
            vis[i*j] = 1;
    }
}

线性筛法

Eratothenes筛法会重复标记合数，如12就会又被2扫到又被3扫到。

线性筛法通过“从小到大产生质因子”的方式标记合数，设数组v记录每个数的最小质因子。

1.一次考虑2~N之间的每一个数i

2.若v[i] = i，说明i是质数

3.扫描不大于v[i]的每个质数p，令v[i*p] = p

每个合数i*p只会被它的最小质因子p晒到一次。

int v[maxn],prime[maxn];
void primes(int n)
{
    memset(v, 0, sizeof v);
    int cnt = 0;
    for(int i = 2; i<=n; i++)
    {
        if(v[i] == 0)
        {
            v[i] = i;
            prime[++cnt] = i;
        }
        for(int j = 1; j<=cnt; j++)
        {
            if(prime[j]>v[i] || prime[j]>n/i)
                break;
            v[i*prime[j]] = prime[j];
        }
    }
}

3.欧几里得算法求gcd

int gcd(int a,int b)
{
    return b?gcd(b,a%b):a;
}

4.扩展欧几里得算法

定理：对于任意整数a,b,存在一对整数x,y,满足ax + by = gcd(a,b).

证明：

在欧几里得算法的最后一步，即b = 0时，显然有x = 1,y = 0，使得a*1+0*0 = gcd(a,0)

若b>0，则gcd(a,b) = gcd(b,a mod b)，假设存在一对整数x，y，满足b*x + (a mod b)*y = gcd(b,a mod b)

因为b*x + (a mod b)*y = b*x + (a - b * (a/b))*y = a*y+ b*(x - (a/b)*y)

所以令x' = y，y' = x - (a/b)*y，得a*x' + b*y' = gcd(a,b)

int exgcd(int a,int b,int &x, int &y)
{
    if(b == 0)
    {
        x = 1;
        y = 0;
        return a;
    }
    int d = exgcd(b,a%b,x,y);
    int z = x;
    x = y;
    y = z-y*(a/b);
    return d;
}

对于更一般的方程ax + by = c，当且仅当d|c，方程才有解，所以可以先求出ax + by = d的一组特解x0,y0,然后同时乘上c/d就可以了

求a模素数m的逆元x，即a*x = 1(mod m)，即m | a*x-1，所以只要解一下线性同余方程a*x + m*y = d就可以了（d = 1）

ll Inv(ll a)
{
    ll d,x,y;
    d=exgcd(a,mod,x,y);
    if(d==1)
      //此处将解平移到1~mod之间
      return (x%mod+mod)%mod;
    return -1;
}

5.质因数分解

常用的是结合试除法和Eratothenes筛法的分解

int p[maxn],c[maxn];
void divide(int n)
{
    cnt = 0;
    for(int i = 2; i<=sqrt(n); i++)
    {
        if(n%i == 0)
        {
            p[++cnt] = i;
            c[cnt] = 0;
            while(n%i == 0)
            {
                n /= i;
                c[cnt]++;
            }
        }
    }
    if(n>1)
    {
        p[++cnt] = n;
        c[cnt] = 1;
    }
}

6.欧拉函数

普通计算欧拉函数方法

int euler(int n)
{
    int m = sqrt(n) + 0.5;
    int ans = n;
    for (int i = 2; i <= m; ++i)
    if (n%i == 0)
    {
        ans = ans / i * (i - 1);
        while (n%i == 0) n /= i;
    }
    if (n > 1) ans = ans / n *(n - 1);
    return ans;
}

Eratosthenes筛法

复杂度 O(NlogN)。

筛法基于 对任意x的整数倍数，2x，3x，… 都不是质数。从2开始，从小到大扫描每个数X，把它的倍数2x，3x，…,floor(N/x)*x标记为合数。当扫描到一个数时，如果它尚未标记，则它一定不能被2~x-1之间的任何数整除，该数就是质数。

又因为实际上小于等于x^2的x的倍数在扫描更小的数时就已经被标记过了，因此，我们可以对Eratosthens筛法进行优化。对于每个数x，我们只需要从x^2开始，把x^2，(x+1)*x，…，floor（N/x）*x标记为合数即可。

int n;
ll phi[maxn+5];
void euler(int n)
{
    for(int i = 2; i<=n; i++)
        phi[i] = i;
    for(int i = 2; i<=n; i++)
    {
        if(phi[i] == i)
        {
            for(int j = i; j<=n; j += i)
                phi[j] = phi[j]/i*(i-1);
        }
    }
}

线性筛

复杂度O（N）

1.若p|n且p^2|n，则phi(n) = phi(n/p)*p

2.若p|n但p^2不能整除n，则phi(n) = phi(n/p)*(p-1)

在线性筛法中，每个合数n只会被它的最小质因子p筛一次。我们恰好可以在此时执行以上两条判断，从phi(n/p)递推到phi(n)

int v[maxn+5],prime[maxn+5];
ll phi[maxn+5];
void euler(int n)
{
    memset(v,0,sizeof v);
    cnt = 0;
    for(int i = 2; i<=n; i++)
    {
        //i是质数
        if(v[i] == 0)
        {
            v[i] = i;
            prime[cnt++] = i;
            phi[i] = i-1;
        }
        //给当前的数i乘上一个质因子
        for(int j = 0; j<cnt; j++)
        {
            //i有比prime[j]更小的质因子或超出n的范围
            if(prime[j]>v[i] || prime[j]>n/i)
                break;
            //prime[j]是i*prime[j]的最小质因子
            v[i*prime[j]] = prime[j];
            phi[i*prime[j]] = phi[i]*(i%prime[j] ? prime[j]-1:prime[j]);
        }
    }



}

7、矩阵快速幂求斐波那契数列

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
ll MOD = 1e9+7;
ll a1,a2,a3,n;
struct sq{
    ll sq[2][2];
    void mes(){
        memset(sq,0,sizeof(sq));
    }
};
sq operator *(sq a,sq b){
    sq tmp;
    tmp.mes();
    for(int i=0;i<2;i++){
        for(int j=0;j<2;j++){
            for(int k=0;k<2;k++){
                tmp.sq[i][j]=(tmp.sq[i][j]+a.sq[i][k]*b.sq[k][j])%MOD;
            }
        }
    }
    return tmp;
}
ll fi(ll n){
    sq ans,m;
    ans.mes();
    m.mes();
    m.sq[0][0]=m.sq[0][1]=m.sq[1][0]=1;
    m.sq[1][1]=0;
    ans.sq[0][0]=ans.sq[1][1]=1;
    ans.sq[0][1]=ans.sq[1][0]=0;
    while(n){
        if(n&1){
            ans=ans*m;
        }
        m=m*m;
        n>>=1;
    }
    return ans.sq[0][1]%MOD;
}

8.快速求组合数

#include<iostream>
#include<cstring>
#include<algorithm>
#include<cstdio>
#include<string>
#include<cmath>
#include<vector>
#include<stack>
#include<deque>
#include<map>
#include<set>
#include<queue>
typedef long long ll;
inline void sd(int &n) { scanf("%d", &n); }
const double pi = acos(-1);
const int inf = 0x3f3f3f3f;
using namespace std;
const int maxn = 1e6 + 6;

const ll MOD = 1e9 + 7;

ll fac[maxn];
ll inv[maxn];

ll Qpow(ll x, ll n){
	ll ret = 1;
	ll tmp = x % MOD;
	
	while(n){
		if(n & 1) ret = (ret * tmp) % MOD;
		tmp = (tmp * tmp) % MOD;
		n >>= 1;
	}
	return ret;
}

void init(){
	fac[0] = 1;
	for(int i = 1;i < maxn;++i){
		fac[i] = fac[i - 1] * i % MOD;
	}
	inv[maxn - 1] = Qpow(fac[maxn - 1], MOD - 2);
	for(int i = maxn - 2;i >= 0;i--){
		inv[i] = inv[i + 1] * (i + 1) % MOD;
	}
}

ll C(ll a, ll b){
	if(b > a){
		return 0;
	}
	if(b == 0) return 1;
	return fac[a] * inv[b] % MOD * inv[a - b] % MOD;
}