清北学堂五一数学营笔记

Southern_Dynasty · 2023-04-29 16:25:45 · 个人记录

矩阵

https://www.luogu.com.cn/blog/531930/shou-xie-yi-ge-ju-zhen-lei。

P4159 [SCOI2009] 迷路

如果边权是 1 的话，答案就是邻接矩阵的 k 次方的第 1 行第 n 列。

这个题有边权且很小（0\le w\le 10），考虑拆成 1 边权。拆边运算量太大，考虑拆点。

我们把一个点拆成 9 个。记 id(i,j) 表示点 i 拆出来的第 j+1 个点，其中 0\le j\le 8，满足 id(i,j) 到 i 的距离为 j。id(i,0) 就是点 i 本身。

我们把 id 看做一个 n\times 9 的二维数组，这样我们就可以给每一个位置一个编号：id(i,j)=9(i-1)+j。容易发现点的范围是 [0,9n-1]。

建边我们“向内”建，由 id(i,j) 向 id(i,j-1) ；连权值为 1 的边，这样从 id(i,j) 往前走 j 条边就可以到达 id(i,0)。

然后是原图的建边。如果 w(i,j)\not=0，那么由 id(i,0) 向 id(j,w(i,j)-1) 连边。也就是先从 id(i,0) 经过一条边到达 id(j,w(i,j)-1)，然后再经过 w(i,j)-1 条边到达 id(j,0)，一共 w(i,j) 条边，和题目相符。

最后我们就得到了一个大小为 9n\times9n 的邻接矩阵，直接用矩阵快速幂计算该矩阵的 t 次方即可。答案就是该矩阵的第 id(1,0) 行第 id(n,0) 列的数。

由于点的范围是 [0,9n-1] 的，所以矩阵乘法的时候直接从 0 枚举到 9n-1 就可以了。最后不要忘记取模。

#include<bits/stdc++.h>
//#pragma GCC optimize("Ofast")
#define gt getchar
#define pt putchar
#define y1 y233
//typedef unsigned int uint;
typedef long long ll;
typedef unsigned long long ull;
//typedef __int128 lll;
//typedef __uint128_t ulll;
const int mod=2009;
using namespace std;
inline bool __(char ch){return ch>=48&&ch<=57;}
inline int read(){
    int x=0;bool sgn=0;char ch=gt();
    while(!__(ch)&&ch!=EOF){sgn|=(ch=='-');ch=gt();}
    while(__(ch)){x=(x<<1)+(x<<3)+(ch-48);ch=gt();}
    return sgn?-x:x;
}
template<class T>
inline void print(T x){
    static char st[70];short top=0;
    if(x<0)pt('-');
    do{st[++top]=x>=0?(x%10+48):(-(x%10)+48),x/=10;}while(x);
    while(top)pt(st[top--]);
}
template<class T>
inline void printsp(T x){
    static char st[70];short top=0;
    if(x<0)pt('-');
    do{st[++top]=x>=0?(x%10+48):(-(x%10)+48),x/=10;}while(x);
    while(top)pt(st[top--]);pt(32);
}
template<class T>
inline void println(T x){
    static char st[70];short top=0;
    if(x<0)pt('-');
    do{st[++top]=x>=0?(x%10+48):(-(x%10)+48),x/=10;}while(x);
    while(top)pt(st[top--]);pt(10);
}
inline void put_str(string s){
    int siz=s.size();
    for(int i=0;i<siz;++i) pt(s[i]);
    printf("\n");
}
int n,tot,t;
inline int id(int x,int y){
    return (x-1)*9+y;
}
struct Matrix{
    int val[105][105];
    Matrix(){
        memset(val,0,sizeof(val));
    }
    inline void reset(){
        for(int i=0;i<tot;++i) val[i][i]=1;
    }
}tmp;
inline Matrix operator*(const Matrix &a,const Matrix &b){
    Matrix c;
    for(int i=0;i<tot;++i){
        for(int k=0;k<tot;++k){
            for(int j=0;j<tot;++j){
                c.val[i][j]=(c.val[i][j]+a.val[i][k]*b.val[k][j]%mod)%mod; 
            }
        }
    }
    return c;
}
inline Matrix ksm(Matrix a,int b){
    Matrix c;
    c.reset();
    while(b){
        if(b&1)c=c*a;
        a=a*a,b>>=1;
    }
    return c;
}
signed main(){
    n=read(),t=read(),tot=9*n;
    for(int i=1;i<=n;++i){
        for(int j=1;j<=8;++j) tmp.val[id(i,j)][id(i,j-1)]=1;
        for(int j=1;j<=n;++j){
            int x;
            scanf("%1d",&x);
            if(x) tmp.val[id(i,0)][id(j,x-1)]=1;
        }
    }
    tmp=ksm(tmp,t);
    println(tmp.val[id(1,0)][id(n,0)]);
    return 0;
}

一个小语法

如果你在主函数中定义了一个全局变量 k，然后在某个循环中有定义了一次，如：

int k=1;
signed main(){
    for(int k=2;k<=3;++k) println(k);
    return 0;
}

这时候就会输出内部循环中的 k。但是只要改成 ::k 就可以访问到外部的 k。

数论

判断质数

一、试除法

先判断 1,0，然后直接枚举 [2,x] 的数试除即可。O(x)

显而易见的优化是，质因子总是成对出现的，因此只需要枚举较小的那个质因子（不超过 \sqrt{x}）。试除 [2,\sqrt{x}] 的数，O(\sqrt{x})。

同样的，我们只需要试除质数。可以先 O(\sqrt{x}) 筛出 [2,\sqrt{x}] 的质数然后试除即可。根据素数密度就知道复杂度约为 O(\dfrac{\sqrt{x}}{\ln x})。

