与程序竞赛有关的数学知识点

容斥原理

集合的并集

A1,A2...,An 是有限集合,则 
|A1A2An|=ni=1|Ai|ni=1j>i|AiAj|ni=1j>ik>j|AiAjAk|+(1)n|A1A2An|

Sylvester公式

给定集合N和具有性质i的集合 A1,A2...,An ,则 
|A¯¯¯1A¯¯¯2A¯¯¯n|=|N|(ni=1|Ai|ni=1j>i|AiAj|+(1)n|A1A2An|)

[列题] 区间(S,E]中与n互质的元素个数

Ai 为n的第i个质因子 pi 的倍数的集合且 Ai(S,E],i=1,2,3k ,那么 |Ai|=[Epi][Spi],|AiAj|=[Epipj][Spipj] ,然后套用上面公式即可解决问题了。关键是将这个公式表达出来,有没有发现一共有 2n 项多项式,而且组合数量为奇数时负,组合数量为偶数时符号为正。

LL POIAE(int n,int S,int E){
    //质因子分解
    int p[32],u=0;
    for(int i=2;i*i<=n;i++)
        if(n%i==0){
            p[u++]=i;
            while(n%i==0) n/=i;
        }
    if(n>1) p[u++]=n;

    LL qua=0;
    for(int i=1;i<1<<u;i++){
        int cnt=0,com=1;
        for(int k=0;k<u;k++)
            if(i>>k&1) cnt++,com*=p[k];
        qua+=cnt&1?E/com-S/com:-E/com+S/com;
    }
    return E-S-qua;
}
 
 
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19

欧拉函数φ(n)

欧拉函数φ(n) 是小于等于n且与n互素的正整数的个数,假设 n=pα11pα22pα33pαkk ,则有 

φ(n)=n(11p1)(11p2)(11p3)(11pk)

证明:

Ai 为1到n之间 pi 的倍数的集合,i=1,2,3,…,k,则 
φ(n)= |A¯¯¯1A¯¯¯2A¯¯¯k|  
= |N|(ki=1|Ai|ki=1j>i|AiAj|+(1)k|A1A2Ak|)  
= n(ki=1npiki=1j>inpipj+ki=1j>ih>jnpipjph+(1)knp1p2p3pk)  
= n(11p1)(11p2)(11p3)(11pk)

性质:

  • 欧拉函数是积性函数,即是说若m,n互质 φ(mn)=φ(m)φ(n)  。
  • 若n是质数p的k次幂, φ(n)=φ(pk)=pkpk1=(p1)pk1 ,因为除了p的倍数外,其他数都跟n互质。

欧拉函数φ(n) c++实现代码:

我们为了方便根据以上性质将其变形得到

φ(n)=(p11)pα111(p11)pα212...(pk1)pαk1k

如此便可写出φ(n)的 o(n)  的代码了。

LL Eular(int n){
    LL ans=1;
    for(int i=2;i*i<=n;i++)
        if(n%i==0){
            ans*=i-1;
            n/=i;
            while(n%i==0) n/=i,ans*=i;
        }
    if(n>1) ans*=n-1;
    return ans;
}
 
 
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11

欧拉函数值表筛法

int e[MAX_N];
void Euler(){
    for(int i=0;i<MAX_N;i++) e[i]=i;
    for(int i=2;i<MAX_N;i++)
        if(e[i]==i)
            for(int k=i;k<MAX_N;k+=i) e[k]=e[k]/i*(i-1);
}
 
 
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7

[列1]区间[1,n]与n的最大公约数的和

i=1ngcd(i,n)=d|ndφ(n/d)

[列2]

UVA11426 求  ni=1nj=i+1gcd(i,j)

int phi[MAX_N];LL res[MAX_N];
void solve(){//筛法的魅力
    for(int i=1;i<MAX_N;i++) phi[i]=i;
    for(int i=2;i<MAX_N;i++)
        if(phi[i]==i)
            for(int j=i;j<MAX_N;j+=i)
                phi[j]=phi[j]/i*(i-1);
    for(int i=1;i<MAX_N;i++)
        for(int k=i;k<MAX_N;k+=i)
            res[k]+=1LL*i*phi[k/i];
    for(int i=1;i<MAX_N;i++)
        res[i]+=res[i-1]-i;
}
 
 
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13

