数学总结

发布日期: 2022-06-06

更新日期: 2024-06-26

文章字数: 6.8k

OI模板-数学

注意，本文不同于 OI数学总结-数论或其他总结类文章

本文倾向于代码实现，只是因为数论的内容与代码联系紧密，并且数论内容较多，因此显得较为相似。

最大公约数 gcd

最大公约数指能够整除多个不全为零的整数的最大正整数

常用性质： $a\times b = \gcd(a,b) \times \mathrm{lcm}(a,b)$

int gcd(int a,int b)
{
    return b==0?a:gcd(b,a%b);
} // 注意这份代码对于(+,-)的情况会求出负数

扩展欧几里德算法(Exgcd)

板子题： P5656 【模板】二元一次不定方程 (exgcd)

若 $b \ne 0$ ，扩展欧几里得算法求出的可行解必有 $|x| \le b ,|y| \le a$ 。

int exgcd(int a,int b,int &x,int &y)
{
    int d=a;
    if(!b)x=1,y=0;
    else d=exgcd(b,a%b,y,x),y-=a/b*x;
    return d;
}
void solve()
{
    int a,b,c,x,y;
    read(a);read(b);read(c);
    int d=exgcd(a,b,x,y);
    if(c%d!=0){puts("-1");return;} // 无解
    x*=c/d;y*=c/d;
    int p=b/d,q=a/d,k;
    if(x<0)k=ceil((1.0-x)/p),x+=p*k,y-=q*k; // 将x提高到最小正整数
    else k=(x-1)/p,x-=p*k,y+=q*k; // 将x降低到最小正整数
    if(y>0) // 存在x>0,y>0的解（正整数）
    {
        write((y-1)/q+1);   pc(' '); // 将y减到1的方案数即为解的个数
        write(x);           pc(' '); // 当前的x为最小正整数
        write((y-1)%q+1);   pc(' '); // 将y取到最小正整数
        write(x+(y-1)/q*p); pc(' '); // 将x提升到最大
        write(y);           pc(' '); // 特解即为y的最大值
    }else // 仅存在整数解
    {
        write(x);                       pc(' '); // 当前的x
        write(y+q*(int)ceil((1.0-y)/q));pc(' '); // 将y提高到正整数
    }
    pc('\n');
}

逆元

如果一个线性同余方程组 $ax \equiv 1 \pmod{b}$ ，则 $x$ 称为 $a \bmod b$ 的逆元，记作 $a^{-1}$

模意义下的乘法逆元不唯一。如果 $b$ 不为素数，可能不存在逆元。

模板题：P3811 【模板】乘法逆元

求解逆元有多重方法，我们一个个来讲。

第一种方法：$O(n)$ 预处理，$O(1)$ 查询

求的是 $1,2,\cdots ,n$ 的模 $p$ 意义下的乘法逆元，$p$ 一般为质数。

为什么说一般呢？因为「 $p$ 为合数」会导致有些时候不存在逆元。

如果正好存在逆元也是可以继续做的。

原理：

$\begin{aligned} i^{-1} &\equiv \begin{cases} 1,&\text{if }i=1,\\\\ -\left\lfloor\dfrac{p}{i}\right\rfloor \left(p\bmod{i}\right)^{-1}, &\text{otherwises.} \end{cases} \pmod{p} \end{aligned}$

注意空间复杂度为 $O(n)$

inv[1] = 1;
for(int i = 2; i <= n; i++)
    inv[i] = (p - p / i) * inv[p % i] % p;

第二种方法：费马小定理 $O(\log p)$ 查询 $x$ 的逆元

原理：$x \equiv a^{p-2}\pmod{p}$ （根据费马小定理），$p$ 需要是质数

cout << qpow(a,p-2) << '\n'; // 快速幂

第三种方法：$\text{Exgcd}~O(\log p)$ 查询 $x^{-1}$

求的是 $a$ 在模 $p$ 意义下的逆元，要求模数 $p$ 与 $a$ 互质。

int exgcd(int a,int b,int &x,int &y)
{
    int d = a;
    if(!b) { x = 1,y = 0; }
    else { d = exgcd(b, a % b, y, x), y -= a / b * x; }
    return d;
}
int solve(int a,int p)
{
	int x,y;
    exgcd(a,p,x,y);
    return (x = (x % p + p) % p);
}

第四种方法（特殊）：$\mathcal{O}(n)$ 预处理阶乘的逆元

因为 $n!=(n-1)! \times n$ ，所以 $((n-1)!)^{-1}=(n!)^{-1} \times n$ 。模数 $p$ 为质数。

for(int i=1; i<=n; i++) fac[i] = fac[i-1] * i % p;
invf[n] = qpow(fac[n], p-2); // 费马小定理
for(int i=n-1; i>=1; i--) invf[i] = invf[i+1] * (i+1) % p;

请注意，如果 $p \le n$ ，这种方法就不可行了，因为 $n! \bmod p =0$ ，会出现不存在逆元的情况

这种时候一般采用上文线性求逆元的方法，只要求 $1\sim p-1$ 的就行了，即

// luogu P2675
for(int i=1; i<=n; i++) fac[i] = fac[i-1] * i % p;
for(int i=2; i<p; i++) inv[i] = (p-p/i) * inv[p%i] % p; // 注意是1 ~ p-1

