数论-模板 (待更新)

质因数分解

在求解数论问题时, 经常会把一个整数分解为多个质数的乘积(形式唯一),
还有一些题目需要从它的质因数分解形式来考虑.

质因数分解:

$$x = \prod_{i=1}^np_i^{a_i}(p_i为不可分解的质因数)$$

eg:
$12 = 2 * 2 * 3 = 2^2*3$
注: 不含的质因子可认为其指数是0

实现1:

void get(int x)
{
    for(int i = 2; i * i <= x; ++i)//从2开始枚举因子
        while(x % i == 0)//含有质因子i时, 可以把x分解中i的指数消去
        {
            x /= i;   
            ++num[i];    //标记当前质因子的指数
        }
    if(x > 1) ++num[x]; //剩下的质因子指数为0/1
}

实现2:(用线性筛筛出预选区间的质因子再分解)

const int maxn = 1010;
int cnt;
int prime[maxn];
bool isPrime[maxn];
void Euler_Seive()
{
    memset(isPrime, true, sizeof isPrime);
    isPrime[1] = 0;
    for(int i = 2; i < maxn; ++i)
    {
        if(isPrime[i]) prime[++cnt] = i;
        for(int j = 1; j <= cnt && i * prime[j] < maxn; ++j)
        {
            isPrime[i * prime[j]] = 0;
            if(i % prime[j] == 0) break;
        }
    }
}
void get(int x)
{
    Euler_Seive();
    for(int i = 1; prime[i] * prime[i] <= x; ++i)
        while(x % prime[i] == 0)
        {
            x /= prime[i];   
            ++num[prime[i]];    
        }
    if(x > 1) ++num[x]; 
}

整除

若a % b = 0, 则a能被b整除或者b是a的因子(约数)记为b | a

质因子角度考虑a能被b整除:

a中所有质因子且指数不小于其在b中的指数(不包含的质因子认为其指数为0)

x的约数不能含有其他质因子且其$p_i$指数范围为$[0, a_i]$
即:
$任意d | x, d = \prod_{i=1}^np_i^{t_i}(t_i <= a_i)$

一个正整数x的约数个数f(x):

$x = \prod_{i=1}^np_i^{a_i}$
$f(x) = \prod_{i=1}^n(a_i + 1)$

---

最大公约数和最小公倍数

a, b的最大公约数记gcd(a, b)

$gcd(a, b) = Max(D)(任意d \in D,d | a, d | b)$
eg:

a = 12, b = 4

D = {1, 2, 4}

gcd(12, 4) = Max(D) = 4

---

a, b的最小公约数记lcm(a, b)

$lcm(a, b) = Min(M)(任意d \in M,a | d, b | d)$

eg:

a = 12, b = 4

M = {12, 24, 36….}

lcm(12, 4) = Min(M) = 12

---

质因子角度考虑gcd(a, b)、lcm(a, b)

设:

$a = \prod_{i=1}^np_i^{a_i}$

$b = \prod_{i=1}^np_i^{b_i}$

注:展开形式按质数全体来看, 不含的质因子指数为0

则有定义可得

$任意d \in D, d = \prod_{i=1}^np_i^{d_i}(d_i<=a_i\bigwedge d_i<=b_i)$

即:

$任意d \in D, d = \prod_{i=1}^np_i^{d_i}(d_i<=Min(a_i,b_i))$

$gcd(a, b) = Max(D) = \prod_{i=1}^np_i^{Min(a_i,d_i)}$

同理:

$任意d \in M, d = \prod_{i=1}^np_i^{d_i}(d_i>=a_i\bigwedge d_>=b_i)$

即:

$任意d \in M, d = \prod_{i=1}^np_i^{d_i}(d_i>=Max(a_i,b_i))$

$lcm(a, b) = Min(M) = \prod_{i=1}^np_i^{Max(a_i,d_i)}$

同时:

a * b = gcd(a, b) * lcm(a, b)

---

实现gcd(a, b)和lcm(a, b)

两个不全为0的自然数a, b, 求gcd(a, b)

设d = gcd(a, b), 则$d | a \bigwedge d | b$

则d | (a - b)

推广至d | (a - k * b)

且可证gcd(a - k * b, b) = d

又因为a % b = a - [a / b] * b

所以gcd(a, b) = gcd(a % b, b) = gcd(b, a % b)