代码太简单，不写了。

二、Miller - Rabin

https://www.luogu.com.cn/paste/0ffg5egd

分解质因数

一、朴素做法

直接枚举 [2,\sqrt{x}] 的所有数 a 然后除就可以了。如果 a 是合数，那么 x 已经被 a 的某个更小的质因子除过了，所以 x 必然不是 a 的倍数。因此这样找到的 a 必然是质数。最后，由于 x 最多只有一个 >\sqrt{x} 的质因子，所以特判一下即可。时间复杂度 O(\sqrt{x}) 同样可以通过先筛质数来优化复杂度。

二、Pollard's rho 算法

还没过。

UPD：没讲

CF45G - Prime Problem

纯纯的找规律 & 打表题。

首先我们需要知道哥德巴赫猜想，并且知道它在本题的数据范围内已经得到了验证（也可以自己手动验证）。

首先计算出数列的和 sum=\dfrac{n(n+1)}{2}。初始都归在同一组。

• 如果 sum 为质数（试除法），那么直接输出，结束；

• 如果 sum 为偶数，根据哥猜，可以拆成两个质数之和。通过打表可以发现必然有一个不超过 n。考虑每一个 i，验证即可；

• 如果 sum 为奇数，由于 n\ge2，所以必然能拆出一个 3 来，然后 sum 就变成偶数，此时再考虑上一种情况即可。注意不要考虑上 3。

#include<bits/stdc++.h>
//#pragma GCC optimize("Ofast")
#define gt getchar
#define pt putchar
#define y1 y233
//typedef unsigned int uint;
typedef long long ll;
typedef unsigned long long ull;
//typedef __int128 lll;
//typedef __uint128_t ulll;
const int N=6005;
using namespace std;
inline bool __(char ch){return ch>=48&&ch<=57;}
inline int read(){
    int x=0;bool sgn=0;char ch=gt();
    while(!__(ch)&&ch!=EOF){sgn|=(ch=='-');ch=gt();}
    while(__(ch)){x=(x<<1)+(x<<3)+(ch-48);ch=gt();}
    return sgn?-x:x;
}
template<class T>
inline void print(T x){
    static char st[70];short top=0;
    if(x<0)pt('-');
    do{st[++top]=x>=0?(x%10+48):(-(x%10)+48),x/=10;}while(x);
    while(top)pt(st[top--]);
}
template<class T>
inline void printsp(T x){
    static char st[70];short top=0;
    if(x<0)pt('-');
    do{st[++top]=x>=0?(x%10+48):(-(x%10)+48),x/=10;}while(x);
    while(top)pt(st[top--]);pt(32);
}
template<class T>
inline void println(T x){
    static char st[70];short top=0;
    if(x<0)pt('-');
    do{st[++top]=x>=0?(x%10+48):(-(x%10)+48),x/=10;}while(x);
    while(top)pt(st[top--]);pt(10);
}
inline void put_str(string s){
    int siz=s.size();
    for(int i=0;i<siz;++i) pt(s[i]);
    printf("\n");
}
int n,sum,ans[N];
bool vis[N];
inline bool is_prime(int x){
    if(x<2)return 0;
    for(int i=2;i*i<=x;++i) if(x%i==0)return 0;
    return 1;
}
inline void output(){
    for(int i=1;i<=n;++i) printsp(ans[i]);
    printf("\n");
}
signed main(){
    n=read(),sum=n*(n+1)/2;
    for(int i=1;i<=n;++i) ans[i]=1;
    if(is_prime(sum)){
        output();
        return 0;
    }
    if(sum%2==1&&!is_prime(sum-2))sum-=3,ans[3]=3,vis[3]=1;
    for(int i=2;i<=n;++i){
        if(vis[i])continue;
        if(is_prime(i)&&is_prime(sum-i)){
            ans[i]=2;
            break;
        }
    }
    for(int i=1;i<=n;++i) printsp(ans[i]);
    return 0;
}

逆元

若 b\times c \equiv 1 \pmod p，那么就称 c 为 b 在模 p 意义下的乘法逆元，类似于倒数。借助逆元，我们就可以完成模意义下的除法操作了。

那么怎么求逆元呢？

费马小定理：a^{p-1}\equiv1\pmod p，满足 p 为质数，\gcd(a,p)=1。

我们发现 a\times a^{p-2}=a^{p-1}，所以：

然后这个 a^{p-2} 不就是 a 的逆元了吗！

快速幂计算即可，时间复杂度 O(\log p)。

所以，在 p 为质数，a 于 p 互质的情况下，a\div b \bmod c=a\times b^{p-2} \bmod c。

欧拉函数：\phi(p) 表示 [1,p] 中有多少个数与 p 互质。

欧拉定理：a^{\phi(p)}\equiv1\pmod p，满足 \gcd(a,p)=1

这个时候我们会发现如果一个数 p 是素数，那么 \phi(p)=p-1（除 p 外都与它互质），这样就是费马小定理了。所以费马小定理实际上是欧拉定理的特殊情况。

这样 a 的逆元就是 a^{\phi(p)-1}。

注意：逆元存在的充分条件是 \gcd(a,p)=1，其他情况下不一定存在逆元。

最后的问题就是如何算出 \phi(n) 了。

一、暴力

枚举 [1,n-1] 的所有数并和 n 取 \gcd，时间复杂度 O(n\log n)。

二、优化

暴力的复杂度太高了，我们考虑优化。

• 当 n= 一个质数 p 时，显然 \phi(n)=p-1。