[列3]

[1,n]与n互素的和 

gcd(i,n)=1i<ni=nφ(n)/2


莫比乌斯函数

莫比乌斯函数 μ(n) ,在狄利克雷卷积的乘法中与恒等函数互为逆元, μ(1)=1 ,对于无平方因子数 n=ti=1pi μ(n)=(1)t ,对于有平方因子数 n μ(n)=0

o(n)  的板子

int mu(int n) {
    int cnt = 0;
    for(int i = 2; i * i <= n; ++i) {
        if(n%i == 0) {
            n /= i;++cnt;
            if(n%i == 0) return 0;
        }
    }
    if(n > 1) ++cnt;
    return cnt&1?-1:1;
}
 
 
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11

同余

两个整数 a,b,若它们除以正整数 m所得的余数相等,则称 a,b对于模 m同余,记作 ab(modm)

其性质这里不赘述了,可以参看 维基


费马小定理

费马小定理是数论中的一个定理:假如a是一个整数,p是一个质数,那么  apa 是p的倍数,可以表示为 

apa(modp)

如果a不是p的倍数,这个定理也可以写成 
ap11(modp)

[列题]计算 2100 除以13的余数 
2100212×8+4(mod13)(212)824(mod13)1816(mod13)16(mod13)3(mod13)  
则答案为3


欧拉定理

gcd(a,n)=1 ,则

aφ(n)1(mod n)


扩展欧几里得算法

给予二整数a、b,必存在有整数x、y使得ax + by = gcd(a,b) 
对于这个问题,扩展欧几里得算法能够快速解出使上式成立(x,y),复杂度o(log(n))

int exGcd(int a,int b,int &x,int &y){
    if (b==0){
        x=1,y=0;
        return a;
    }
    int q=exGcd(b,a%b,y,x);
    y-=a/b*x;
    return q;
}
 
 
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9

通过上面代码可以得到一组特解,记为 (x0,y0) ,那么通解可以表示为

x=x0bky=y0+akkZ

例题 POJ1061


模逆元

ab1(mod m) ,称a和b对于模数m来说互为逆元

定理:a存在模m的乘法逆元的充要条件是gcd(a,m) = 1

证明:

  1. 充分性:由欧拉定理知 aφ(m)1(mod m) ,那么必有 b=aφ(m)1  使 ab1(mod m)
  2. 必要性:既然 ab1(mod m) ,即有 abmk=1=kgcd(a,m),kZ,kZ 成立,即有 gcd(a,m)=1

[例] 已知 a, m, 求b

将定义式改写成

ab1=mk,kZabmk=1,kZ  
利用扩展欧几里得算法不难求得一组解 (b,k)  使得

ab+mk=gcd(a,m)

结合定理得 
gcd(a,m)1 ,那么b无解 
gcd(a,m)=1 ,那么 b=bgcd(a,m)=b

代码

那么也能写出代码了,若无答案返回-1,若有返回最小的正解

int mod_inverse(int a,int m){
    int x,y;
    if(exGcd(a,m,x,y)==1) return (m+x%m)%m;
    return -1
}
 
 
  • 1
  • 2
  • 3
  • 4
  • 5

o(n)得出[1,n]的(mod p)逆元

该算法要求所mod p为素数且n< p,记  a  的逆元表示为  a1  
k=p%i ,假设以求得 k1(modp) ,则

ii1kk11(modp)  
ii1(p%i)(p%i)1(modp)  
ii1(p[p/i]i)(p%i)1(modp)  
ii1((i1)p+p[p/i]i)(p%i)1(modp)  
ii1i(p[p/i])(p%i)1(modp)  
由于p为素数,则 gcd(i,p)=1 ,则 
i1(p[p/i])(p%i)1(modp)