欧拉函数 $\varphi$

欧拉函数 $\varphi(n)$ 表示小于等于 $n$ 和 $n$ 互质的数的个数。

设 $n=\prod_{i=1}^{s}p_i^{k_i}$

$\varphi(n) = n \prod\limits_{i=1}^s\left(\dfrac{p_i-1}{p_i}\right) = \prod(p_i-1)\times p_i^{k_i-1}$

常用性质：$\varphi(1) = 1,~\varphi(p) = p-1,~p\in\mathbb{P}$

模板题：poj2407 。套路地，欧拉函数也有多种求解方法。

第一种方法：$O(\sqrt{n})$ 求 $\varphi(n)$

这种做法就是根据欧拉函数的定义，模拟计算。

int Euler_phi(int n)
{
    int ans=n;
    for(int i=2; i<=n/i; 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;
}
signed main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);cout.tie(0);
    // freopen("check.in","r",stdin);
    // freopen("check.out","w",stdout);
    int n;
    while(cin>>n&&n!=0)
        cout << Euler_phi(n) << endl;

    return 0;
}

第二种方法：线性筛 $O(n)$ 求解

可能是模板题：~~poj2478 （这个题要改成 $\varphi$ 的前缀和）~~

int phi[N],prime[N],pcnt;
bool ck[N];
void Euler()
{
    ck[1]=1;
    phi[1]=1;
    for(int i=2; i<=N-5; i++)
    {
        if(!ck[i])
        {
            prime[++pcnt]=i;
            phi[i]=i-1;
        }
        for(int j=1; j<=pcnt&&i*prime[j]<=N-5; j++)
        {
            int pos=i*prime[j];
            ck[pos]=1;
            if(i%prime[j])
            {
                phi[pos]=phi[i]*phi[prime[j]];
            }
            else
            {
                phi[pos]=phi[i]*prime[j];
                break;
            }
        }
    }
}

扩展欧拉定理

$\begin{aligned} a^b \equiv \begin{cases} a^b,&b< \varphi(m)\\\\ a^{b \ \bmod \ \varphi(m)+\varphi(m)},&b\ge \varphi(m) \end{cases} \pmod{m} \end{aligned}$

例题：

求

$a^b\bmod m$

$1\le a,m\le 10^{12},1\le b\le 10^{20000000}$（这个是加强版的，朴素版只要int就行）