• 当 n=p^2 时，考虑用总数减去不互质的。显然只有 p 的倍数不互质，[1,n] 中一共有 p 个数是 p 的倍数，所以 \phi(n)=n-p=p^2-p=p(p-1)。

• 当 n=p^k 时，考虑将 [1,n] 分为若干长度为 p 的段，一共有 \dfrac{n}{p} 段。每段只有最后一个数（p 的倍数）与 n 不互质，其他 p-1 个数都互质，所以 \phi(p)=\dfrac{n}{p}(p-1)。稍微整理一下（把 p 移到 p-1 下面）得 \phi(n)=n(1-\dfrac{1}{p})。

• 加一个质因子：当 n=p_1^kp_2^k 时，还是按照上面的方法分段，容斥一下得 \phi(n)=n-\dfrac{n}{p_1}-\dfrac{n}{p_2}+\dfrac{n}{p_1p_2}。因式分解一下：先提 n，得 \phi(n)=n(1-\dfrac{1}{p_1}-\dfrac{1}{p_2}+\dfrac{1}{p_1p_2})；分组分解，得 \phi(n)=n(1-\dfrac{1}{p_1})(1-\dfrac{1}{p_2})。

• 然后就可以得出规律了：将 n 质因数分解为 \sum_{i=1}^kp_i^{a_i}，那么答案就是 n\prod_{i=1}^k(1-\dfrac{1}{p_i})。

所以直接分解就可以了，一个小技巧：先除后乘防止溢出。在求 lcm 中也有用。

inline int get_phi(int n){
    int phi=n;
    for(int i=2;i*i<=n;++i){
        if(n%i==0){
            phi=phi/i*(i-1);
            while(n%i==0) n/=i;
        }
    }
    if(n!=1) phi=phi/n*(n-1);
    return phi;
}

这样时间复杂度就是 O(\sqrt{n}) 了。你想先筛一下也行。

求 [1,n] 逆元

假设我们现在要求 ax+by=\gcd(a,b)。令 d=\gcd(a,b)， 由于 \gcd(a,b)=\gcd(b,a \bmod b)，所以我们考虑把 b 和 a\bmod b 也写成类似的形式：

现在我们要从这个式子推出 x,y。把 a\bmod b 拆开，得：

拆项，得：

把带 a 的凑到一起，带 b 的凑到一起，得：

这个时候，我们发现这个式子和原来的式子形式一样了！因此：

因此只要递归算出 bx+(a\bmod b)\times y = d 的解即可。

那递归的终止条件呢？当 b=0 时，\gcd(a,b)=a，因此令 x=1,y=0 即可。

代码实现时可以用一个三个元素的结构体维护 x，y 和 \gcd，这样就可以顺便求出 \gcd(a,b) 了。

代码。

用途？

现在有一个问题：x,y,a,b 为整数，那么 ax+by 能表示出的最小正整数是多少？

答案是 \gcd(a,b)，现在考虑证明：

设 \gcd(a,b)=d，那么 a 可以表示为 a'd，b 可以表示为 b'd，代回原式：

提 d，得：

所以得数一定是 \gcd(a,b) 的倍数，所以 \gcd(a,b) 就是最小正整数。这也说明了，只有得数是 \gcd(a,b) 的倍数的时候，方程才有解。这就是裴蜀定理。

现在要求 ax+by=c ，满足 c 是 \gcd(a,b) 的倍数。怎么做？

• 把得到的 x,y 同乘 \dfrac{c}{\gcd(a,b)} 即可。

现在只求出了一组解，求通解怎么办？

• 假设 x_0,y_0 是一开始 exgcd 的出来的一组解，那么通解可以写成 x=x_0+k\times\dfrac{b}{\gcd(a,b)},y=y_0-k\times\dfrac{a}{\gcd(a,b)} 的形式，其中 k 为任意整数。

而且有结论：上述算法取得的一组解 (x,y) 满足 |x|\le|\frac{b}{2\gcd(a,b)}|,|y|\le|\frac{a}{2\gcd(a,b)}|，证明可以考虑归纳。具体见《深进》P298-P299。

P1082 [NOIP2012 提高组] 同余方程

这个题既然保证有解，那么 \gcd(a,b)=1。

做法一：容易发现这个 x 就是 a 的逆元。由于 b 不一定是质数，因此用欧拉定理即可。

做法二：我们改写一下它的形式。既然 ax\bmod b=1，那么 ax 一定可以写成 by+1 的形式。移个项，变为 ax-by=1，那么直接上 exgcd 就行了。

但是这个题要求最小正整数解。考虑通解公式，容易发现 x\equiv x_0\pmod \dfrac{b}{\gcd(a,b)}，由于 \gcd=1，所以 x\equiv x_0\pmod b，容易发现答案就是 x_0\bmod b，所以 (%mod+mod)%mod 防止负数即可。

CRT/EXCRT

用途是解如下同余方程组

如果不要求 a 两两互质就是 EXCRT，否则就是普通 CRT。

对于普通 CRT，有一个小构造：记 M=\prod_{i=1}^{n}a_i，m_i=\dfrac{M}{a_i},t_i\equiv\dfrac{1}{m_i}\pmod {a_i}，那么答案为 \sum_{i=1}^{n}b_it_im_i。

证明：对于第 i 个方程组，对于 j\not=i，m_j 都是 a_i 的倍数，所以 b_jt_jm_j=0；而对于 i，这三项都不是 a_i 的倍数，所以 b_it_im_i=b_i\times1=b_i。

对于 a 不互质的情况，m_i 与 a_i 也不一定互质，也就不一定存在 t_i（逆元），这时就需要用 EXCRT 了。

EXCRT 在这里。

筛法

一、暴力