取[0,p)中的解 
i1=(p[p/i])(p%i)1%p

代码为

int inv[MAXN];
inv[1] = 1;
for (int i = 2; i<MAXN; i++)
    inv[i] = inv[MOD%i]*(MOD-MOD/i)%MOD;
 
 
  • 1
  • 2
  • 3
  • 4

孙子定理(中国剩余定理)

用现代数学的语言来说明的话,中国剩余定理给出了以下的一元线性同余方程组: 

(S):xa1(modm1)xa2(modm2)xan(modmn)

有解的判定条件,并用构造法给出了在有解情况下解的具体形式。

中国剩余定理说明:假设整数 m1,m2,...,mn 其中任两数互质,则对任意的整数: a1,...,an ,方程组 (S) 有解,并且通解可以用如下方式构造得到:

  1. 设  M=m1×m2××mn=i=1nmi 是整数 m1,m2,...,mn 的乘积,并设 Mi=M/mi,i{1,2,,n} ,即 Mi 是除了 mi 以外的n − 1个整数的乘积。
  2. ti=M1i mi 的逆元:  tiMi1(modmi),i{1,2,,n}
  3. 方程组 (S) 的通解形式为: x=a1t1M1+a2t2M2++antnMn+kM=kM+ni=1aitiMi,kZ. 在模  M 的意义下,方程组  (S) 只有一个解: 
    x=i=1naitiMi.

证明

对于  j{1,2,,n},ji,gcd(mi,mj)=1

考察乘积 aitiMi 可知:

aitiMiai1ai(modmi),  
j{1,2,,n},ji,aitiMi0(modmj).

所以  x=a1t1M1+a2t2M2++antnMn  满足:

i{1,2,,n},x=aitiMi+jiajtjMjai+ji0ai(modmi).

到此以得证 (S) 的一个解 x=ni=1aitiMi.

通解

假设  x1 x2 都是方程组  (S) 的解,那么:

i{1,2,,n},x1x20(modmi). ,即  mi|x1x2

即  ni=1mi|x1x2 ,即  M|x1x2

到此说明方程组  (S)  的任何两个解之间必然相差 M 的整数倍,则通解可表示为 

{kM+i=1naitiMi;kZ}.

代码

LL CRT(int *m,int *a,int n){
    LL M=1,ans=0,x,t_i;
    for (int i=1;i<=n;i++) M*=m[i];
    for (int i=1;i<=n;i++){
        LL M_i=M/m[i];
        exGcd(m[i],M_i,x,t_i);     //求 M_i 模 m[i] 的逆元 t_i
        ans=(ans+a[i]*t_i*M_i)%M;
    }
    return ans;
}
 
 
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10

矩阵

数学上,一个m×n的矩阵是一个由m行n列元素排列成的矩形阵列。矩阵里的元素可以是数字、符号或数学式。

大小相同(行数列数都相同)的矩阵之间可以相互加减,具体是对每个位置上的元素做加减法。矩阵的乘法则较为复杂。两个矩阵可以相乘,当且仅当第一个矩阵的列数等于第二个矩阵的行数。矩阵的乘法满足结合律和分配律,但不满足交换律。

矩阵乘法

矩阵相乘最重要的方法是一般矩阵乘积。它只有在第一个矩阵的列数(column)和第二个矩阵的行数(row)相同时才有定义。一般单指矩阵乘积时,指的便是一般矩阵乘积。若A为  m×n 矩阵,B为 n×p 矩阵,则他们的乘积AB(有时记做A · B)会是一个  m×p 矩阵。其乘积矩阵的元素如下面式子得出: 

(AB)ij=r=1nairbrj=ai1b1j+ai2b2j++ainbnj.

矩阵快速幂

类比快速幂不难的到矩阵快速幂,若有矩阵a,则求  an  代码如下