#include <bits/stdc++.h>
using namespace std;
#define int long long
#define INF 0x3f3f3f3f3f3f3f3f
#define gc() getchar()
int a,m,b,flag;
int Euler_phi(int n)
{
    int ans=n;
    for(int i=2; i<=n/i; 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;
}
int qpow(int a,int b,int p)
{
    int ans=1,base=a%p;
    while(b)
    {
        if(b&1)ans=(__int128)ans*base%p;
        base=(__int128)base*base%p;
        b>>=1;
    }
    return ans;
}
signed main()
{
    // freopen("check.in","r",stdin);
    // freopen("check.out","w",stdout);
    scanf("%lld%lld",&a,&m);
    a%=m;
    int phi=Euler_phi(m);
    char ch=gc();
    while(!isdigit(ch))ch=gc();
    while(isdigit(ch))
    {
        b=(b<<1)+(b<<3)+(ch^48);
        if(b>=phi)flag=1,b%=phi;ch=gc();
    }
    if(b>=phi)flag=1,b%=phi;
    if(flag)b+=phi;
    printf("%lld\n",qpow(a,b,m));
    return 0;
}

Lucas定理

对于质数 $p$ ，有

$\dbinom{n}{m} \bmod p = \dbinom{\left\lfloor{n/p}\right\rfloor}{\left\lfloor{m/p}\right\rfloor} \times \dbinom{n \bmod p}{m \bmod p} \bmod p$

$\binom{\left\lfloor{n/p}\right\rfloor}{\left\lfloor{m/p}\right\rfloor}$ 采用递归处理，后面的直接 $O(p)$ 预处理组合数

模板题：P3807 【模板】卢卡斯定理/Lucas 定理

#include <bits/stdc++.h>
using namespace std;
#define int long long
#define INF 0x3f3f3f3f3f3f3f3f
#define N (int)(2e5+15)
int Q,n,m,p;
int fac[N];
int qpow(int a,int b,int p)
{
	int ans=1,base=a;
	while(b)
	{
		if(b&1)ans=ans*base%p;
		base=base*base%p;
		b>>=1;
	}
	return ans;
}
// n 选 m
int C(int n,int m,int p)
{
	if(n<m)return 0;
	return fac[n]*qpow(fac[m],p-2,p)%p*qpow(fac[n-m],p-2,p)%p;
}
int lucas(int n,int m,int p)
{
	if(!m)return 1;
	return lucas(n/p,m/p,p)*C(n%p,m%p,p)%p;
}
signed main()
{
	cin >> Q;
	while(Q--)
	{
		cin >> n >> m >> p;
		fac[1]=fac[0]=1;
		for(int i=2; i<=n+m; i++)
			fac[i]=fac[i-1]*i%p;
		cout << lucas(n+m,m,p) << '\n'; // C(n+m,n)
	}
	return 0;
}

中国剩余定理

参考中国剩余定理 & 扩展

中国剩余定理(CRT)可求解如下形式的一元线性同余方程组

$\begin{cases} x \equiv a_1\ \pmod{p_1} \\[6pt] x\equiv a_2\ \pmod{p_2} \\[6pt] \quad\vdots \\[6pt] x \equiv a_n\ \pmod{p_k}\end{cases}$

其中 $p_1,p_2,\cdots,p_k$ 两两互质。

板子题：P1495 【模板】中国剩余定理（CRT）/ 曹冲养猪

时间复杂度 $\mathcal{O}(n \log V)$

代码：

#include <bits/stdc++.h>
using namespace std;
#define int long long
#define INF 0x3f3f3f3f3f3f3f3f
void up(int &x,int y) { x < y ? x = y : 0; }
void down(int &x,int y) { x > y ? x = y : 0; }
#define N ((int)(15))

int a[N],p[N];
void exgcd(int a, int b, int &x, int &y)
{
	if(!b) { x = 1, y = 0; }
	else exgcd(b, a % b, y, x), y -= a / b * x;
}
int solve(int n, int P)
{
	int res = 0;
	for(int i = 1; i <= n; i++)
	{
		int m = P / p[i], x, y;
		exgcd(m, p[i], x, y); // m * x % p[i] = 1;
		res = (res + a[i] * m % P * x % P) % P;
	}
	return (res % P + P) % P;
}
signed main()
{
	ios::sync_with_stdio(0);
	cin.tie(0);cout.tie(0);
	// freopen("check.in","r",stdin);
	// freopen("check.out","w",stdout);
	int n, P = 1; cin >> n;
	for(int i = 1; i <= n; i++) {
		cin >> p[i] >> a[i]; P = P * p[i];
	}
	cout << solve(n, P) << '\n';
	return 0;
}

扩展中国剩余定理

ExCRT 用于解决 $p_i$ 不一定互质的情况

板子题：P4777 【模板】扩展中国剩余定理（EXCRT）

时间复杂度 $\mathcal{O}(n \log V)$

代码：

// 2024年02月09日 03:38:48
#include <bits/stdc++.h>
using namespace std;
#define int long long
#define INF 0x3f3f3f3f3f3f3f3f
void up(int &x,int y) { x < y ? x = y : 0; }
void down(int &x,int y) { x > y ? x = y : 0; }
#define N ((int)(1e5 + 15))

int exgcd(int a, int b, int &x, int &y)
{
    int d = a;
    if(!b) { x = 1, y = 0; }
    else { d = exgcd(b, a % b, y, x), y -= a / b * x;}
    return d;
}
int a[N], p[N];
int ExCRT(int n)
{
    int u, v, p1 = p[1], res = a[1];
    for(int i = 2; i <= n; i++)
    {
        int p2 = p[i], c = (a[i] - res % p2 + p2) % p2;
        int g = exgcd(p1, p2, u, v); assert(c % g == 0);
        int l = p2 / g; u = (__int128_t) u * (c / g) % l;
        res += u * p1; p1 *= l; res = (res % p1 + p1) % p1;
    }
    return (res % p1 + p1) % p1;
}
signed main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);cout.tie(0);
    // freopen("check.in","r",stdin);
    // freopen("check.out","w",stdout);
    int n; cin >> n;
    for(int i = 1; i <= n; i++) cin >> p[i] >> a[i];
    cout << ExCRT(n) << '\n';
    return 0;
}

莫比乌斯函数 $\mu$

$\mu$ 为莫比乌斯函数，定义为

$\begin{aligned} \mu(n) = \begin{cases} 1,&n=1\\[8pt] (-1)^k,&k\text{ 为 }n \text{ 本质不同质因子个数}\\[8pt] 0,&n\text{含有平方因子（其实就是剩余的情况）} \end{cases} \end{aligned}$

一般线性筛 $\mathcal{O}(n)$ 求 $\mu$ ，用于莫比乌斯反演，参考 OI数学总结-数论

// 2023年07月04日 20:43:15
void Mobius(int n = N - 10)
{
    mu[1] = 1;
    for(int i = 2; i <= n; i++)
    {
        if(!ck[i]) { prime[++pcnt] = i; mu[i] = -1; }
        for(int j = 1; j <= pcnt && i * prime[j] <= n; j++)
        {
            int pos = i * prime[j]; ck[pos] = true;
            if(i % prime[j]) { mu[pos] = -mu[i]; } else { mu[pos] = 0; break; }
        }
    }
}

排列数

一般排列数我们都利用 $A^m_n = m! \times C^m_n$ 计算，下面我们讲一个特殊的计算方法。

排列 $A_n^m$ 直接算可以 $O(n)$

例如求解 $A_{n-m+1}^{m} \bmod p$

int res=1;
for(int i=n-m+1; i>=n-2*m+2; i--)
    res=res*i%p;
cout << res << '\n'

组合数

从 $n$ 个不同元素中取 $k(k \le n)$ 个不同元素的方案数

$\binom{n}{k} = \dfrac{n!}{k!(n-k)!}$

高中一般记作 $\mathrm{C}_n^{k}$ ，不是很推荐使用。

特别地，规定当 $k > n$ 时， $\mathrm{A}_n^{k} = \mathrm{C}_n^{k}=0$

$O(n^2)$ 递推