暴力判断每个数是否是质数即可。试除法可以做到 O(n\sqrt{n})，Miller - Rabin 可以做到 O(n\log n)。

二、埃氏筛

把每个数的倍数筛掉。

not_prime[1]=1;
for(int i=2;i<=n;++i){
    for(int j=i+i;j<=n;j+=i){
        not_prime[j]=1;
    }
}
for(int i=2;i<=n;++i) if(!not_prime[i]) prime[++cnt]=i;

复杂度是多少？\sum_{i=2}^{n}\dfrac{n}{i}。补一个 \dfrac{n}{1}，提出 n，得 n\sum_{i=1}^{n}\dfrac{1}{i}。有调和级数可得约为 O(n\log n)。

优化：只需要筛质数的倍数。

not_prime[1]=1;
for(int i=2;i<=n;++i){
    if(!not_prime[i]){
        for(int j=i+i;j<=n;j+=i){
            not_prime[j]=1;
        }
    }
}
for(int i=2;i<=n;++i) if(!not_prime[i]) prime[++cnt]=i;

不加证明地给出算法的复杂度：约为 O(n\log\log n)。（是 \log 套 \log，不是 \log 乘 \log）。

其实还可以优化：设 j=ki，当 k<i 时，j 已经被 k 筛过了，因此从 j=i^2 开始筛即可。

not_prime[1]=1;
for(int i=2;i<=n;++i){
    if(!not_prime[i]){
        for(int j=i*i;j<=n;j+=i){
            not_prime[j]=1;
        }
    }
}
for(int i=2;i<=n;++i) if(!not_prime[i]) prime[++cnt]=i;

三、线性筛

埃氏筛的缺点在于，每个数可能被它的质因子筛了多次。所以它的时间复杂度必然不可能达到 O(n)。

线性筛的思想就是，让每个数被它的最小质因子筛掉。

过程：维护 bool 数组，判断某个数是否是质数；再维护一张质数表，里面存放 [2,n] 的所有质数（省略 1）。从小到大扫描 i，如果 i 不是质数，存放到素数表里。然后，取出素数表里的数，把该质数 \times i 筛掉（超过 n 就 break）。

const int N=1e6+5;
int n,prime[N],cnt,T;
bool not_prime[N];
inline void do_prime(){
    not_prime[0]=not_prime[1]=1,cnt=0;
    for(int i=2;i<=N-5;++i){
        if(!not_prime[i])prime[++cnt]=i;
        for(int j=1;j<=cnt&&i*prime[j]<=N-5;++j){
            not_prime[i*prime[j]]=1;
        }
    }
}

可惜，这样还是 O(n\log\log n) 的。

但是我们只需要在第 9 行下面添加如下代码：

if(i%prime[j]==0)break;

这样就 O(n) 了。

为什么？

当 i\bmod prime_j=0 时，那么 i 的最小质因子就是 prime_j（否则前面已经 break 掉了）。而如果 j+1，那么被筛掉的数是 prime_{j+1}\times i，由于 prime_{j+1}>prime_j，所以 prime_j 也是 prime_{j+1}\times i 的最小质因子。这样它就不是由它的最小质因子筛掉了，需要 break 掉。 这样就可以保证每个数都是由它的最小质因子筛掉，时间复杂度 O(n)。

代码：

const int N=1e6+5;
int n,prime[N],cnt,T;
bool not_prime[N];
inline void do_prime(){
    not_prime[0]=not_prime[1]=1,cnt=0;
    for(int i=2;i<=N-5;++i){
        if(!not_prime[i])prime[++cnt]=i;
        for(int j=1;j<=cnt&&i*prime[j]<=N-5;++j){
            not_prime[i*prime[j]]=1;
            if(i%prime[j]==0)break;
        }
    }
}

这玩意还能用来筛积性函数。

积性函数

定义：对于一个函数 f，若 f(a)f(b)=f(ab)，且 \gcd(a,b)=1，那么称 f 为积性函数。如果没有 \gcd(a,b)=1 的限制，那么 f 就是完全积性函数。

举个例子：\phi （欧拉函数）是积性函数。

回忆一下公式：设

那么

假设 a,b 互质，所以 a,b 没有公共的质因子，所以

所以

而

所以 \phi(ab)=\phi(a)\phi(b)，所以 \phi 是积性函数。

现在我们要线性求出 [1,n] 中所有数的 \phi 值。

首先，\phi(1)=1。并且当 p 为质数时，\phi(p)=p-1。这样我们就只需要考虑计算合数的欧拉函数值。

由于我们会筛掉 i\times prime_j，所以直接更新 \phi(i\times prime_j)=\phi(i)\times\phi(prime_j) 就可以了……吗？

并不行，因为两者不一定互质。什么时候不互质？i\bmod prime_j=0 时。所以我们重新算一下 \phi 就好了。等于多少呢？等于 prime_j\times\phi(i)。

考虑证明，令 i=\sum_{j=1}^{n}p_j^{a_j}，那么

由于 prime_j 是 i 的质因子，所以乘上它之后只会让它的指数 +1，而不会多出一个新的质因子。所以：

容易发现后面的部分就是 \phi(i)，所以：

这样我们就可以筛出 [1,n] 所有数的欧拉函数值了。

const int N=1e6+5;
int n,prime[N],cnt,T,phi[N];
bool not_prime[N];
inline void do_prime(){
    not_prime[0]=not_prime[1]=phi[1]=1,cnt=0;
    for(int i=2;i<=N-5;++i){
        if(!not_prime[i])prime[++cnt]=i,phi[i]=i-1;
        for(int j=1;j<=cnt&&i*prime[j]<=N-5;++j){
            not_prime[i*prime[j]]=1;
            phi[i*prime[j]]=phi[i]*phi[prime[j]];
            if(i%prime[j]==0){
                phi[i*prime[j]]=prime[j]*phi[i];
                break;
            }
        }
    }
}