实现:

int gcd(int a, int b) 
{ 
    //辗转相除出现0时gcd(a, 0) = a
    return b == 0 ? a : gcd(b, a % b); 
}
int lcm(int a, int b)
{
    //上面结论, 先除防溢出
    return a / gcd(a, b) * b;
}

---

扩展欧几里得(求解ax + by = gcd(a, b))

不予介绍具体原理

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

---

矩阵快速幂

const int maxn = 10;
struct Matrix
{
    int mat[maxn][maxn];
    int row, col;
    Matrix(int _row = maxn, int _col = maxn)
    {
        row = _row;
        col = _col;
    }
    void zero()
    {
        memset(mat, 0, sizeof mat);
    }
    void E()
    {
        for(int i = 0; i < maxn; ++i) mat[i][i] = 1;
    }
    Matrix operator * (const Matrix &rhs) const
    {
        Matrix res(row, rhs.col);
        res.zero();
        for(int i = 0; i < row; ++i)
            for(int j = 0; j < rhs.col; ++j)
                for(int k = 0; k < col; ++k)
                    res.mat[i][j] += mat[i][k] * rhs.mat[k][j];
        return res;
    }
    Matrix operator ^ (int b) const
    {
        Matrix res(row, col);
        res.E();
        Matrix a(row, col);
        memcpy(a.mat, mat, sizeof mat);
        while(b)
        {
            if(b & 1) res = res * a;
            a = a * a;
            b >>= 1;
        }
        return res;
    }
}

---

容斥定理

在对一个集合多个约束计数时, 简单的相加可能会出现重叠, 这时就需要容斥计数.(求加和并把重叠的部分消去)

概率论教材中中有提到:
$|A_1 \bigcup A_2| = |A_1| + |A_2| - |A_1 \bigcap A_2|$

这就是最简单的容斥.(n = 2时)

可用数学归纳法证明容斥定理:

$$|\bigcup _{i = 1} ^ n A_i| = \sum_{i=1}^n (-1)^{i-1} \sum_{1<=A_1<=A_2···<=A_i<=n}|A_1 \bigcap A_2 ··· \bigcap A_i |$$

---

容斥定理常用预处理

1. 容斥计数时枚举的约束一定是当前计数约束全集的非空子集, 我们可以用状态表示, 若n = 5, 用$[1, 2^5)表示某一非空子集$

2. 预处理子集(适用于搜索中):
   就是求一个数表示状态的子状态

   eg: 7 = 111(B)表示待求的$|A_1 \bigcup A_2 \bigcup A_3|$(约束从1开始计数)

   则约束全集是$\{A_1, A_2, A_3\}$, 那么其非空子集集合是

   $\{A_1, A_2, A_3, A_1 \bigcap A_2, A_1 \bigcap A_3, A_2 \bigcap A_3, A_1 \bigcap A_2 \bigcap A_3\}$

   其状态表示分别书1 = 001B, 2 = 010, 4 = 100B, 3 = 011B, 5 = 101B, 6 = 110B, 7 = 111B.

   可以预处理出7的子状态

int bitcnt[maxn];//二进制1的个数
vector<int> vec[maxn];
inline int lowbit(int x)
{
    return x & (-x);
}
void init()
{
    bitcnt[0] = 0; 
    for(int i = 1; i < max_m; ++i) bitcnt[i] = bitcnt[i ^ lowbit(i)] + 1;
    for(int i = 1; i < max_m; ++i)
    {
        vector<int> bit_pos;//二进制各1的位置
        for(int j = 0; j < 15; ++j) if(i & (1 << j)) bit_pos.push_back(j);
        for(int j = 0; j < 1 << bitcnt[i]; ++j)
        {
            int x = 0;
            for(int k = 0; k < bitcnt[i]; ++k)
                if(j & (1 << k)) x |= 1 << bit_pos[k];
            vec[i].push_back(x);//包含空集0, 视题目而定
        }
    }
}

---

例题:
求[1, r]区间内和x互质的数有多少

区间有r个数, 可以转换成求区间内和x不互质的数的个数.

对于任意a属于[1, r], 若gcd(a, x) != 1则gcd(a, x)必定包含x的某些质

因子(指数至少为1)即a必定包含x的某些质因子.

设$x = \prod_{i=1}^np_i^{a_i}$

则令$A_i$表示[1, r]中包含x第i个质因子的数的个数