for(int i=0; i<=1000; i++)
{
    C[i][0]=1;
    for(int j=1; j<=i; j++)
        C[i][j]=(C[i-1][j]+C[i-1][j-1])%p;
}

$O(p)$ Lucas定理

详见 [OI数学总结-数论]

P3807 【模板】卢卡斯定理/Lucas 定理

#include <iostream>
#include <string>
#include <vector>
#include <algorithm>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <iomanip>
#include <random>
using namespace std;
#define int long long
#define INF 0x3f3f3f3f3f3f3f3f
#define N (int)(2e5+15)
int Q,n,m,p;
int fac[N];
int qpow(int a,int b,int p)
{
	int ans=1,base=a;
	while(b)
	{
		if(b&1)ans=ans*base%p;
		base=base*base%p;
		b>>=1;
	}
	return ans;
}
// n 选 m
int C(int n,int m,int p)
{
	if(n<m)return 0;
	return fac[n]*qpow(fac[m],p-2,p)%p*qpow(fac[n-m],p-2,p)%p;
}
int lucas(int n,int m,int p)
{
	if(!m)return 1;
	return lucas(n/p,m/p,p)*C(n%p,m%p,p)%p;
}
signed main()
{
	cin >> Q;
	while(Q--)
	{
		cin >> n >> m >> p;
		fac[1]=fac[0]=1;
		for(int i=2; i<=n+m; i++)
			fac[i]=fac[i-1]*i%p;
		cout << lucas(n+m,m,p) << '\n'; // C(n+m,n)
	}
	return 0;
}

快速幂

时间复杂度 $O(\log n)$

空间复杂度 $O(1)$

// 2024年05月28日 10:48:24
int qpow(int a, int b)
{
    int r = 1;
    for(; b; b >>= 1, a = a * a % mod)
        if(b & 1) r = r * a % mod;
    return r;
}

矩阵快速幂

注意点：

$O(n^3\log M)$
初始状态
临时变量不可开过大

新的板子(upd.20220730) P1962 斐波那契数列，此时 $n=2,~M$ 则为代码中的 $n$ 。

#include <iostream>
#include <string>
#include <vector>
#include <algorithm>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <iomanip>
#include <random>
using namespace std;
#define int long long
#define INF 0x3f3f3f3f3f3f3f3f
typedef vector< vector<int> > mat;
#define N (int)()
const int mod=1e9+7;
void add(int &x,int y){x+=y; if(x>=mod)x-=mod;}
// (n*m) x (m*p) = (n*p)
mat mul(mat &A,mat &B)
{
    int n=A.size(),m=A[0].size(),p=B[0].size();
    mat C(n,vector<int>(p));
    for(int k=0; k<m; k++)
        for(int i=0; i<n; i++)
            for(int j=0; j<p; j++)
                add(C[i][j],A[i][k]*B[k][j]%mod);
    return C;
}
mat qpow(mat A,int k)
{
    int n=max(A.size(),A[0].size()); 
    mat B(n,vector<int>(n));
    for(int i=0; i<n; i++) B[i][i]=1;
    while(k)
    {
        if(k&1) B=mul(B,A);
        A=mul(A,A); k>>=1;
    }
    return B;
}
void print(mat &A)
{
    int n=A.size(),m=A[0].size();
    for(int i=0; i<n; i++)
        for(int j=0; j<m; j++)
            cout << A[i][j] << " \n"[j==m-1];
}
signed main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);cout.tie(0);
    // freopen("check.in","r",stdin);
    // freopen("check.out","w",stdout);
    int n; cin >> n;
    if(n<=2) cout << "1";
    else
    {
        mat A(2,vector<int>(2)),B(2,vector<int>(1));
        A[0][0]=A[0][1]=A[1][0]=1;
        B[0][0]=B[1][0]=1;
        A=qpow(A,n-2); A=mul(A,B);
        cout << A[0][0] << '\n';
    }
    return 0;
}

旧的板子

const int p=2009;
int n,t;
struct mat
{
    int n,g[205][205];
    void clear()
    {
        memset(g,0,sizeof(g));
    }
    mat operator*(const mat &o)const
    {
        static mat tmp;
        tmp.n=n;tmp.clear();
        for(int i=1; i<=n; i++)
            for(int k=1; k<=n; k++)
            {
                int r=g[i][k]%p;
                for(int j=1; j<=n; j++)
                    tmp.g[i][j]=(tmp.g[i][j]+r*o.g[k][j]%p)%p;
            }
        return tmp;
    }
}M;

mat qpow(mat a,int b)
{
    static mat ans,base=a;
    ans.n=a.n;ans.clear();
    for(int i=1; i<=ans.n; i++)
        ans.g[i][i]=1;
    while(b)
    {
        if(b&1)ans=ans*base;
        base=base*base;
        b>>=1;
    }
    return ans;
}

判断素数

判断素数朴素版

时间复杂度 $O(Q\sqrt{n})$

#include <bits/stdc++.h>
using namespace std;
#define int long long
int x;
bool ck(int x)
{
    if(x<2)return 0;
    for(int i=2; i<=x/i; i++)
    {
        if(x%i==0)
            return 0;
    }
    return 1;
}
signed main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);cout.tie(0);
    while(cin >> x)
        cout << (ck(x)?'Y':'N') << endl;
    return 0;
}