BSGS 算法

Baidu Search Google Search 算法

原题链接。

一、暴力

直接枚举 x，并迭代计算 a^x 并检验即可。

时间复杂度为 O(\text{答案})，明显不靠谱。

代码还是放一下吧：

inline int solve(int a,int b,int p){
    // 求解 a^x % p = b
    int val=1;
    for(int x=0;x;++x){
        if(val==b)return x;
        val=1ll*val*a%p;
        if(val==1)return -1;// 出现循环，无解 
    } 
}

为什么循环一定从 1 开始呢？因为费马小定理（p 是质数）。同理，我们还可以知道循环节长度为 p-1。因此，加上无解判断的代码时间复杂度实际上是 O(p) 的。

所以我们改一下：

inline int solve(int a,int b,int p){
    // 求解 a^x % p = b
    int val=1;
    for(int x=0;x<p-1;++x){
        if(val==b)return x;
        val=1ll*val*a%p;
    }
    return -1;
}

这样我们就可以喜提 0 分了。

二、正解

由于循环节长度只有 p-1，所以我们只需要考虑 p-1 个数。

我们将这些数分组，假设每组有 s 个数。

暴力扫描每一个组，复杂度还是一样的。

现在第一组还是暴力，但是我们需要考虑用更快的方法求出后面的组。

很明显，第二组的每个数都是第一组对应数的 a^s 倍。换句话说，把第二组中的数全部除 a^s，就得到了第一组。

我们可以发现，如果第二组存在 b，那么第一组一定存在 \dfrac{b}{a^s}！

换句话说，如果第 i 组存在 b，那么第一组一定存在 \dfrac{b}{a^{s(i-1)}}。

这样我们只需要在第一组里面查找了。怎么查找？哈希或者二分。

容易发现这是一个类似分块的操作。我们需要扫描第一组，用哈希维护，时间复杂度 O(s)；然后需要扫描每一组找答案，时间复杂度 O(\dfrac{p}{s})。根据均值不等式很容易知道 s=\sqrt{p} 最优，这样时间复杂度就是 O(\sqrt{p})。

我代码里偷懒用了 set，还写了快速幂，时间复杂度是 O(\sqrt{p}\log p) 的，很逊。

代码：

#include<bits/stdc++.h>
//#pragma GCC optimize("Ofast")
#define gt getchar
#define pt putchar
#define y1 y233
//typedef unsigned int uint;
typedef long long ll;
typedef unsigned long long ull;
//typedef __int128 lll;
//typedef __uint128_t ulll;
using namespace std;
inline bool __(char ch){return ch>=48&&ch<=57;}
inline ll read(){
    ll x=0;bool sgn=0;char ch=gt();
    while(!__(ch)&&ch!=EOF){sgn|=(ch=='-');ch=gt();}
    while(__(ch)){x=(x<<1)+(x<<3)+(ch-48);ch=gt();}
    return sgn?-x:x;
}
template<class T>
inline void print(T x){
    static char st[70];short top=0;
    if(x<0)pt('-');
    do{st[++top]=x>=0?(x%10+48):(-(x%10)+48),x/=10;}while(x);
    while(top)pt(st[top--]);
}
template<class T>
inline void printsp(T x){
    static char st[70];short top=0;
    if(x<0)pt('-');
    do{st[++top]=x>=0?(x%10+48):(-(x%10)+48),x/=10;}while(x);
    while(top)pt(st[top--]);pt(32);
}
template<class T>
inline void println(T x){
    static char st[70];short top=0;
    if(x<0)pt('-');
    do{st[++top]=x>=0?(x%10+48):(-(x%10)+48),x/=10;}while(x);
    while(top)pt(st[top--]);pt(10);
}
inline void put_str(string s){
    int siz=s.size();
    for(int i=0;i<siz;++i) pt(s[i]);
    printf("\n");
}
ll a,b,p;
inline ll ksm(ll a,ll b,ll p){
    ll res=1;
    a%=p;
    while(b){
        if(b&1)res=res*a%p;
        a=a*a%p,b>>=1;
    }
    return res;
}
inline ll solve(ll a,ll b,ll p){
    // 求解 a^x % p = b
    ll s=sqrt(p),val=1;
    set<ll>S;
    for(int i=0;i<s;++i){
        S.insert(val);
        val=val*a%p;
    }
    for(int i=0;1ll*i*s<=p;++i){
        // 查找第 i 组（从 0 开始，这样就查找 b/a^si 即可） 
        ll inv=ksm(a,i*s,p);
        ll tmp=b*ksm(inv,p-2,p)%p;
        if(S.count(tmp)){
            // inv 就是第 i 行的第一个数 
            for(int j=i*s;;++j){
                // 暴力查找第 i 组中满足条件的数 
                if(inv==b)return j;
                inv=inv*a%p;
            }
        }
    }
    return -1;
}
signed main(){
    p=read(),a=read(),b=read();
    int ans=solve(a,b,p);
    if(ans==-1)printf("no solution\n");
    else println(ans);
    return 0;
}

容易发现，每次的循环的 tmp 只是在上一次的基础上乘上 \dfrac{1}{a^s}。所以快速幂预处理出 \dfrac{1}{a^s} 就可以去掉那只 \log。

当 i=0 的时候情况稍微有点复杂，特判一下即可。

提供一个哈希表板子。

由于 Hash 表的时间复杂度是期望线性的，所以时间复杂度 O(\sqrt{p}+\log p)=O(\sqrt{p})。

