第六章 数学(二)

3 欧拉函数

互质：公因数只有$1$的两个非零自然数，叫做互质数。

欧拉函数的定义：$1\sim N$中与$N$互质的数的个数被称为欧拉函数，记为$\varphi (N)$。

由算数基本定理中，$N=P_1^{{\alpha }_1}\cdot P_2^{{\alpha }_2}\cdot \cdot \cdot P_k^{{\alpha }_k}$，则$$\varphi (N)=N\times \frac{P_1-1}{P_1}\times \frac{P_2-1}{P_2}\times ...\times \frac{P_k-1}{P_k}$$

证明 (容斥原理)：$N$的质因子只有$P_1、P_2、...、P_k$

1. 从$1\sim N$中去掉$P_1、P_2、...、P_k$的所有倍数，会多去掉一些；
2. 加上$P_i\times P_j$的倍数；
3. 减去所有$P_i\times P_j\times P_k$的倍数；
4. …

上面步骤形成的公式为：
$$\begin{array}{rl}N& -\frac{N}{P_1}-\frac{N}{P_2}-...-\frac{N}{P_k}\\ & +\frac{N}{P_1{P}_2}+\frac{N}{P_1{P}_3}+...\\ & -\frac{N}{P_1{P}_2{P}_3}-\frac{N}{P_1{P}_2{P}_4}-...\\ & +\frac{N}{P_1{P}_2{P}_3{P}_4}+...\\ & ...\end{array}$$将$\varphi (N)=N\times \frac{P_1-1}{P_1}\times \frac{P_2-1}{P_2}\times ...\times \frac{P_k-1}{P_k}$展开即得上式。

欧拉函数的应用
欧拉定理：若$a$与$n$互质，那么$$a^{\varphi (n)}\equiv 1(modn)$$
例如$a=5、n=6$，那么有$\varphi (6)=2$，则$5^{2}mod6=1$ 。

欧拉定理推论(费马小定理)：若$a$与$p$互质且$p$为质数的时候，
$$(a^{\varphi (p)}=a^{p-1})\equiv 1(modp)$$

欧拉定理证明：设$1\sim n$中和$n$互质的数为$a_1,a_2,...,a_{\varphi (n)}$，因为$a$与$n$互质，那么有$a{a}_1,a{a}_2,...,a{a}_{\varphi (n)}$同样与$n$互质，且这些数两两互不相同（证明见下面）。因此有$$\begin{array}{rl}& a^{\varphi (n)}\cdot (a_1\cdot a_2\cdot ...\cdot a_{\varphi (n)})\equiv (a_1\cdot a_2\cdot ...\cdot a_{\varphi (n)})(modn)\\ & a^{\varphi (n)}\equiv 1(modn)\end{array}$$

证明$a{a}_1,a{a}_2,...,a{a}_{\varphi (n)}$这些数两两互不相同：若有$a{a}_i\equiv a{a}_j(modn)$，那么就有$a(a_i-a_j)\equiv 0(modn)$，因为$a$与$n$互质，也就有$a_i\equiv a_j$，矛盾。

3.1 朴素法求欧拉函数

ACwing 873

分解质因数时间复杂度$O(\sqrt{N})$，这个题有$n$个数，所以这个题的时间复杂度为$O(n\sqrt{N})$

#include <iostream>
using namespace std;

int phi(int x) {
    int res = x;
    for (int i = 2; i <= x / i; i++)
        if (x % i == 0) {
            res = res / i * (i - 1);
            while (x % i == 0) x /= i;
        }
    if (x > 1) res = res / x * (x - 1);
    return res;
}

int main() {
    int n;
    cin >> n;
    while (n--) {
        int x; cin >> x;
        cout << phi(x) << endl;
    }
    return 0;
}

3.2 筛法求欧拉函数

ACwing 874

根据线性筛法求质数的方法，来求欧拉函数。

代码中几个地方解释：

1. 第18行：当$i$为一个质数时，它的欧拉函数为$i-1$；
2. 第23行：$imod{P}_j=0$，这个时候$primes[j]$为$i$的最小质因子，所以$P_j\cdot i$的质因数和$i$的质因数相等，那么有欧拉函数$$\begin{array}{rl}& \varphi (i)=i\times \frac{P_1-1}{P_1}\times \frac{P_2-1}{P_2}\times ...\times \frac{P_k-1}{P_k}\\ & \varphi (P_j\cdot i)=P_j\times i\times \frac{P_1-1}{P_1}\times \frac{P_2-1}{P_2}\times ...\times \frac{P_k-1}{P_k}\end{array}$$因此有$$\varphi (P_j\cdot i)=P_j\times \varphi (i)$$
3. 第27行：$imod{P}_j\ne 0$，这时$P_j$是$i\times P_j$的最小质因子，但$P_j$不是$i$的最小质因子。那么有欧拉函数$$\begin{array}{rl}& \varphi (i)=i\times \frac{P_1-1}{P_1}\times \frac{P_2-1}{P_2}\times ...\times \frac{P_k-1}{P_k}\\ & \varphi (P_j\cdot i)=P_j\times i\times \frac{P_1-1}{P_1}\times \frac{P_2-1}{P_2}\times ...\times \frac{P_k-1}{P_k}\times \frac{P_j-1}{P_j}\end{array}$$因此有$$\begin{array}{rl}\varphi (P_j\cdot i)& =P_j\times \varphi (i)\times \frac{P_j-1}{P_j}\\ & =(P_j-1)\times \varphi (i)\end{array}$$