MillerRabin

$O(Q\log a_i)$

#include <bits/stdc++.h>
using namespace std;
#define int long long
#define test_time 10
mt19937 rd(time(0));
int mul(int a,int b,int p)
{
    int res=(__int128)a*b%p;
    return res;
}
int qpow(int a,int b,int p)
{
    int ans=1,base=a;
    while(b)
    {
        if(b&1)ans=mul(ans,base,p);
        base=mul(base,base,p);
        b>>=1;
    }
    return ans;
}
bool MillerRabin(int n)
{
    if(n<3||n%2==0)return n==2;
    int a=n-1,b=0;
    while(a%2==0)a>>=1,++b;
    for(int i=1,j; i<=test_time; i++)
    {
        int x=rd()%(n-2)+2;
        int v=qpow(x,a,n);
        if(v==1)continue;
        for(j=0; j<b; j++)
        {
            if(v==n-1)break;
            v=mul(v,v,n);
        }
        if(j>=b)return 0;
    } 
    return 1;
}
signed main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);cout.tie(0);
    int x;
    while(cin >> x)
        cout << (MillerRabin(x)?'Y':'N') << endl;
    return 0;
}

线性筛

时间复杂度近似$O(n)$

空间复杂度 $O(n)$

#include <bits/stdc++.h>
using namespace std;
#define int long long
#define gc() getchar()
#define pc(a) putchar(a)
#define INF 0x3f3f3f3f3f3f3f3f
#define MAXN (int)(1e8+5)
#define MAXM (int)(2e6+5)
template<typename T>inline void read(T &k)
{char ch=gc();T x=0,f=1;while(!isdigit(ch)){if(ch=='-')f=-1;ch=gc();}
while(isdigit(ch)){x=(x<<1)+(x<<3)+(ch^48);ch=gc();}k=x*f;}
template<typename T>inline void write(T k)
{if(k<0){k=-k;pc('-');}if(k>9)write(k/10);pc(k%10+'0');}
int n,p;
bool ck[MAXN];
int a[MAXM],tot;
void prime()
{
	for(int i=2; i<=n+2; i++)
	{
		if(!ck[i])a[++tot]=i;
		for(int j=1; j<=tot&&a[j]*i<=n+2; j++)
		{
			ck[a[j]*i]=1;
			if(i%a[j]==0)break;
		}
	}
}
signed main()
{
	read(n);read(p);
	prime();
	for(int i=1,k; i<=p; i++)
	{
		read(k);
		write(a[k]);pc('\n');
	}
	return 0;
}

求欧拉函数

暴力 $O(Q\sqrt{n})$ poj2407

int Euler_phi(int n)
{
    int ans=n;
    for(int i=2; i<=n/i; 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;
}
signed main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);cout.tie(0);
    // freopen("check.in","r",stdin);
    // freopen("check.out","w",stdout);
    int n;
    while(cin>>n&&n!=0)
        cout << Euler_phi(n) << endl;

    return 0;
}

线性筛求欧拉函数 $O(n)$

~~poj2478 （这个题要改成 $\varphi$ 的前缀和）~~

int phi[N],prime[N],pcnt;
bool ck[N];
void Euler()
{
    ck[1]=1;
    phi[1]=1;
    for(int i=2; i<=N-5; i++)
    {
        if(!ck[i])
        {
            prime[++pcnt]=i;
            phi[i]=i-1;
        }
        for(int j=1; j<=pcnt&&i*prime[j]<=N-5; j++)
        {
            int pos=i*prime[j];
            ck[pos]=1;
            if(i%prime[j])
            {
                phi[pos]=phi[i]*phi[prime[j]];
            }
            else
            {
                phi[pos]=phi[i]*prime[j];
                break;
            }
        }
    }
}

杜教筛

参考 P4213 【模板】杜教筛

$T(n)=\sum_{i=1}^{\left\lfloor\frac{n}{m}\right\rfloor} \sqrt{\frac{n}{i}}=\mathcal{O}\left(\frac{n}{\sqrt{m}}\right)$

当 $m = n^{\frac{2}{3}}$ 时，$T(n) = \mathcal{O}(n^{\frac{2}{3}})$ 。

// 2024年06月26日 00:32:57
#include <bits/stdc++.h>
using namespace std;
#define int long long
#define INF 0x3f3f3f3f3f3f3f3f
void up(int &x, int y) { x < y ? x = y : 0; }
void down(int &x, int y) { x > y ? x = y : 0; }
#define rep(i, a, b) for(int i = (a), i##END = (b); i <= i##END; i++)
#define Rep(i, a, b) for(int i = (a), i##END = (b); i >= i##END; i--)
#define N ((int)(5e6 + 15))