代码：

#include<bits/stdc++.h>
//#pragma GCC optimize("Ofast")
#define gt getchar
#define pt putchar
#define y1 y233
//typedef unsigned int uint;
typedef long long ll;
typedef unsigned long long ull;
//typedef __int128 lll;
//typedef __uint128_t ulll;
const int N=1e5+5;
using namespace std;
inline bool __(char ch){return ch>=48&&ch<=57;}
inline ll read(){
    ll x=0;bool sgn=0;char ch=gt();
    while(!__(ch)&&ch!=EOF){sgn|=(ch=='-');ch=gt();}
    while(__(ch)){x=(x<<1)+(x<<3)+(ch-48);ch=gt();}
    return sgn?-x:x;
}
template<class T>
inline void print(T x){
    static char st[70];short top=0;
    if(x<0)pt('-');
    do{st[++top]=x>=0?(x%10+48):(-(x%10)+48),x/=10;}while(x);
    while(top)pt(st[top--]);
}
template<class T>
inline void printsp(T x){
    static char st[70];short top=0;
    if(x<0)pt('-');
    do{st[++top]=x>=0?(x%10+48):(-(x%10)+48),x/=10;}while(x);
    while(top)pt(st[top--]);pt(32);
}
template<class T>
inline void println(T x){
    static char st[70];short top=0;
    if(x<0)pt('-');
    do{st[++top]=x>=0?(x%10+48):(-(x%10)+48),x/=10;}while(x);
    while(top)pt(st[top--]);pt(10);
}
inline void put_str(string s){
    int siz=s.size();
    for(int i=0;i<siz;++i) pt(s[i]);
    printf("\n");
}
ll a,b,p;
inline ll ksm(ll a,ll b,ll p){
    ll res=1;
    a%=p;
    while(b){
        if(b&1)res=res*a%p;
        a=a*a%p,b>>=1;
    }
    return res;
}
template<class P,class Q>
struct Hash_Table{
    static const int mod=5e5+9;
    struct Node{
        P key;Q val;
        int nxt;
    }node[N];
    int head[mod+5],cnt;
    inline int get_hash(P x){return (x%mod+mod)%mod;}
    inline void clear(){
        for(int i=1;i<=cnt;++i){
            head[get_hash(node[i].key)]=0,node[i].key=node[i].nxt=0;
            Q tmp;
            swap(node[i].val,tmp);
        }
        cnt=0;
    }
    inline void insert(P key){
        int val=get_hash(key);
        node[++cnt].key=key;
        node[cnt].nxt=head[val];
        head[val]=cnt;
    }
    inline Q &operator[](const P key){
        int u=get_hash(key);
        for(int i=head[u];i;i=node[i].nxt){
            P v=node[i].key;
            if(v==key)return node[i].val;
        }
        insert(key);
        return node[cnt].val;
    }
    inline bool count(const P key){
        int u=get_hash(key);
        for(int i=head[u];i;i=node[i].nxt){
            P v=node[i].key;
            if(v==key)return 1;
        }
        return 0;
    } 
};
Hash_Table<ll,int>mp;
inline ll BSGS(ll a,ll b,ll p){
    // 求解 a^x % p = b
    ll s=sqrt(p),val=1;
    for(int i=0;i<s;++i) mp[val]=i,val=val*a%p;
    if(mp.count(b))return mp[b];
    ll inv=ksm(ksm(a,s,p),p-2,p),tmp=b%p;
    for(int i=1;1ll*i*s<=p;++i){
        // 查找第 i 组（从 0 开始，这样就查找 b/a^si 即可） 
        tmp=(tmp*inv)%p;
        if(mp.count(tmp))return mp[tmp]+i*s;
    }
    return -1;
}
signed main(){
    p=read(),a=read(),b=read();
    int ans=BSGS(a,b,p);
    if(ans==-1)printf("no solution\n");
    else println(ans);
    return 0;
}

也可以换成 unorderep_map 或者 gp/cc_hash_table，但是性能会差一点。

组合数学

基础知识

加法原理：两种选择是有关系/矛盾的，不能同时选择（同一阶段）。

乘法原理：两种选择是没有关系的，可以同时选择（不同阶段）。

排列：从 n 个人里选出 m 个人站成一列（考虑顺序, 12 和 21 是两种情况），有几种选法？

对于第一个人，有 n 种选法（都可以选），第二个人可以选 n-1 种（第一个人已经选了）……一直到第 m 个人，有 n-m+1 种选法。

由于每个人的选择是独立的，所以乘法原理：P(n,m)=n\times n-1\times n-2\times ... \times (n-m+1)=\dfrac{n!}{m!}（排列数）

ZHX 语录：小学奥数里用 A 是因为怕我们不认识 P。

组合：从 n 个人里选出 m 个人，不考虑顺序（12 和 21 是一种情况）。要排除重复的方案数，有 \binom{n}{m}=\dfrac{n!}{m!(n-m)!}

组合

定义式：\binom{n}{m}=\dfrac{n!}{m!(n-m)!}

现在来考虑一些性质：

• 递推式：\binom{n}{m}=\binom{n-1}{m}+\binom{n-1}{m}（DP 思想，考虑第 n 个选不选，同一阶段，加法原理）。

• 当 n\ge1 时，\sum_{i=0}^{n}(-1)^i\binom{n}{i}=0。

相当于证明选偶数个物品的方案数等于选奇数个物品的方案数的方案数。

考虑递推式：\binom{n}{m}=\binom{n-1}{m}+\binom{n-1}{m}。那么，选奇数个物品的方案数实际上就是

弱智档，输出 0 即可。

• 部分分二：n,m\le 1000，p 没有限制