则$A_i = \lfloor \frac{r}{p_i} \rfloor, A_i \bigwedge A_j = \lfloor \frac{r}{p_ip_j}\rfloor$(取整不可省去)

根据容斥定理可解

code:

int prime[1010];
int cnt;
void get(int x)
{
    for(int i = 2; i * i <= x; ++i)
        if(x % i == 0)
        {
            while(x % i == 0) x /= i;
            prime[cnt++] = i;
        }
    if(x > 1) prime[cnt++] = x;
}
int cal(int r)
{
    int m = 0;
    for(int i = 1; i < 1 << cnt; ++i)       //子状态
    {
        int bitcnt = 0;
        int tmp = 1;
        for(int j = 0; j < cnt; ++j)
            if(i & (1 << j))
            {
                ++bitcnt;
                tmp *= prime[j];
            }
        m += bitcnt & 1 ? r / tmp : r / tmp * -1;
    }
    return r - m;
}

---

欧拉函数

$对于一个正整数x, x的欧拉函数指不超过x的并且和x互质的正整数的个数, 记为\phi(x)$

设$x = \prod_{i=1}^np_i^{a_i}$

有上文容斥定理例题可得

$$\phi(x) = x-(\lfloor\frac{x}{p_1}\rfloor + \lfloor\frac{x}{p_2} \rfloor+ ···\lfloor\frac{x}{p_n}\rfloor - \lfloor\frac{x}{p_1*p_2}\rfloor - \lfloor\frac{x}{p_1*p_3}\rfloor - \lfloor\frac{x}{p_2*p_3}\rfloor···+(-1)^{n-1}* \lfloor\frac{x}{p_1*p_2*···p_n} \rfloor )$$

由于$p_1···p_n$都是x的质因子, 必定能整除, 取整符号可以略去

$$\phi(x) = x-\frac{x}{p_1} -\frac{x}{p_2} - ···\frac{x}{p_n} + \frac{x}{p_1*p_2} + \frac{x}{p_1*p_3} + \frac{x}{p_2*p_3}···+(-1)^n* \frac{x}{p_1*p_2*···p_n}$$

因此$\phi(x) = x * \prod_{i=1}^n(1-\frac{1}{p_i})$(可用排列组合证明)

特殊:$\phi(1) = 1$

---

欧拉函数性质

1. 若正整数x是质数, 则$\phi(x) = x- 1$(公式可证)

2. 对于任意正整数x, 满足: $x = \sum_{d | x} \phi(d)$(质因数分解角度可证)

3. 积性函数: $若a,b互质, \phi(ab) = \phi(a)*\phi(b)$

---

欧拉函数实现

//求正整数x的欧拉函数值
int getPhi(int x)
{
    int res = x;
    for(int i = 2; i * i <= x; ++i)
        if(x % i == 0)
        {
            while(x % i == 0) x /= i;
            res -= res / i;// (1 - 1 / p_i)
        }
    if(x > 1) res -= res / x;
    return res;
}
//线性筛之后求正整数x的欧拉函数(略)

//线性筛求一个区间[1, maxn]内数的欧拉函数
int prime[maxn + 1];
bool isPrime[maxn + 1];
int phi[maxn + 1];
int cnt;
void Euler_Seive()
{
    memset(isPrime, true, sizeof isPrime);
    isPrime[1] = 0;
    phi[1] = 1;
    for(int i = 2; i <= maxn; ++i)
    {
        if(isPrime[i])
        {
            prime[++cnt] = i;
            phi[i] = i - 1;//质数欧拉函数
        }
        for(int j = 1; j <= cnt && i * prime[j] <= maxn; ++j)
        {
            isPrime[i * prime[j]] = 0;
            if(i % prime[j] == 0)
            {
                phi[i * prime[j]] = phi[i] * prime[j];//根据欧拉函数公式,每个质因子只需处理一次, i包含第j个质因子时不需再处理
                break;
            }
            phi[i * prime[j]] = phi[i] * (prime[j] - 1);//新质因子,i和prime[j]互质,应用欧拉函数积性
        }
    }
}
//枚举筛法
void init1()
{
    memset(phi, 0, sizeof phi);
    phi[1] = 1;
    for(int i = 2; i <= maxn; ++i)
    {
        if(phi[i] == 0)
            for(int j = i; j <= maxn; j += i)
            {
                if(phi[j] == 0) phi[j] = j;
                phi[j] = phi[j] / i * (i - 1);
            }
    }
}

