学习笔记-数学

· · 算法·理论

数学

通常以组合数学，数论及杂项等形式出现。

快速幂

用于在 O(\log n) 时间内求 a^n 或是 a^n\bmod p 的解的情况。有分治和倍增两种写法，倍增更优。核心代码如下：

int fast_pow(int x,int p)
{
    int res=1;
    while(p){
        if(p&1) res=res*x%mod;
        x=x*x%mod, p>>=1;
    }
    return res;
}

int quick_pow(int x,int p)
{
    if(!p) return 1;
    int res=quick_pow(x,p>>1);
    res=res*res%mod;
    if(p&1) res=res*x%mod;
    return res;
}

欧几里得算法

欧几里得算法用于在 O(\log \min(a,b)) 复杂度内求解 \gcd(a,b) 的算法，其核心是：

\gcd(a,b)=\gcd(b,a\bmod b)

核心代码如下：

int gcd(int a,int b)
{
    if(b==0) return a;
    return gcd(b,a%b);
}

扩展欧几里德算法（exgcd）

用于求二元线性丢番图方程 ax+by=gcd(a,b) 的一组解。

设

ax_1+by_1=\gcd(a,b) bx_2+(a\bmod b)y_2=\gcd(b,a\bmod b)

可以得到

\begin{cases} x_1=y_2\\ y_1=x_2-\lfloor \frac{a}{b} \rfloor y_2 \end{cases}

递归求解即可。核心代码如下：

void ex_gcd(int a,int b,int &x,int &y)
{
    if(!b) { x=1; y=0; return ; }
    ex_gcd(b,a%b,y,x); y-=a/b*x;
}

用扩展欧几里得算法求出 ax+by=gcd(a,b) 的特解后，求 ax+by=c 的解如下：

  1. 记 d=\gcd (a,b)，ax+by=gcd(a,b) 的一组特解为 x_0,y_0，在 ax_0+bx_0=d 两边同时乘以 \frac{c}{d}，得
\frac{c}{d}\rfloor ax_0+\frac{c}{d} bx_0=c
  1. 得出 ax+by=c 的一组特解 (x_0,y_0) 为
\begin{cases} x_{0}^{'}=\frac{c}{d} x_0 \\ y_{0}^{'}=\frac{c}{d} y_0 \end{cases}

模逆元

对于非零整数a,m，如果存在 b 使得 ab\equiv 1 \mod m，就称 b 是 a 模 m 意义下的逆元。

扩展欧几里得算法

利用扩欧求解 ax+my=gcd(a,m)。得到的 x 模 m 就是 a 的逆元。

快速幂

由费马小定理可得对于 a\perp p 有

a^{p-1}\equiv 1 \mod p

所以 a^{-1} \bmod p 等于 a^{p-2}\bmod p，利用快速幂求解。

递推求一组逆元

对于素数 p，有递推式

i^{-1}\equiv -\lfloor \frac{p}{i} \rfloor (p \bmod i)^{-1} \mod p

模板代码如下： 

#include<bits/stdc++.h>
#define int long long
using namespace std;

int N,P;
int inv[3000005];

signed main()
{
    cin>>N>>P;
    inv[1]=1; cout<<1<<'\n';
    for(int i=2;i<=N;i++){
        inv[i]=(P-P/i)*inv[P%i]%P;
        cout<<inv[i]<<'\n';
    } 
}

矩阵

不会用。

给出矩阵定义，矩阵乘法，快速幂的模板。 

#include<bits/stdc++.h>
#define int long long
using namespace std;

constexpr int Mod=1e9+7;
int N,k;

struct Matrix{
    int mat[105][105];
    Matrix operator * (const Matrix &a) const{
        Matrix b;
        memset(b.mat,0,sizeof b.mat);
        for(int i=1;i<=N;i++) 
            for(int j=1;j<=N;j++) 
                for(int k=1;k<=N;k++)
                    b.mat[i][j]=(b.mat[i][j]+mat[i][k]*a.mat[k][j])%Mod;
        return b;
    }
}A;

Matrix fast_pow(Matrix a,int p){
    Matrix res;
    memset(res.mat,0,sizeof res.mat);
    for(int i=1;i<=N;i++) res.mat[i][i]=1;
    while(p){
        if(p&1) res=res*a;
        a=a*a;
        p>>=1;
    }
    return res;
}

signed main()
{
    cin>>N>>k;
    for(int i=1;i<=N;i++)
        for(int j=1;j<=N;j++)
            cin>>A.mat[i][j];
    Matrix res=fast_pow(A,k);
    for(int i=1;i<=N;i++){
        for(int j=1;j<=N;j++) cout<<res.mat[i][j]<<' ';
        cout<<'\n';
    }
}