递推求组合数即可。

• 部分分三：n,m\le10^6，p 为质数

递推求出阶乘和阶乘逆元，然后套公式即可。

• 部分分四：n\le 10^9,m\le 10^3，p 没有限制

所以，还可以得出

独立事件：如果 A 事件不影响 B 事件发生，那么就称它们为独立事件。实际上，如果 A,B 独立，有：

实际上，对比条件概率的公式，有：

这也可以解释：由于这两个事件独立，所以 B 事件的发生对 A 事件的发生没有影响。

比如，扔两次骰子，这两个事件就是互相独立的。

期望

还是扔骰子，每个数的概率都是 \dfrac{1}{6}，求扔出来的数的期望。

期望，实际上就是每个事件的权值（比如扔骰子的期望） \times 概率再加起来，所以扔骰子的期望值就是：

那如果是扔出来的数的平方的期望呢？同理，为：

接下来，介绍一下期望的性质：

• 期望的和 = 和的期望。比如我们有两颗骰子：第一颗是正常的，每种情况概率都是 \dfrac{1}{6}；第二颗灌了铅，扔出 6 的概率是 \dfrac{95}{100}，其他都是 \dfrac{1}{100}，那么：

实际上，这在任何情况下都适用：

这也是对的。这样，就可以大大降低运算量。

严谨证明？不会

习题

一、课件 P111

某小区有男 51 人女 49 人，男人色盲的概率为 \dfrac{2}{100}，女人色盲的概率为 \dfrac{0.25}{100}，现在抽到一个色盲，问这个人是男人的概率。

很明显，这是求 P(\text{这个人是男人}|\text{这个人是色盲})，根据条件概率的公式，有：

剩下懒得约分了。

二、课件 P117 - 小葱过河

问题在于比较 (\dfrac{99}{100})^{100} 和 (\dfrac{999}{1000})^{1000} 的大小。

发现指数可以约掉，问题变为比较 \dfrac{99}{100} 和 (\dfrac{999}{1000})^{10} 的大小。

考虑放缩法，由于 \dfrac{999}{1000}>\dfrac{998}{999}>...>\dfrac{990}{991}，所以

约分，得：

所以，答案是不过河走那 1000 个石头。

三、课件 Question 11 - 小胡抛硬币

操作二和操作三之间的那个点非常关键，我们考虑枚举它的坐标。

显然，最坏情况下会走无限步，所以 x,y 都有可能达到无限。

• 首先，我们需要走到 (x,0)，这意味着我们需要连续抛 x 次反面，并再抛一次正面停下来。这部分的概率为 (1-p)^xp

• 然后，我们需要走到 (x,y)，这意味着我们需要连续抛 y 次反面，并再抛一次正面停下来。这部分的概率为 (1-p)^yp

• 最后，我们需要原路返回。很明显，我们一共需要抛 x+y 次，且需要抛 x 次正面， y 次反面，方案数为 \binom{x+y}{x}。这部分的概率为 q^x(1-q)^y\binom{x+y}{x}。这实际上就是二项式分布。

把上面三个部分乘起来（不是一个阶段），然后枚举 x,y，有：

我们考虑用求和变形技巧化简式子。

首先，把能合并的合并，有：

然后我们介绍第四个求和变形的技巧：

• 换元法：我们令 t=x+y，然后枚举 t 和 x。有：

分离无关变量：

这是什么？二项式定理！

所以：

剩下的是什么？是等比数列。

我们可以使用等比数列求和公式：

所以：

分子分母取相反数，得：

在极限意义下，(1-p)^{\infty} 可以看做 0，所以：

问题来了：当 p=0 时，\dfrac{1}{p} 无意义，怎么办？

实际上，当 p=0 时，在使用等比数列求和之前就可以变为 0，所以还是对的。

四、BZOJ 2396

可以在 Dark BZOJ 上交。网址。

梯子在这里。

我们可以构造一个 n\times 1 的矩阵（向量） D，很明显，如果 A\times B=C，那么：

首先，C\times D 是 O(n^2) 的，对于左边的式子，由于矩阵有结合律，所以我们可以先 O(n^2) 算出 B\times D，再 O(n^2) 算出 A\times B\times D，可以接受。

但是反过来对吗？

明显不对，比如 D 是零矩阵。

但是我们多随几次，正确率就很高了。

代码（我随了十次）：