int N;
struct matrix{LL m[MAX_N][MAX_N];};
matrix multi(matrix a,matrix b){
    matrix tmp;
    for(int i=1;i<=N;i++)
        for(int j=1;j<=N;j++){
            tmp.m[i][j]=0;
            for(int k=1;k<=N;k++)
                tmp.m[i][j]+=a.m[i][k]*b.m[k][j];
            tmp.m[i][j]%=MOD;
        }
    return tmp;
}
matrix fast_mod(matrix a,int n){
    matrix res;
    for(int i=1;i<=N;i++)
        for(int j=1;j<=N;j++)
            res.m[i][j]=(i==j);
    while(n){
        if(n&1) res=multi(res,a);
        a=multi(a,a);
        n>>=1;
    }
    return res;
}
 
 
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25

线性递推关系与矩阵乘法

设数列  {hn}  满足  k  阶常系数线性递推关系: 

hn=C1hn1+C2hn2+C3hn3++Ckhnk+bn

bn  可以为常数,也可以是关于n的函数

若  bn  为常数,则构造转移矩阵 

M=C110000C201000C300100Ck100010Ck00000bn00001(k+1)×(k+1)

与初始向量 
X=hk1hk2hk3h2h1h01(k+1)×1

易见 
Y=MX=C110000C201000C300100Ck100010Ck00000bn00001hk1hk2hk3h2h1h01=hkhk1hk2h3h2h11

事实上我们可以发现,对于任意的 k 阶常系数线性递推关系,我们总可以构造一个 k×k 或 (k+1)×(k+1) 的转移矩阵 M, 对于初始值向量 
X=hk1hk2hk3h2h1h0k×1hk1hk2hk3h1h0bn(k+1)×1

使得  Y=Mnk+1XY 第一行第一列的元素恰好为  hn

利用矩阵乘法计算递推数列的某一项代码

使用上面的两个知识点就不难敲出代码了

#include<stdio.h>
#include<string.h>
#include<queue>
#include<vector>
#define INF 0x3f3f3f3f
#define MAX_N 10
#define MOD 1000000007
using namespace std;
typedef long long LL;
int N;LL b_n=0,C[MAX_N],h[MAX_N];
struct matrix{LL m[MAX_N][MAX_N];};
matrix multi(matrix a,matrix b){
    matrix tmp;
    for(int i=1;i<=N;i++)
        for(int j=1;j<=N;j++){
            tmp.m[i][j]=0;
            for(int k=1;k<=N;k++)
                tmp.m[i][j]+=a.m[i][k]*b.m[k][j];
            tmp.m[i][j]%=MOD;
        }
    return tmp;
}
matrix fast_mod(matrix a,int n){
    matrix res;
    for(int i=1;i<=N;i++)
        for(int j=1;j<=N;j++)
            res.m[i][j]=(i==j);
    while(n){
        if(n&1) res=multi(res,a);
        a=multi(a,a);
        n>>=1;
    }
    return res;
}
void init(matrix &res,matrix &H){
    for(int i=1;i<=N;i++) res.m[1][i]=C[i];
    res.m[1][N]=b_n;
    for(int i=2;i<=N;i++)
        for(int j=1;j<=N;j++)
            res.m[i][j]=(i==j+1);
    res.m[N][N-1]=0;res.m[N][N]=1;

    for(int i=1;i<=N;i++){
        H.m[i][1]=h[N-i];
        for(int j=2;j<=N;j++) H.m[i][j]=0;
    }
    H.m[N][1]=1;
}
void slove(int k,int n){
    matrix res,H;
    init(res,H);

    res=fast_mod(res,n-k+1);
    res=multi(res,H);
    LL ans=res.m[1][1];

    printf("%lld\n",ans);
}

int main()
{
    int k,n;
    while(scanf("%d%d",&n,&k)!=EOF){
        for(int i=1;i<=k;i++) scanf("%lld",&C[i]);
        for(int i=1;i<=k;i++) scanf("%lld",&h[i]);
        N=k+1;
        slove(k,n);
    }
    return 0;
}
 
 
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 35
  • 36
  • 37
  • 38
  • 39
  • 40
  • 41
  • 42
  • 43
  • 44
  • 45
  • 46
  • 47
  • 48
  • 49
  • 50
  • 51
  • 52
  • 53
  • 54
  • 55
  • 56
  • 57
  • 58
  • 59
  • 60
  • 61
  • 62
  • 63
  • 64
  • 65
  • 66
  • 67
  • 68
  • 69
  • 70