时间复杂度：$O(n)$

#include <iostream>

using namespace std;

typedef long long LL;

const int N = 1000010;

int primes[N], cnt; // 存储质数
int euler[N]; // 存储每个数的欧拉函数
bool st[N];

void get_eulers(int n) {
    euler[1] = 1;
    for (int i = 2; i <= n; i++) {
        if (!st[i]) {
            primes[cnt++] = i;
            euler[i] = i - 1; 
        }
        for (int j = 0; primes[j] <= n / i; j++) {
            int t = primes[j] * i;
            st[t] = true;
            if (i % primes[j] == 0) { 
                euler[t] = euler[i] * primes[j];
                break;
            }
            euler[t] = euler[i] * (primes[j] - 1);
        }
    }
}

int main() {
    int n; cin >> n;
    get_eulers(n);
    LL res = 0;
    for (int i = 1; i <= n; i++) res += euler[i];
    cout << res << endl;
    return 0;
}

3.3 可见的点

ACwing 201

假设存在一个点 $(x_0,y_0)$，该点与坐标原点 $(0,0$ 之间的直线斜率为 $\frac{y_0}{x_0}$，且直线方程为 $y=\frac{y_0}{x_0}x$。如果该点为可见点，即该点为从坐标原点沿直线出发所能到达的第一个点。

算法思路：

$(x_0,y_0)$ 为直线上的第一个点 $⟺$ $x_0、y_0$ 互质

证明：

• 假设 $x_0、y_0$ 不互质，那么存在 $(x_0,y_0)=d>0$，即存在点 $(\frac{x_0}{d},\frac{y_0}{d})$ 存在直线上且位于 $(x_0,y_0)$ 左下方，那么 $(x_0,y_0)$ 便不是第一个可见点。
• 假设 $x_0、y_0$ 互质，且存在点 $({x_0}^{\prime },{y_0}^{\prime })$ 在直线上，${x_0}^{\prime }<x_0、{y_0}^{\prime }<y_0$，那么有 $\frac{{y_0}^{\prime }}{{x_0}^{\prime }}=\frac{y_0}{x_0}$ 。因为 $x_0、y_0$ 互质，也就是说分数 $\frac{y_0}{x_0}$ 已经最简，因此不能能存在一个更小的 $({x_0}^{\prime },{y_0}^{\prime })$ 在点 $x_0、y_0$ 左下方，即如果直线上 $x_0、y_0$ 互质，它就一定是直线上的第一个可见点。

因此该问题转换为了：求解$0\le x、y\le N$中有多少对 $(x,y)$ 满足 $x、y$ 互质？

由图关于直线 $y=x$ 对称性，可以只考虑直线 $y=x$ 下面(不包括直线上的唯一一个点 $(1,1)$ )的所有情况。枚举每一个横坐标 $x=2,2,3,...,N$ ，可以得到与每一个横坐标互质的纵坐标个数为 $\varphi (x)$，所以在直线下方与所有横纵标互质的纵坐标个数为 ${\sum }_{x=2}^N\varphi (x)$，因此整个图的可见点个数为$$2\sum _{x=2}^N\varphi (x)+1$$

#include <cstdio>
using namespace std;

const int N = 1010;

int primes[N], cnt;
bool st[N];
int phi[N];

void init(int n) {
    phi[1] = 1;
    for (int i = 2; i <= n; i++) {
        if (!st[i]) {
            primes[cnt++] = i;
            phi[i] = i - 1;
        }
        for (int j = 0; primes[j] * i <= n; j++) {
            st[i * primes[j]] = true;
            if (i % primes[j] == 0) {
                phi[i * primes[j]] = phi[i] * primes[j];
                break;
            }
            phi[i * primes[j]] = phi[i] * (primes[j] - 1);
        }
    }
}

int main() {
    init(N - 1);
    int n, m; scanf("%d", &m);
    for (int T = 1; T <= m; T++) {
        scanf("%d", &n);
        int res = 1;
        for (int i = 1; i <= n; i++) res += phi[i] * 2;
        printf("%d %d %d\n", T, n, res);
    }
    return 0;
}

3.4 最大公约数

ACwing 220

算法思路：

因为$(x,y)=p$，其中$x、y\in [1,N]$且$p$为质数，所以有$(\frac{x}{p},\frac{y}{p})=1$，记$x^{\prime }=\frac{x}{p}、y^{\prime }=\frac{y}{p}$，那么原问题就转换为：有多少个$(x^{\prime },y^{\prime })$，其中$1\le x^{\prime }、y^{\prime }\le \lfloor \frac{N}{p}\rfloor$，满足$(x^{\prime },y^{\prime })=1$，即$x^{\prime }、y^{\prime }$互质。

对于$p$可以枚举$[1,N]$的所有质数$p$，到这里这个问题就和上一题的处理方法基本一样。只是这里$\lfloor \frac{N}{p}\rfloor$是变化的，为了快速求解$\varphi (1)+\varphi (2)+...+\varphi (\lfloor \frac{N}{p}\rfloor )$，可以在筛法求解欧拉函数之后求一个$\varphi$的前缀和数组$s[]$，其中定义$\varphi (1)=0$。

#include <cstdio>

typedef long long LL;
const int N = 1e7 + 10;

int primes[N], cnt;
bool st[N];
int phi[N];
LL s[N]; // phi前缀和数组 假定phi(1) = 0

void init(int n) {
    for (int i = 2; i <= n; i++) {
        if (!st[i]) {
            primes[cnt++] = i;
            phi[i] = i - 1;
        }
        for (int j = 0; primes[j] * i <= n; j++) {
            st[primes[j] * i] = true;
            if (i % primes[j] == 0) {
                phi[i * primes[j]] = phi[i] * primes[j];
                break;
            }
            phi[i * primes[j]] = phi[i] * (primes[j] - 1);
        }
    }

    for (int i = 1; i <= n; i++) s[i] = s[i - 1] + phi[i];
}

int main() {
    int n; scanf("%d", &n);
    init(n);
    LL res = 0;
    for (int i = 0; i < cnt; i++) {
        int p = primes[i];
        res += s[n / p] * 2 + 1;
    }
    printf("%lld\n", res);
    return 0;
}

4 快速幂 (欧拉降幂)

下面所有的 $log$ 均以2为底

用于在时间复杂度为$O(lo{g}^k)$以内，快速求出 $a^{k}modp$ 的结果，其中$1\le a,p,k\le 10^9$。

原理：
先预处理出$lo{g}^k$个值：$a^{2^0}modp、a^{2^1}modp、...、a^{2^{lo{g}^k}}modp$，其中每个数的求解方法为
$$\begin{array}{rl}& a^{2^0}=a^1\\ & a^{2^1}=a^{2^0\cdot 2}=(a^{2^0}{)}^2\\ & a^{2^2}=(a^{2^1}{)}^2\\ & ...\\ & a^{2^{lo{g}^k}}=(a^{2^{lo{g}^{k-1}}}{)}^2\end{array}$$所以从前往后平方$k$次就可以求出这$lo{g}^k$个数。

然后将$a^k$拆开成$$\begin{array}{rl}{a}^k& =a^{2^{x_1}}\cdot a^{2^{x_2}}\cdot ...\cdot a^{2^{x_t}}\\ & =a^{2^{x_1}+2^{x_2}+...+2^{x_t}}\end{array}$$即将$k$表示成二进制形式$$k=2^{x_1}+2^{x_2}+...+2^{x_t}$$对于其中的每个数，可以从预处理出来的数中找到。

时间复杂度$O(lo{g}^k)$

4.1 快速幂

ACwing 875

#include <iostream>
#include <algorithm>

using namespace std;

typedef long long LL;

LL qmi(int a, int b, int p) {
    LL res = 1 % p;
    while (b) {
        if (b & 1) res = res * a % p;
        a = (LL) a * a % p;
        b >>= 1;
    }
    return res;
}

int main() {
    int n; scanf("%d", &n);
    while (n--) {
        int a, b, p;
        scanf("%d%d%d", &a, &b, &p);
        printf("%lld\n", qmi(a, b, p));
    }
    return 0;
}

4.2 快速幂求逆元 (将除法变成乘法)

ACwing 876

乘法逆元： 若整数$b，m$互质，并且对于任意的整数$a$，如果满足 $b|a$ ($\frac{a}{b}$为一个整数)，则存在一个整数$x$，使得 $$a/b\equiv a\times x(modm)$$则称$x$为$b$模$m$的乘法逆元，记为$b^{-1}(modm)$。对上式变形有$$b\cdot b^{-1}\equiv 1(modm)$$

对于题目就转换为，找一个数$x$使 $b\cdot x\equiv 1(modp)$，且由题可知$p$为质数，即$p\ge 2$。如果$b$和$p$互为质数，由费马小定理可知$$b^{p-1}=b\cdot b^{p-2}\equiv 1(modp)$$所以对于一个质数$b$而言，它的乘法逆元即为$b^{p-2}$。

$b$存在乘法逆元的充要条件是$b$与模数$m$互质。当模数$m$为质数时，$b^{m-2}$即为$b$的乘法逆元。

对于本题就是求一个快速幂$a^{p-2}(modp)$。

#include <iostream>
#include <algorithm>

using namespace std;

typedef long long LL;

LL qmi(int a, int b, int p) {
    LL res = 1;
    while (b) {
        if (b & 1) res = res * a % p;
        a = a * (LL) a % p;
        b >>= 1;
    }
    return res;
}

int main() {
    int n; scanf("%d", &n);
    while (n--) {
        int a, p; scanf("%d%d", &a, &p);
        if (a % p == 0) puts("impossible");
        else printf("%lld\n", qmi(a, p - 2, p));
    }
    return 0;
}

4.3 序列的第k个数

ACwing 1289

算法思路：

1. 除了全等数列(数列中的每个数值都相等)，一个数列的前三项会不会既是等差数列，又是等比数列？
   假设数列前三项为$a、b、c$：
   如果该数列为等差数列，则应该满足$a+c=2b$；
   如果该数列为等比数列，则应该满足$ac=b^2$；
   将$a=2b-c$带入$ac=b^2$可得$$\begin{array}{rl}& c(2b-c)=b^2\\ \Rightarrow & b^2-2bc+c^2=0\\ \Rightarrow & (b-c{)}^2=0\\ \Rightarrow & b=c\end{array}$$即有$a=c=b$。因此，如果一个数列既为等差数列又为等比数列，那么它一定是一个全等数列。也就是说，除了全等数列之外，不可能存在一个数列既是等差数列，又是等比数列。
2. 情况1：如果该数列为一个等差数列，则公差$d=b-a$，第$k$项为$a+d(k-1)$；
3. 情况2：如果该数列为一个等比数列，则公比为$p=\frac{b}{a}$，第$k$项为$a{p}^{k-1}=a(\frac{b}{a}{)}^{k-1}$，其中因为序列为整数数列，所以可以知道 $\frac{b}{a}$ 一定是一个整数。

#include <iostream>
#include <algorithm>

using namespace std;

typedef long long LL;

const int mod = 200907;

int qmi(int a, int k) {
    int res = 1;
    while (k) {
        if (k & 1) res = (LL) res * a % mod;
        a = (LL) a * a % mod;
        k >>= 1;
    }
    return res;
}

int main() {
    int n; cin >> n;
    while (n--) {
        int a, b, c, k; cin >> a >> b >> c >> k;
        if (a + c == b * 2) cout << (a + (b - a) * (LL) (k - 1)) % mod << endl;
        else cout << (LL) a * qmi(b / a, k - 1) % mod << endl;
    }
    return 0;
}

4.4 越狱

ACwing 1290

算法思路：

发生越狱的情况 = 每个犯人的所有信仰(全集) - 不发生越狱的情况

发生越狱的情况 = $m^n-m(m-1{)}^{n-1}$，使用快速幂即可求解。

注意最终结果取余后应该为正数！

#include <cstring>
#include <iostream>
#include <algorithm>

using namespace std;

typedef long long LL;

const int mod = 100003;

int qmi(int a, LL k) {
    int res = 1;
    while (k) {
        if (k & 1) res = (LL) res * a % mod;
        a = (LL) a * a % mod;
        k >>= 1;
    }
    return res;
}

int main() {
    int m; LL n; cin >> m >> n;
    cout << (qmi(m, n) - (LL) m * qmi(m - 1, n - 1) % mod + mod) % mod << endl;
    return 0;
}

4.5 矩阵乘法

矩阵乘法通常与快速幂结合考察

矩阵乘法基本模板：假设有矩阵$P_{A\times B}、Q_{B\times C}$，两个矩阵相乘$$P_{A\times B}\times Q_{B\times C}=R_{A\times C}$$代码如下：

for (int i = 1; i <= A; i ++ )
    for (int j = 1; j <= C; j ++ )
        for (int k = 1; k <= B; k ++ )
            R[i][j] += P[i][k] * Q[k][j];

矩阵乘法重要性质：结合律，即$A\cdot B\cdot C=A\cdot (B\cdot C)$

注意：矩阵乘法没有交换律

该性质的应用：假设需要求解$A\cdot \begin{array}{c}n\\ \overbrace{x\cdot x\cdot x\cdots x}\end{array}=A\cdot (\begin{array}{c}n\\ \overbrace{x\cdot x\cdot x\cdots x}\end{array})=A\cdot x^n$，这里可以考虑快速幂。因此可以使用矩阵结合律和快速幂在$O(lo{g}^n)$的时间复杂度内解决类似于$A\cdot \begin{array}{c}n\\ \overbrace{x\cdot x\cdot x\cdots x}\end{array}$的问题。

4.5.1 斐波那契前 n 项和

ACwing 1303

首先构造行向量(列向量也可以)$F_n=[f_n,f_{n+1}],F_{n+1}=[f_{n+1},f_{n+2}]$，观察存在矩阵$A$使得$F_{n+1}=F_n\cdot A$。显然有$$[f_n,f_{n+1}]\left[\begin{array}{cc}0& 1\\ 1& 1\end{array}\right]=[f_{n+1},f_{n+2}]$$因此可得矩阵$$A=\left[\begin{array}{cc}0& 1\\ 1& 1\end{array}\right]$$那么就有$$F_n=F_0\cdot A^n$$矩阵乘法没有交换律，但是有结合律，所以可以将$A^n$拆开，如$$A^n=A^1\times A^2\times A^4\times A^8\cdots$$

而这个题目是求前$f_n$项之和，那么构造向量$F_n=[f_n,f_{n+1},S_n]$，$F_{n+1}=[f_{n+1},f_{n+2},S_{n+1}]$，其中$S_i$表示前$i$项之和。假设存在矩阵$A$使得$F_n\cdot A=F_{n+1}$，那么可以有$$[f_n,f_{n+1},S_n]\left[\begin{array}{ccc}0& 1& 0\\ 1& 1& 1\\ 0& 0& 1\end{array}\right]=[f_{n+1},f_{n+2},S_{n+1}]$$因此可以得到矩阵$$A=\left[\begin{array}{ccc}0& 1& 0\\ 1& 1& 1\\ 0& 0& 1\end{array}\right]$$于是就有$F_n=F_0\cdot A^n=F_1\cdot A^{n-1}$。对于$A^{n-1}$使用快速幂计算。

#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;

typedef long long LL;
const int N = 3;

int n, m;

void mul(int c[], int a[], int b[][N]) {
    int temp[N] = {0};
    for (int i = 0; i < N; i++)
        for (int j = 0; j < N; j++)
            temp[i] = (temp[i] + (LL) a[j] * b[j][i]) % m;

    memcpy(c, temp, sizeof temp);
}

void mul(int c[][N], int a[][N], int b[][N]) {
    int temp[N][N] = {0};
    for (int i = 0; i < N; i++)
        for (int j = 0; j < N; j++)
            for (int k = 0; k < N; k++)
                temp[i][j] = (temp[i][j] + (LL) a[i][k] * b[k][j]) % m;

    memcpy(c, temp, sizeof temp);
}

int main() {
    cin >> n >> m;

    int f1[N] = {1, 1, 1};
    int a[N][N] = {
            {0, 1, 0},
            {1, 1, 1},
            {0, 0, 1}
    };

    // 求A^n
    n--;
    while (n) {
        if (n & 1) mul(f1, f1, a);  // res = res * a
        mul(a, a, a);  // a = a * a
        n >>= 1;
    }

    cout << f1[2] << endl;

    return 0;
}

4.5.2 佳佳的斐波那契

ACwing 1304

这个题是上一题的扩展。

算法思路：

由上一题可知 $F_n=[f_n,f_{n+1},S_n]$，这道题给出 $T_n=1\cdot F_1+2\cdot F_2+\cdots n\cdot F_n$

注意：要使用$A\cdot \begin{array}{c}n\\ \overbrace{x\cdot x\cdot x\cdots x}\end{array}$解决问题，就必须要保证$A$中是没有变量的。而 $T_n={\sum }_{i=1}^{n}i\cdot F_i$ 中$F_i$前面的$i$是变量，因此要想办法消掉。

做如下变换：$$\begin{array}{rlr}n\cdot S_n-T_n& =(n-1)F_1+(n-2)F_2+\cdots F_{n-1}& \text{①}\\ (n+1)\cdot S_{n+1}-T_{n+1}& =n{F}_1+(n-1)F_2+\cdots F_n& \text{ ②}\\ ②-①& =S_n& \end{array}$$因此，令$P_n=n\cdot S_n-T_n$，则$$\begin{array}{rl}& T_n=n\cdot S_n-P_n\\ & P_n=P_{n-1}+S_{n-1}\\ & S_n=S_{n-1}+f_n\\ & f_n=F_{n-1}+f_{n-2}\end{array}$$到此，重新考虑本题，定义$$\{\begin{array}{ll}& F_n=[f_n,f_{n+1},S_n,P_n]\\ & F_{n+1}=[f_{n+1},f_{n+2},S_{n+1},P_{n+1}]\end{array}$$其中$P_n=n\cdot S_n-T_n$，则有$$[f_n,f_{n+1},S_n,P_n]\left[\begin{array}{cccc}0& 1& 0& 0\\ 1& 1& 1& 0\\ 0& 0& 1& 1\\ 0& 0& 0& 1\end{array}\right][f_{n+1},f_{n+2},S_{n+1},P_{n+1}]$$因此有$$A=\left[\begin{array}{cccc}0& 1& 0& 0\\ 1& 1& 1& 0\\ 0& 0& 1& 1\\ 0& 0& 0& 1\end{array}\right]$$则$F_n=F_1\cdot A^{n-1}$，$F_1=[1,1,1,0]$。

时间复杂度：快速幂$O(lo{g}^n)$，矩阵乘法$O(n^3)$，总共时间复杂度为$O(4^{3}lo{g}^n)$。

Tips：在写代码会遇到$A_{1\times 4}\times B_{4\times 4}$和$C_{4\times 4}\times B_{4\times 4}$两种情况，为了节省代码，可以将$A_{1\times 4}$扩充成$A_{4\times 4}$，其中第二行到第四行的所有值为$0$。比如对于$F_1=[1,1,1,0]$有
$$F_1=\left[\begin{array}{cccc}1& 1& 1& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right]$$

#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;

typedef long long LL;
const int N = 4;

int n, m;

void mul(int c[][N], int a[][N], int b[][N]) { // c = a * b
	// c 和 a 有可能相同，因此定义一个备份数组
	// static表示变量只会被定义一次，以后每次访问函数mul它就不需要重新定义
    static int t[N][N]; 
    // 这里不能写 sizeof c，传入的c为一个指针，不能通过指针获得数组长度
    memset(t, 0, sizeof t); 

    for (int i = 0; i < N; i++)
        for (int j = 0; j < N; j++)
            for (int k = 0; k < N; k++)
                t[i][j] = (t[i][j] + (LL) a[i][k] * b[k][j]) % m;

    memcpy(c, t, sizeof t);
}

int main() {
    cin >> n >> m;

    // {fn, fn+1, sn, pn}
    // pn = n * sn - tn
    int f1[N][N] = {1, 1, 1, 0};
    int a[N][N] = {
            {0, 1, 0, 0},
            {1, 1, 1, 0},
            {0, 0, 1, 1},
            {0, 0, 0, 1},
    };

	// 快速幂
    int k = n - 1;
    while (k) {
        if (k & 1) mul(f1, f1, a);  // f1 = f1 * a
        mul(a, a, a);  // a = a * a
        k >>= 1;
    }

    cout << (((LL) n * f1[0][2] - f1[0][3]) % m + m) % m << endl; // Tn = n x Sn - Pn

    return 0;
}

4.5.3 GT考试

ACwing 1305

类似于ACwing 1052 设计密码，那一个题$n\in [1,50]$，这个题$n\in [1,10^9]$。本题和设计密码不包含的子串只有一个，如果不包含的子串有多个，则可以参考ACwing 1053 修复DNA。

算法思路：DP状态机

记不吉利数字为$S$。

集合：$f(i,j)$表示所有长度为$i$、不包含$S$，且末尾部分与$S$前缀匹配的最大长度为$j$的所有字符串的集合。

属性：数量

状态计算：

$$f[i+1][k]=a[0][k]\ast f[i][0]+a[1][k]\ast f[i][1]+\cdots +a[m-1][k]\ast f[i][m-1]$$

发现$f[i]$与$f[i+1]$的关系是线性的，所以可以用向量来计算，$F[i+1]=F[i]\ast A$。其中$A$向量中$a[i][j]$表示在总长度相邻的两个串中，将末尾部分与不吉利数字重复部分的长度从$i$变成$j$的方案数，这是一个定值，可以预处理。

比如$$\begin{array}{rl}& f(i+1,0)=a_{00}f(i,0)+a_{10}f(i,1)+\cdots \\ & f(i+1,1)=a_{01}f(i,1)+a_{11}f(i,1)+\cdots \\ & \cdots \end{array}$$那么可以记$F_{i+1}=[f(i+1,0),f(i+1,1),\cdots ]$，$F_i=[f(i,0),f(i,1),\cdots ]$，其所有系数记为矩阵$A$，则有$F_{i+1}=F_i\cdot A$，且$A$的取值与串的位数$i$无关，所以可以推出$F_n=F_{n-1}\cdot A=F_1\cdot A^{n-1}$。

时间复杂度：

快速幂$O(lo{g}^n)$，其中$n\in [1,10^9]$，矩阵乘法$20^3$，因此总共时间复杂度为$O(20^{3}lo{g}^n)$。

#include <cstring>
#include <iostream>
#include <algorithm>

using namespace std;

const int N = 25;

int n, m, mod;
char str[N];
int ne[N];
int a[N][N];

void mul(int c[][N], int a[][N], int b[][N])  // c = a * b
{
    static int t[N][N];
    memset(t, 0, sizeof t);

    for (int i = 0; i < m; i++)
        for (int j = 0; j < m; j++)
            for (int k = 0; k < m; k++)
                t[i][j] = (t[i][j] + a[i][k] * b[k][j]) % mod;

    memcpy(c, t, sizeof t);
}

int qmi(int k) {
    int f0[N][N] = {1};
    while (k) {
        if (k & 1) mul(f0, f0, a);  // f0 = f0 * a
        mul(a, a, a);  // a = a * a
        k >>= 1;
    }

    int res = 0;
    for (int i = 0; i < m; i++) res = (res + f0[0][i]) % mod;
    return res;
}

int main() {
    cin >> n >> m >> mod;
    cin >> (str + 1);

    // kmp
    for (int i = 2, j = 0; i <= m; i++) {
        while (j && str[j + 1] != str[i]) j = ne[j];
        if (str[j + 1] == str[i]) j++;
        ne[i] = j;
    }

    // 初始化A[i][j]
    for (int j = 0; j < m; j++)
        for (int c = '0'; c <= '9'; c++) {
            int k = j;
            while (k && str[k + 1] != c) k = ne[k];
            if (str[k + 1] == c) k++;
            if (k < m) a[j][k]++;
        }
    
    // F[n] = F[0] * A^n
    cout << qmi(n) << endl;

    return 0;
}

5 扩展欧几里得算法

裴蜀定理：对于任意一对正整数$a、b$，那么一定存在非零整数$x、y$，使得$ax+by=(a,b)$(最大公约数)，且$(a,b)$是$a、b$能凑出来的最小正整数。

证明：
设$ax+by=d$，因为$a、b$是$(a,b)$的倍数，所以$d$也一定是$(a,b)$的倍数，那么最小的正整数就应该是$1$倍的情况下，也就是$(a,b)$。

对于求解$x、y$就需要使用扩展欧几里得算法。

5.1 扩展欧几里得算法

ACwing 877

对于代码中的几点解释：

1. 当$b=0$的时候，$(a,0)=a$，要求一个$x、y$使得$ax+by=(a,0)=a$，显然有$x=1,y=0$；
2. 当递归函数进行到第12行的时候，就应该存在$$by+(amodb)x=(b,a)=d$$因为$amodb=a-\lfloor \frac{a}{b}\rfloor \cdot b$，带入上式后整理得$$ax+b(y-\lfloor \frac{a}{b}\rfloor \cdot x)=d$$所以就有第12行的代码，$x$不变，$y=y-\lfloor \frac{a}{b}\rfloor \cdot x$。

#include <iostream>
#include <algorithm>

using namespace std;

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

int main() {
    int n;
    scanf("%d", &n);
    while (n--) {
        int a, b;
        scanf("%d%d", &a, &b);
        int x, y;
        exgcd(a, b, x, y);
        printf("%d %d\n", x, y);
    }
    return 0;
}

5.2 线性同余方程

ACwing 878

对等式做变形有：存在一个整数$y$，使得$$\begin{array}{rl}ax& \equiv b(modm)\\ & =my+b\\ ax& -my=b\\ ax& +m{y}^{\prime }=b\end{array}$$其中，$y^{\prime }=-y$。

原问题就转换为了求一个$x、y^{\prime }$使$ax+m{y}^{\prime }=b$有解，且有解的条件为$(a,m)|b$。使用扩展欧几里得算法可以求出$ax+m{y}^{\prime }=(a,m)=d$，将系数扩大$b/d$倍，就可以求出$ax+m{y}^{\prime }=b$的解。

#include <iostream>
#include <algorithm>

using namespace std;

typedef long long LL;

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

int main() {
    int n; scanf("%d", &n);
    while (n--) {
        int a, b, m; scanf("%d%d%d", &a, &b, &m);
        int x, y;
        int d = exgcd(a, m, x, y);
        if (b % d) puts("impossible");
        else printf("%d\n", (LL) b / d * x % m);
    }
    return 0;
}

5.3 同余方程

ACwing 203

算法思路：

$$\begin{array}{rl}& ax\equiv 1(modb)\\ \Rightarrow & \exists y\in Z，ax+by=1\end{array}$$通过扩展欧几里得算法可以求出一组解$(x_0,y_0)$，因为$(a,b)=d=1$，故其所有解为$$\{\begin{array}{l}x=x_0+k\frac{b}{d}=x_0+kb\\ y=y_0-k\frac{a}{d}=y_0-ka\end{array}$$因为 $ax\equiv 1(modb)$，因此$x=x_0+kb\ne 0$，$x$的最小正整数解就为 $x_{0}modb$。

#include <cstdio>
using namespace std;

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

int main() {
    int a, b; scanf("%d%d", &a, &b);
    int x, y; exgcd(a, b, x, y);
    printf("%d\n", (x % b + b) % b);
    return 0;
}

5.4 青蛙约会

ACwing 222

算法思路：

假定两只青蛙 $A、B$ 初始坐标分别为 $a、b$，不妨设 $a<b$，若它们跳 $x$ 次后相遇。那么$A$要追$B$，每跳一次 $A$ 追赶 $B$的长度为 $(m-n)$ 米，那么跳了 $y$ 圈之后(每一圈长度为$L$米)，相遇时存在如下关系：$$\begin{array}{rl}& (m-n)\cdot x=b-a+y\cdot L\\ \Rightarrow & (m-n)\cdot x-L\cdot y=b-a\\ \Rightarrow & (m-n)\cdot x+L\cdot y^{\prime }=b-a\end{array}$$其中 $y^{\prime }=-y$。因为式子中 $(m-n)、L、b-a$ 均为常数，仅存在变量 $x、y$，所以可以使用扩展欧几里得算法先求解方程$$(m-n)\cdot x+L\cdot y^{\prime }=d$$其中$d=(m-n,L)$。得到该方程一个解 $(x_0,y_0)$，那么也就可以得到该方程的所有解$$\{\begin{array}{l}x=x_0-k\frac{L}{d}\\ y=y_0+k\frac{m-n}{d}\end{array}$$因为最终求解最小正整数解，所以 $x=x_0-k\frac{L}{d}\ne 0$，这里和上面一题最后一样，$x=x_{0}mod\frac{L}{d}$ ，最后将$(x,y)$扩大$\frac{b-a}{d}$倍即可。

#include <cstdio>
#include <algorithm>
using namespace std;

typedef long long LL;

LL exgcd(LL a, LL b, LL &x, LL &y) {
    if (!b) {
        x = 1, y = 0;
        return a;
    }
    LL d = exgcd(b, a % b, y, x);
    y -= a / b * x;
    return d;
}

int main() {
    LL a, b, m, n, L;
    scanf("%lld%lld%lld%lld%lld", &a, &b, &m, &n, &L);
    LL x, y;
    LL d = exgcd(m - n, L, x, y);
    if ((b - a) % d) puts("Impossible");
    else {
        x *= (b - a) / d;
        LL t = abs(L / d);
        printf("%lld\n", (x % t + t) % t);
    }
    return 0;
}

5.5 最幸运的数字

ACwing 202

算法思路：

假设有$x$个$8$，则有$$\begin{array}{rl}& \begin{array}{c}x\\ \overbrace{8\cdots 8}\end{array}\Rightarrow 8\times \begin{array}{c}x\\ \overbrace{1\cdots 1}\end{array}\Rightarrow 8\times \frac{\begin{array}{c}x\\ \overbrace{9\cdots 9}\end{array}}{9}\Rightarrow 8\times \frac{10^x-1}{9}\end{array}$$题目中问是否存在 $x$ 个 $8$ 能整除 $L$，那么就转换为$$L|8\frac{10^x-1}{9}\Rightarrow 9L|8(10^x-1)$$令 $d=(L,8)$ 后有$$\frac{9L}{d}|10^x-1$$

解释一下这一步：记 $A=10^x-1$，则原式变成 $$9L|8A$$因为 $gcd(8,9)=1$，令 $gcd(8,9L)=gcd(8,L)=d$，于是有$$\frac{9L}{d}|\frac{8}{d}\cdot A$$又因为$$gcd(\frac{9L}{d},\frac{8}{d})=gcd(\frac{L}{d},\frac{8}{d})=1$$故$$\frac{9L}{d}|A$$

因为 $\frac{9L}{d}$ 为一个常数，令 $C=\frac{9L}{d}$，则有$$\begin{gathered} C|10^x-1 \\ \Rightarrow 10^x\equiv 1(modC) \end{gathered}$$

欧拉定理推论(费马小定理)：若$a$与$p$互质且$p$为质数的时候
$$(a^{\varphi (p)}=a^{p-1})\equiv 1(modp)$$

这个题要求解最小的 $x$，先给出结论：满足方程 $10^x\equiv 1(modC)$ 最小的正整数 $x$ 一定满足 $x|\varphi (C)$。

结论证明：
假设 $x$ 是满足方程 $10^x\equiv 1(modC)$ 最小的正整数，但是 $x$ 并不能整除 $\varphi (C)$，那么必然满足 $\varphi (C)=q\cdot x+r$，其中 $0<r<x$。因为由 $10^x\equiv 1⟺10^{qx}\equiv 1$，$10^{\varphi (C)}\equiv 1⟺10^{q\cdot x+r}\equiv 1$，那么可以得到 $10^r\equiv 1$ 且 $0<r<x$，这里的找到了一个更小的正整数满足方程 $10^x\equiv 1(modC)$，与假设矛盾。

时间复杂度：$O(\sqrt{n})$，主要是求约数。

#include <cstdio>
#include <iostream>
#include <algorithm>
using namespace std;

typedef long long LL;

LL qmul(LL a, LL k, LL b) {
    LL res = 0;
    while (k) {
        if (k & 1) res = (res + a) % b;
        a = (a + a) % b;
        k >>= 1;
    }
    return res;
}

// 因为 max(L) = 2e9，C = 9 * L / d，C 会达到 10e10 级别
// 在qmi中，如果再平方，那么达到 10e20 会爆 long long
// 因此要使用 龟速乘qmul 来代替乘法
LL qmi(LL a, LL k, LL b) {
    LL res = 1;
    while (k) {
        if (k & 1) res = qmul(res, a, b);
        a = qmul(a, a, b);
        k >>= 1;
    }
    return res;
}

LL get_euler(LL C) {
    LL res = C;
    for (LL i = 2; i <= C / i; i++)
        if (C % i == 0) {
            while (C % i == 0) C /= i;
            res = res / i * (i - 1); // 先除后乘，防止溢出
        }
    if (C > 1) res = res / C * (C - 1);
    return res;
}

int main() {
    int T = 1; // case
    LL L;
    while (scanf("%lld", &L), L) {
        int d = 1;
        while (L % (d * 2) == 0 && d * 2 <= 8) d *= 2; // d = gcd(L, 8)
        LL C = 9 * L / d;
        LL phi = get_euler(C);

        LL res = 1e18;
        if (C % 2 == 0 || C % 5 == 0) res = 0; // 10 和 C 不互质无解（费马小定理）
        else {
            for (LL d = 1; d * d <= phi; d++) // 枚举 phi 的所有约数 d 和 phi / d
                if (phi % d == 0) {
                    if (qmi(10, d, C) == 1) res = min(res, d);
                    if (qmi(10, phi / d, C) == 1) res = min(res, phi / d);
                }
        }

        printf("Case %d: %lld\n", T++, res);
    }

    return 0;
}

如果不适用龟速乘，可以看编译器是否支持__int128，如果支持，可以省去龟速乘代码，直接在快速幂中修改，如下所示：

LL qmi(LL a, LL k, LL b) {
    LL res = 1;
    while (k) {
        if (k & 1) res = (__int128)res * a % b;
        a = (__int128) a * a % b;
        k >>= 1;
    }
    return res;
}