#include<bits/stdc++.h>
//#pragma GCC optimize("Ofast")
#define gt getchar
#define pt putchar
#define y1 y233
//typedef unsigned int uint;
typedef long long ll;
typedef unsigned long long ull;
//typedef __int128 lll;
//typedef __uint128_t ulll;
const int N=1005;
using namespace std;
inline bool __(char ch){return ch>=48&&ch<=57;}
inline ll read(){
    ll x=0;bool sgn=0;char ch=gt();
    while(!__(ch)&&ch!=EOF){sgn|=(ch=='-');ch=gt();}
    while(__(ch)){x=(x<<1)+(x<<3)+(ch-48);ch=gt();}
    return sgn?-x:x;
}
template<class T>
inline void print(T x){
    static char st[70];short top=0;
    if(x<0)pt('-');
    do{st[++top]=x>=0?(x%10+48):(-(x%10)+48),x/=10;}while(x);
    while(top)pt(st[top--]);
}
template<class T>
inline void printsp(T x){
    static char st[70];short top=0;
    if(x<0)pt('-');
    do{st[++top]=x>=0?(x%10+48):(-(x%10)+48),x/=10;}while(x);
    while(top)pt(st[top--]);pt(32);
}
template<class T>
inline void println(T x){
    static char st[70];short top=0;
    if(x<0)pt('-');
    do{st[++top]=x>=0?(x%10+48):(-(x%10)+48),x/=10;}while(x);
    while(top)pt(st[top--]);pt(10);
}
inline void put_str(string s){
    int siz=s.size();
    for(int i=0;i<siz;++i) pt(s[i]);
    printf("\n");
}
int n;
struct Matrix{
    int n,m;
    ll val[N][N];
    Matrix(int _n=0,int _m=0){
        n=_n,m=_m;
        memset(val,0,sizeof(val));
    }
    inline void Reader(){
        for(int i=0;i<n;++i) for(int j=0;j<m;++j) val[i][j]=read();
    }
    inline void Output(){
        for(int i=0;i<n;++i){
            for(int j=0;j<m;++j) printsp(val[i][j]);
            printf("\n");
        }
    }
    inline bool operator==(const Matrix &b)const{
        for(int i=0;i<n;++i) for(int j=0;j<m;++j) if(val[i][j]!=b.val[i][j]) return 0;
        return 1;
    }
}A,B,C;
inline Matrix operator*(const Matrix &a,const Matrix &b){
    Matrix c(a.n,b.m);
    for(int i=0;i<c.n;++i){
        for(int k=0;k<a.m;++k){
            for(int j=0;j<c.m;++j){
                c.val[i][j]+=a.val[i][k]*b.val[k][j];
            }
        }
    }
    return c;
}
mt19937 Rand(time(0));
inline bool check(Matrix a,Matrix b,Matrix c){
    // 检查 a * b 是否 = c
    Matrix Vector(c.m,1); 
    for(int i=0;i<Vector.n;++i) Vector.val[i][0]=Rand();
    return a*(b*Vector)==c*Vector;
}
inline void solve(){
    A.n=C.n=A.m=B.n=B.m=C.m=n;
    A.Reader(),B.Reader(),C.Reader();
    for(int i=1;i<=10;++i){
        if(!check(A,B,C)){
            printf("No\n");
            return;
        }
    }
    printf("Yes\n");
}
signed main(){
    while(~scanf("%d",&n)) solve();
    return 0;
}

注意 check 函数的最后，要加括号，不然复杂度就伪了。

这题在 luogu 上有个原：P10102。

不过这题只需要随一次就行了，题解说错误率是 \frac{1}{\bmod}，没看懂，好深刻。

贴个码：

#include<bits/stdc++.h>
#include<ext/pb_ds/assoc_container.hpp>
#include<ext/pb_ds/tree_policy.hpp>
#include<ext/pb_ds/hash_policy.hpp>
#define gt getchar
#define pt putchar
#define fst first
#define scd second
#define SZ(s) ((int)s.size())
#define all(s) s.begin(),s.end()
#define pb push_back
#define eb emplace_back
typedef long long ll;
typedef double db;
typedef long double ld;
typedef unsigned long long ull;
typedef unsigned int uint;
const int N=3005;
const int mod=998244353;
using namespace std;
using namespace __gnu_pbds;
typedef pair<int,int> pii;
template<class T,class I> inline void chkmax(T &a,I b){a=max(a,(T)b);}
template<class T,class I> inline void chkmin(T &a,I b){a=min(a,(T)b);}
inline bool __(char ch){return ch>=48&&ch<=57;}
template<class T> inline void read(T &x){
    x=0;bool sgn=0;static char ch=gt();
    while(!__(ch)&&ch!=EOF) sgn|=(ch=='-'),ch=gt();
    while(__(ch)) x=(x<<1)+(x<<3)+(ch&15),ch=gt();
    if(sgn) x=-x;
}
template<class T,class ...I> inline void read(T &x,I &...x1){
    read(x);
    read(x1...);
}
template<class T> inline void print(T x){
    static char stk[70];short top=0;
    if(x<0) pt('-');
    do{stk[++top]=x>=0?(x%10+48):(-(x%10)+48),x/=10;}while(x);
    while(top) pt(stk[top--]);
}
template<class T> inline void printsp(T x){
    print(x);
    putchar(' ');
}
template<class T> inline void println(T x){
    print(x);
    putchar('\n');
}
int T,n,a[N][N],b[N][N],c[N][N],d[N];
int res1[N],res2[N],res3[N];
inline void add(int &a,int b){
    a+=b;
    if(a>=mod) a-=mod; 
}
mt19937 Rand(time(0));
inline void solve(){
    read(n);
    for(int i=1;i<=n;++i) for(int j=1;j<=n;++j) read(a[i][j]);
    for(int i=1;i<=n;++i) for(int j=1;j<=n;++j) read(b[i][j]);
    for(int i=1;i<=n;++i) for(int j=1;j<=n;++j) read(c[i][j]);
    for(int i=1;i<=n;++i) d[i]=Rand()%mod,res1[i]=res2[i]=res3[i]=0;
    for(int i=1;i<=n;++i){
        for(int j=1;j<=n;++j){
            add(res1[i],1ll*d[j]*a[j][i]%mod);
        }
    }
    for(int i=1;i<=n;++i){
        for(int j=1;j<=n;++j){
            add(res2[i],1ll*res1[j]*b[j][i]%mod);
        }
    }
    for(int i=1;i<=n;++i){
        for(int j=1;j<=n;++j){
            add(res3[i],1ll*d[j]*c[j][i]%mod);
        }
    }
    for(int i=1;i<=n;++i) if(res2[i]!=res3[i]) return printf("No\n"),void();
    printf("Yes\n");
}
signed main(){
    read(T);
    while(T--) solve();
    return 0;
}

显然 O(n^2)。

但是原来的代码过不了，说明不用封装的时候别 xjb 封装净增常数。