char ck[N]; int pcnt, p[N]; int phi[N], mu[N];
unordered_map<int, int> sum_phi, sum_mu;
void init(const int _ = N - 5)
{
    phi[1] = mu[1] = 1;
    rep(i, 2, _)
    {
        if(!ck[i]) { p[++pcnt] = i; phi[i] = i - 1; mu[i] = -1; }
        for(int j = 1; j <= pcnt && i * p[j] <= _; j++)
        {
            const int pos = i * p[j]; ck[pos] = true;
            if(i % p[j] == 0) { phi[pos] = phi[i] * p[j]; break; }
            else { phi[pos] = phi[i] * phi[p[j]]; mu[pos] = -mu[i]; }
        }
    }
    rep(i, 1, _) { phi[i] += phi[i - 1]; mu[i] += mu[i - 1]; }
}
int get_sum_phi(int n)
{
    if(n <= N - 5) return phi[n];
    if(sum_phi[n]) return sum_phi[n];
    int res = 1ll * n * (n + 1) / 2;
    for(int l = 2, r; l <= n; l = r + 1)
    {
        r = n / (n / l);
        res -= 1ll * (r - l + 1) * get_sum_phi(n / l);
    }
    return sum_phi[n] = res;
}
int get_sum_mu(int n)
{
    if(n <= N - 5) return mu[n];
    if(sum_mu[n]) return sum_mu[n];
    int res = 1ll;
    for(int l = 2, r; l <= n; l = r + 1)
    {
        r = n / (n / l);
        res -= 1ll * (r - l + 1) * get_sum_mu(n / l);
    }
    return sum_mu[n] = res;
}
signed main()
{
    ios::sync_with_stdio(0);
    cin.tie(0); cout.tie(0);
    // freopen("check.in","r",stdin);
    // freopen("check.out","w",stdout);
    int qwq; cin >> qwq; init();
    while(qwq--) { 
        static int x; cin >> x;
        cout << get_sum_phi(x) << ' ' << get_sum_mu(x) << '\n';
    }
    return 0;
}

BSGS

传送门：洛谷P3846 [TJOI2007] 可爱的质数/【模板】BSGS 题解

求解 $a^x\equiv b\pmod{p}$ ，其中 $a,b,p$ 均给出， $p$ 为质数，$a,b < p$ 。

// 2024年02月21日 17:33:40
#include <bits/stdc++.h>
using namespace std;
#define int long long
#define INF 0x3f3f3f3f3f3f3f3f
void up(int &x,int y) { x < y ? x = y : 0; }
void down(int &x,int y) { x > y ? x = y : 0; }
#define N ((int)())

int P;
int qpow(int a, int b)
{
    int r = 1;
    while(b)
    {
        if(b & 1) r = r * a % P;
        b >>= 1; a = a * a % P;
    }
    return r;
}
map<int,int> mp;
int BSGS(int a, int b)
{
    int t = sqrt(P) + 1;
    for(int i = 0; i < t; i++) mp[b * qpow(a, i) % P] = i;
    a = qpow(a, t);
    if(!a) return (b == 0 ? 1 : -1);
    for(int i = 1; i <= t; i++)
    {
        int val = qpow(a, i);
        int j = mp.find(val) == mp.end() ? -1 : mp[val];
        if(j >= 0 && i * t - j >= 0) return i * t - j;
    }
    return -1;
}
signed main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);cout.tie(0);
    // freopen("check.in","r",stdin);
    // freopen("check.out","w",stdout);
    int a,b; cin >> P >> a >> b;
    int res = BSGS(a,b);
    if(res == -1) cout << "no solution\n";
    else cout << res << '\n';
    return 0;
}

二次剩余

#include <bits/stdc++.h>
using namespace std;
#define int long long
#define INF 0x3f3f3f3f3f3f3f3f
int Q,n,p;
int i2;
struct num
{
    int r,i;
};
num operator*(num a,num b)
{
    return {((a.r*b.r%p+i2*a.i%p*b.i%p)+p)%p,
            ((a.i*b.r%p+a.r*b.i%p)+p)%p};
}
bool operator==(num a,num b)
{
    return a.r==b.r&&a.i==b.i;
}
num qpow(num a,int b)
{
    num ans={1,0},base=a;
    while(b)
    {
        if(b&1)ans=ans*base;
        base=base*base;
        b>>=1;
    }
    return ans;
}
bool ck(int x)
{
    return qpow({x,0},(p-1)>>1)==(num){1,0};
}
mt19937 rd(time(0));
void solve(int &x1,int &x2)
{
    int a=rd()%p;
    while(!a||ck((a*a+p-n)%p))
        a=rd()%p;
    i2=(a*a+p-n)%p;
    x1=qpow((num){a,1},(p+1)>>1).r;
    x2=p-x1;
}
signed main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);cout.tie(0);
    // freopen("check.in","r",stdin);
    // freopen("check.out","w",stdout);
    int x1,x2;
    cin >> Q;
    while(Q--)
    {
        cin >> n >> p;
        if(qpow((num){n,0},(p-1)>>1).r==p-1)
            cout << "Hola!" << endl;
        else if(!n)cout << 0 << endl;
        else
        {
            solve(x1,x2);
            if(x1>x2)swap(x1,x2);
            if(x1==x2)cout << x1 << endl;
            else cout << x1 << " " << x2 << endl;
        }
    }
    return 0;
}

数论分块

具体参考数论分块学习笔记

例题1 洛谷P2261 [CQOI2007]余数求和

例题2 P2424 约数和

例题1代码：

#include <bits/stdc++.h>
using namespace std;
#define int long long
int n,k;
signed main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);cout.tie(0);
    cin >> n >> k;int res=n*k;
    for(int l=1,r; l<=min(k,n); l=r+1)
    {
        r=min(k/(k/l),n);
        res-=(k/l)*(r-l+1)*(l+r)/2;
    }
    cout << res << endl;
    return 0;
}