//根据性质(约数函数和为本身)
void init2()
{
    for(int i = 1; i <= maxn; ++i) phi[i] = i;
    for(int i = 1; i <= maxn / 2; ++i)
        for(int j = i + i; j <= maxn; j += i)
            phi[j] -= phi[i];
}

---

莫比乌斯函数

定义：

$\mu(1) = 1$

对于n > 1, 若n是质数$\mu(n) = -1$

若n含有平方因子, $\mu(n) = 0$

若n含是a个质因子的乘积, $\mu(n) = (-1)^a$

---

莫比乌斯函数性质

1. 积性函数: 若a, b互质$\mu(ab)=\mu(a)*\mu(b)$(自证)
2. $\sum_{d|n}\mu(d) = [n = 1] = n == 1 ? 1 : 0$质因子排列组合可证

---

筛莫比乌斯函数

//线性筛求一个区间[1, maxn]内数的莫比乌斯！函数
int prime[maxn + 1];
bool isPrime[maxn + 1];
int mobius[maxn + 1];
int cnt;
void Euler_Seive()
{
    memset(isPrime, true, sizeof isPrime);
    isPrime[1] = 0;
    mobius[1] = 1;
    for(int i = 2; i <= maxn; ++i)
    {
        if(isPrime[i])
        {
            prime[++cnt] = i;
            mobius[i] =  -1;//质数莫比乌斯函数
        }
        for(int j = 1; j <= cnt && i * prime[j] <= maxn; ++j)
        {
            isPrime[i * prime[j]] = 0;
            if(i % prime[j] == 0)
            {
                mobius[i * prime[j]] = 0;//有平方因子
                break;
            }
            mobius[i * prime[j]] = -mobius;//(-1)^a,积性函数
        }
    }
}

//根据性质(约数莫比乌斯函数和为元函数)
//元函数:f(e) = [e = 1] = e == 1 ? 1 : 0
void init()
{
    memset(mobius, 0, sizeof mobius);//初始化0
    mo[1] = 1;
    for(int i = 1; i <= maxn / 2; ++i)
        for(int j = i + i; j <= maxn; j += i)
            mobius[j] -= mobius[i];
}

---

莫比乌斯反演

上tls干货:

http://blog.csdn.net/skywalkert/article/details/50500009

重点是狄利克雷卷积和莫比乌斯反演, 积性函数求和(涉及到杜教筛)前置技能点满后再看.

例题:HDU 1695

---

莫比乌斯反演技巧

在莫比乌斯反演中必不可少的是预处理一个区间的因子(狄利克雷卷积)的莫比乌斯函数

//先线性筛mobius函数
//然后是预处理因子
//[1, maxn]
vector<int> fac[maxn + 10];
void getFactor()
{
    for(int i = 1; i <= maxn; ++i)
        if(mobius[i] == 0) continue;    //反演中对求和无贡献
        else for(int j = i; j <= maxn; j += i)
                fac[j].push_back(i);
}

---

逆元

若a, p互质即gcd(a, p) = 1

则存在x 使得a * x % p = 1(扩展欧几里得定理)

即$x = a^{-1}, 读作a关于p的逆元$

---

逆元用在哪

取模运算时出现了除法就可以转成乘除数的逆元

eg:

求a / b % p

已知a % b = 0, 且gcd(b, p) = 1;

设$t = b^{-1}, c = a / b$

则b * t % p = 1, 答案为 c * b * t % p = c % p.

---

求逆元

1. 扩展欧几里得(ax + py = 1, gcd(a, p) = 1)解得的x即为a关于p的一个逆元
2. 欧拉定理: $gcd(a, p) = 1 -> a^{\phi(p)} \mod p = 1, a^{-1} = a^{\phi(p) - 1}$用快速幂取模求解
3. 求一个区间[1, maxn]关于p的逆元(p一般是质数)

void getInv()
{
    inv[1] = 1;
    for(int i = 2; i <= maxn; ++i)
        inv[i] = (p - p / i) * inv[p % i] % p;

    //inv[p%i]*(p%i) % p = 1
    //inv[p%i]*(p-p/i*i)%p=1
    //i*(-p/i)*inv[p%i]=1
    //inv[i] = (p-p/i)*inv[p%i]%p
}