例题 51NOD 1126

《挑战》上的矩阵快速幂模板

以斐波那契数列第n项为例

#include<stdio.h>
#include<vector>
using namespace std;

typedef long long LL;
typedef vector<LL> vec;
typedef vector<vec> mat;
const LL MOD = 1000000009;

mat mul(mat &A,mat &B){
    mat C(A.size(),vec(B[0].size()));
    for(int i=0;i<A.size();++i){
        for(int k=0;k<B.size();++k){
            for(int j=0;j<B[0].size();++j)
                C[i][j]=(C[i][j]+A[i][k]*B[k][j])%MOD;
        }
    }
    return C;
}

mat pow(mat A,LL n){
    mat B(A.size(),vec(A.size()));
    for(int i=0;i<A.size();++i) B[i][i]=1;
    while(n){
        if(n&1) B=mul(B,A);
        A=mul(A,A);
        n>>=1;
    }
    return B;
}

int main()
{
    mat A(2,vec(2));
    A[0][0]=A[0][1]=A[1][0]=1;
    A[1][1]=0;
    LL n;
    scanf("%lld",&n);
    A=pow(A,n);
    printf("%lld\n",A[1][0]);
    return 0;
}
 
 
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 29
  • 30
  • 31
  • 32
  • 33
  • 34
  • 35
  • 36
  • 37
  • 38
  • 39
  • 40
  • 41
  • 42

其他

约数个数

仍何一个整数n都可以表示成 n=pa11pa22...pakk ,那么,n的约数个数 

d|n1=(a1+1)(a2+1)...(ak+1)

例题  LOJ1341

a的k次幂前n位数

用该方法要注意精度

ak=10lgak=10klga=10[klga]+{klga}n[10n110{klga}]

例题 LOJ1282

定积分

Simpson`s 3/8 rule

const double eps=1e-8;
inline double f(double x){
    return x*x;//被积函数
}
inline double getAppr(double a,double b){//Simpson`s 3/8 rule
    return (b-a)*(f(a)+3*f((2*a+b)/3)+3*f((a+2*b)/3)+f(b))/8.0;
}
double simpson(double l,double r){//积分区间[l,r]
    double sum=getAppr(l,r);
    double mid=(l+r)/2;
    double suml=getAppr(l,mid);
    double sumr=getAppr(mid,r);
    return fabs(sum-suml-sumr)<eps?sum:simpson(l,mid)+simpson(mid,r);
}
 
 
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14

例题 地狱飞龙


线性筛法筛素数、莫比乌斯函数、欧拉函数

线性筛法(欧拉筛法)可以在 O(n) 的时间内获得 [1,n] 的所有素数。算法保证每个合数都会被它的最小质因子筛掉,所以复杂度是线性的。同时,我们可以利用这一特性,结合积性函数的性质,在 O(n)的时间内筛出一些积性函数的值。

bool isNotPrime[MAXN + 1];
int mu[MAXN + 1], phi[MAXN + 1], primes[MAXN + 1], cnt;
inline void euler()
{
    isNotPrime[0] = isNotPrime[1] = true;
    mu[1] = 1;
    phi[1] = 1;
    for (int i = 2; i <= MAXN; i++) {
        if (!isNotPrime[i]) {
            primes[++cnt] = i;
            mu[i] = -1;
            phi[i] = i - 1;
        }
        for (int j = 1; j <= cnt; j++) {
            int t = i * primes[j];
            if (t > MAXN) break;
            isNotPrime[t] = true;
            if (i % primes[j] == 0) {
                mu[t] = 0;
                phi[t] = phi[i] * primes[j];
                break;
            } else {
                mu[t] = -mu[i];
                phi[t] = phi[i] * (primes[j] - 1);
            }
        }
    }
}
 
 
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26
  • 27
  • 28
  • 0
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值