康托展开

时间复杂度 $O(n\log n)$

#include <bits/stdc++.h>
using namespace std;
#define int long long
#define INF 0x3f3f3f3f3f3f3f3f
#define N (int)(1e6+15)
#define gc() readchar()
#define pc(a) putchar(a)
#define SIZ (int)(1e6+15)
char buf1[SIZ],*p1=buf1,*p2=buf1;
char readchar()
{
    if(p1==p2)p1=buf1,p2=buf1+fread(buf1,1,SIZ,stdin);
    return p1==p2?EOF:*p1++;
}
template<typename T>void read(T &k)
{
    char ch=gc();T x=0,f=1;
    while(!isdigit(ch)){if(ch=='-')f=-1;ch=gc();}
    while(isdigit(ch)){x=(x<<1)+(x<<3)+(ch^48);ch=gc();};
    k=x*f;
}
template<typename T>void write(T k)
{
    if(k<0){k=-k;pc('-');}
    static T stk[66];T top=0;
    do{stk[top++]=k%10,k/=10;}while(k);
    while(top){pc(stk[--top]+'0');}
}
const int mod = 998244353;
int n,a[N],fac[N],tree[N];
#define lowbit(x) (x&-x);
void add(int x,int v)
{
    while(x&&x<=n)
    {
        tree[x]+=v;
        x+=lowbit(x);
    }
}
int sum(int x)
{
    int res=0;
    while(x>=1)
    {
        res+=tree[x];
        x-=lowbit(x);
    }
    return res;
}
signed main()
{
    read(n);
    fac[0]=1;
    for(int i=1; i<=n; i++)
        fac[i]=fac[i-1]*i%mod,add(i,1);
    int ans=1;
    for(int i=1; i<=n; i++)
    {
        read(a[i]);
        ans+=sum(a[i]-1)*fac[n-i]%mod;
        ans%=mod;add(a[i],-1);
    }
    write(ans);pc('\n');
    return 0;
}

逆康托展开

理论复杂度是 $O(n\log n)$ 的，但是实际上用暴力就行了（n!那么大

这个纯属瞎写的。根本用不到，但是复杂度确实是对的

#include <bits/stdc++.h>
using namespace std;
#define int long long
#define INF 0x3f3f3f3f3f3f3f3f
#define N (int)(1e3+15)
int n,a[N],fac[N],Q;
namespace BT
{
    struct node
    {
        int ch[2],num,cnt,s,sz,sd;
    }t[N];
    #define ls(x) t[x].ch[0]
    #define rs(x) t[x].ch[1]
    double alpha=0.7;
    int tot,rt;
    int tmp[N];
    int New(int x)
    {
        t[++tot]={0,0,x,1,1,1,1};
        return tot;
    }
    void push_up(int at)
    {
        t[at].s=t[ls(at)].s+t[rs(at)].s+1;
        t[at].sz=t[ls(at)].sz+t[rs(at)].sz+t[at].cnt;
        t[at].sd=t[ls(at)].sd+t[rs(at)].sd+(t[at].cnt!=0);
    }
    bool CanR(int at)
    {
        double x=t[at].s*alpha;
        return (t[at].cnt)&&(x<=(double)max(t[ls(at)].s,t[rs(at)].s)||
        (double)t[at].sd<=x);
    }
    void CanR_flatten(int &idx,int at)
    {
        if(!at)return;
        CanR_flatten(idx,ls(at));
        if(t[at].cnt)tmp[++idx]=at;
        CanR_flatten(idx,rs(at));
    }
    int CanR_build(int l,int r)
    {
        if(l>=r)return 0;
        int mid=(l+r)>>1;
        int &at=tmp[mid];
        ls(at)=CanR_build(l,mid);
        rs(at)=CanR_build(mid+1,r);
        push_up(at);
        return at;
    }
    void rebuild(int &at)
    {
        int idx=0;
        CanR_flatten(idx,at);
        at=CanR_build(1,idx+1);
    }
    void insert(int x,int &at)
    {
        if(!at){at=New(x);return;}
        if(t[at].num==x)++t[at].cnt;
        else if(x<t[at].num)insert(x,ls(at));
        else insert(x,rs(at));
        push_up(at);
        if(CanR(at))rebuild(at);
    }
    void remove(int x,int &at)
    {
        if(!at)return;
        if(t[at].num==x)
        {
            if(t[at].cnt)--t[at].cnt;
        }else if(x<t[at].num)remove(x,ls(at));
        else remove(x,rs(at));
        push_up(at);
        if(CanR(at))rebuild(at);
    }
    int getval(int x,int at)
    {
        if(!at)return INF;
        if(x<=t[ls(at)].sz)return getval(x,ls(at));
        if(x<=t[ls(at)].sz+t[at].cnt)return t[at].num;
        else return getval(x-t[ls(at)].sz-t[at].cnt,rs(at));
    }
    int uprbd(int x,int at)
    {
        if(!at)return 1;
        if(x==t[at].num&&t[at].cnt)return t[ls(at)].sz+t[at].cnt+1;
        else if(x<t[at].num)return uprbd(x,ls(at));
        else return uprbd(x,rs(at))+t[ls(at)].sz+t[at].cnt;
    }
    int uprbd_gt(int x,int at)
    {
        if(!at)return 0;
        if(x==t[at].num&&t[at].cnt)return t[ls(at)].sz;
        else if(x<t[at].num)return uprbd_gt(x,ls(at));
        else return uprbd_gt(x,rs(at))+t[ls(at)].sz+t[at].cnt;
    }
    int getpre(int x,int at){return getval(uprbd_gt(x,at),at);}
    int getnext(int x,int at){return getval(uprbd(x,at),at);}
}
namespace BIT
{
    int tree[N];
    #define lowbit(x) (x&-x)
    void add(int x,int v)
    {
        while(x&&x<=n)
        {
            tree[x]+=v;
            x+=lowbit(x);
        }
    }
    int sum(int x)
    {
        int res=0;
        while(x>=1)
        {
            res+=tree[x];
            x-=lowbit(x);
        }
        return res;
    }
}
signed main()
{
    using namespace BIT;
    using namespace BT;
    ios::sync_with_stdio(0);
    cin.tie(0);cout.tie(0);
    cin >> n >> Q;fac[0]=1;
    for(int i=1; i<=n; i++)
        fac[i]=fac[i-1]*i,add(i,1);
    while(Q--)
    {
        char op;
        cin >> op;
        if(op=='P')
        {
            int k,r;cin >> k;
            tot=rt=0;--k;
            for(int i=1; i<=n; i++)
                insert(i,rt);
            for(int i=1; i<=n; i++)
            {
                r=k%fac[n-i];k/=fac[n-i];
                int p=getval(k+1,rt);
                // cout << p << endl;
                cout << p << " \n"[i==n];
                remove(p,rt);k=r;
            }
        }else
        {
            int ans=1;
            memset(tree,0,(n+1)*sizeof(int));
            for(int i=1; i<=n; i++) add(i,1);
            for(int i=1,x; i<=n; i++)
            {
                cin >> x;
                ans+=sum(x-1)*fac[n-i];
                add(x,-1);
            }
            cout << ans << endl;
        }
    }
    return 0;
}