---

原根

干货:

http://blog.csdn.net/acdreamers/article/details/8883285

原根用法:

对于一个质数p, [1, p - 1]内的数都可以用p的原根的幂表示.

有时需要从这个角度分析问题

---

BM算法(线性递推式)

#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
#include <vector>
#include <string>
#include <map>
#include <set>
#include <cassert>
using namespace std;
#define rep(i,a,n) for (int i=a;i<n;i++)
#define per(i,a,n) for (int i=n-1;i>=a;i--)
#define pb push_back
#define mp make_pair
#define all(x) (x).begin(),(x).end()
#define fi first
#define se second
#define SZ(x) ((int)(x).size())
typedef vector<int> VI;
typedef long long ll;
typedef pair<int,int> PII;
const ll mod=1000000007;
ll powmod(ll a,ll b) {ll res=1;a%=mod; assert(b>=0); for(;b;b>>=1){if(b&1)res=res*a%mod;a=a*a%mod;}return res;}

int _,n;
namespace linear_seq {
    const int N=10010;
    ll res[N],base[N],_c[N],_md[N];

    vector<int> Md;
    void mul(ll *a,ll *b,int k) {
        rep(i,0,k+k) _c[i]=0;
        rep(i,0,k) if (a[i]) rep(j,0,k) _c[i+j]=(_c[i+j]+a[i]*b[j])%mod;
        for (int i=k+k-1;i>=k;i--) if (_c[i])
            rep(j,0,SZ(Md)) _c[i-k+Md[j]]=(_c[i-k+Md[j]]-_c[i]*_md[Md[j]])%mod;
        rep(i,0,k) a[i]=_c[i];
    }
    int solve(ll n,VI a,VI b) {
        ll ans=0,pnt=0;
        int k=SZ(a);
        assert(SZ(a)==SZ(b));
        rep(i,0,k) _md[k-1-i]=-a[i];_md[k]=1;
        Md.clear();
        rep(i,0,k) if (_md[i]!=0) Md.push_back(i);
        rep(i,0,k) res[i]=base[i]=0;
        res[0]=1;
        while ((1ll<<pnt)<=n) pnt++;
        for (int p=pnt;p>=0;p--) {
            mul(res,res,k);
            if ((n>>p)&1) {
                for (int i=k-1;i>=0;i--) res[i+1]=res[i];res[0]=0;
                rep(j,0,SZ(Md)) res[Md[j]]=(res[Md[j]]-res[k]*_md[Md[j]])%mod;
            }
        }
        rep(i,0,k) ans=(ans+res[i]*b[i])%mod;
        if (ans<0) ans+=mod;
        return ans;
    }
    VI BM(VI s) {
        VI C(1,1),B(1,1);
        int L=0,m=1,b=1;
        rep(n,0,SZ(s)) {
            ll d=0;
            rep(i,0,L+1) d=(d+(ll)C[i]*s[n-i])%mod;
            if (d==0) ++m;
            else if (2*L<=n) {
                VI T=C;
                ll c=mod-d*powmod(b,mod-2)%mod;
                while (SZ(C)<SZ(B)+m) C.pb(0);
                rep(i,0,SZ(B)) C[i+m]=(C[i+m]+c*B[i])%mod;
                L=n+1-L; B=T; b=d; m=1;
            } else {
                ll c=mod-d*powmod(b,mod-2)%mod;
                while (SZ(C)<SZ(B)+m) C.pb(0);
                rep(i,0,SZ(B)) C[i+m]=(C[i+m]+c*B[i])%mod;
                ++m;
            }
        }
        return C;
    }
    int gao(VI a,ll n) {
        VI c=BM(a);
        c.erase(c.begin());
        rep(i,0,SZ(c)) c[i]=(mod-c[i])%mod;
        return solve(n,c,VI(a.begin(),a.begin()+SZ(c)));
    }
};

int main() {
    for (scanf("%d",&_);_;_--) {
        scanf("%d",&n);
        printf("%d\n",linear_seq::gao(VI{2,24,96,416,1536,5504,18944,64000,212992,702464},n-1));
    }
}

---

其它数论算法

1. lucas定理及其扩展
2. BSGS及其扩展
3. 同余方程组及中国剩余定理
4. 杜教筛
5. 组合数相关

---