高斯消元

P3389 【模板】高斯消元法

#include <bits/stdc++.h>
using namespace std;
#define int long long
#define INF 0x3f3f3f3f3f3f3f3f
void up(int &x,int y) { x < y ? x = y : 0; }
void down(int &x,int y) { x > y ? x = y : 0; }
#define N ((int)(150))

const double eps = 1e-10;
int dcmp(int x)
{
    if(abs(x) <= eps) return 0;
    return x > eps ? 1 : -1;
}
int n; double a[N][N];
signed main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);cout.tie(0);
    // freopen("check.in","r",stdin);
    // freopen("check.out","w",stdout);
    cout << fixed << setprecision(2);
    cin >> n;
    for(int i=1; i<=n; i++)
        for(int j=1; j<=n+1; j++)
            cin >> a[i][j];
    for(int i=1; i<=n; i++) // col
    {
        int p = i;
        for(int j=i+1; j<=n; j++)
            if(abs(a[j][i]) > abs(a[p][i])) p=j; // row
        for(int j=1; j<=n+1; j++)
            swap(a[i][j],a[p][j]); // swap p to i for proc
        if(!dcmp(a[i][i])) return cout << "No Solution\n",0;
        for(int j=1; j<=n; j++)
            if(j != i)
            {
                double tmp = a[j][i] / a[i][i];
                for(int k=1; k<=n+1; k++)
                    a[j][k] -= a[i][k] * tmp;
            }
    }
    for(int i=1; i<=n; i++)
        cout << a[i][n+1] / a[i][i] << '\n';
    return 0;
}

q779

https://q779.cn/2022/06/06/oi-mo-ban-shu-xue/

本博客所有文章除特别声明外，均采用 CC BY-NC-ND 4.0 许可协议。转载请注明来源 q779 !

数学总结

洛谷P4037 [JSOI2008]魔兽地图题解

2022-06-07 OI

算法 DP 数据结构

洛谷P1272 重建道路题解

2022-06-06 OI

算法图论 DP 数据结构

OI模板-数学

OI模板-数学

最大公约数 gcd

扩展欧几里德算法(Exgcd)

逆元

第一种方法：$O(n)$ 预处理 ，$O(1)$ 查询

第二种方法：费马小定理 $O(\log p)$ 查询 $x$ 的逆元

第三种方法：$\text{Exgcd}~O(\log p)$ 查询 $x^{-1}$

第四种方法（特殊）：$\mathcal{O}(n)$ 预处理阶乘的逆元

欧拉函数 $\varphi$

第一种方法：$O(\sqrt{n})$ 求 $\varphi(n)$

第二种方法：线性筛 $O(n)$ 求解

扩展欧拉定理

Lucas定理

中国剩余定理

扩展中国剩余定理

莫比乌斯函数 $\mu$

排列数

组合数

$O(n^2)$ 递推

$O(p)$ Lucas定理

快速幂

矩阵快速幂

判断素数

判断素数 朴素版

MillerRabin

线性筛

求欧拉函数

杜教筛

BSGS

二次剩余

数论分块

康托展开

逆康托展开

高斯消元

第一种方法：$O(n)$ 预处理，$O(1)$ 查询

判断素数朴素版