第五周专题（8.8-8.14）：数学

线性代数

A题轮状病毒（递推，DP，矩阵树定理）

这题是可以暴力打表找规律来求通项，或者硬推出 DP 方程，但是作为数学场的第一题，我们还是小小思考一下这题背后的数学性质：矩阵树定理。

定义 /(G/) 是一个 /(n/) 顶点的无向图，那么有度数矩阵 /(D(G)/) 为：

/[D_{ii}(G)=/deg(i),D_{ij}=0(i/not=j)
/]

同时定义 /(e(i,j)/) 为点 /(i,j/) 之间的边数，并定义邻接矩阵 /(A/) 为：（附：不允许自环）

/[A_{ij}(G)=A_{ji}(G)=e(i,j),i/not=j
/]

那么有 Laplace 矩阵（也被称为基尔霍夫矩阵）/(L(G)=D(G)-A(G)/)，而这张图的生成树的个数就等于该矩阵中任何一个 /(n-1/) 阶主子式的行列式值（/(k/) 阶主子式指从矩阵中选取 /(k/) 行 /(k/) 列组成的新的子矩阵，且选取的行列号必须相同，例如选了 1, 2, 7 行，那么列也必须选 1, 2, 7）。

对于本题，我们可以相对轻松的构造出基尔霍夫矩阵：

/[L=/begin{bmatrix}n & -1 & -1 & -1 & /cdots & -1 & -1
//-1 & 3 & -1 & 0 & /cdots & 0 & -1
//-1 & -1 & 3 & -1 & /cdots & 0 & 0
//-1 & 0 & -1 & 3 & /cdots & 0 & 0
// /vdots & /vdots & /vdots & /vdots & /ddots & /vdots & /vdots
//-1 & 0 & 0 & 0 & /cdots & 3 & -1
//-1 & -1 & 0 & 0 & /cdots & -1 & 3
/end{bmatrix}
/]

随后，我们仅需要求出该矩阵的一个 /(n-1/) 阶主子式的行列式即可（建议去掉第一行第一列）。

这题比较古老，所以必须手写高精度，所以我们直接上 python，但是我又不想写高斯消元，所以就直接套 /(O(n)/) 递推式了（具体证明可见OI-Wiki 矩阵树定理，就是手拆行列式后猜线性递推。

n = int(input())
dp = [0] * 110
dp[1], dp[2] = 3, 7
for i in range(3, n + 1):
    dp[i] = dp[i - 1] * 3 - dp[i - 2]
print(dp[n] - 2)

C题球形空间产生器（高斯消元）

在 /(n/) 维空间中存在着一个球体（2维的话就是一个圆），现在给定球面上的 /(n+1/) 个点的坐标，试求出球心坐标。

/(n/leq 10,|x|/leq 2*10^4/)，给定坐标的值均为六位浮点数，输出答案要求精确到三位

我们假定一个球面（球心坐标 /((a_1,a_2,/cdots,a_n)/)，半径为 /(R/)）的方程为

/[f(x_1,x_2,/cdots,x_n)=/sum/limits_{i=1}^n(x_i-a_i)^2=R^2
/]

这个方程看起来是一个 /(n/) 元二次方程，无法正常线性求解，实际上

通项公式

我们可以直接类似圆的求出该球的通项公式：/(/sum/limits_{i=1}^nx_i^2+/sum/limits_{i=1}^na_ix_i+c=0/)，一共 /(n+1/) 个未知项，坐标带入后求解线性方程组即可
作差法

我们假定 /(n=2/)，得到一个圆的方程 /((x-a)^2+(y-b)^2=R^2/)，那么圆上有两个点 /((x_1,y_1),(x_2,y_2)/)，可得

/[/begin{cases}
2x_1a+2y_1b=x_1^2+y_1^2-R^2+a^2+b^2
//
2x_2a+2y_2b=x_2^2+y_2^2-R^2+a^2+b^2
/end{cases}
/]

那么作差可得 /(2(x_1-x_2)a+2(y_1-y_2)b=(x_1^2+y_1^2)-(x_2^2+y_2^2)/)，也就是一个二元一次方程。

给定了三个点的话，就可以得到 /(C_3^2=3/) 个二元一次方程，而其中一个和另外两个线性相关（也就是说它可以通过另外两个方程来得到），所以实际上对于求解有效的仅有两个方程。

类似的，我们将其扩展到 /(n/) 维球上时，也发现 /(n+1/) 个点坐标可以推出 /((n+1)-1=n/) 个线性无关的方程组，足够求出球心了。

记第一个点的坐标为 /((x_1,x_2,/cdots,x_n)/)，第二个点的坐标为 /((y_1,y_2,/cdots,y_n)/)，那么可得方程组

/[2/sum/limits_{i=1}^n(x_i-y_i)a_i=/sum/limits_{i=1}^n(x_i^2-y_i^2)
/]

不管咋样，我们都可以求出这些方程组，然后扔给高斯消元来求解即可，复杂度 /(O(n^3)/)（我个人认为下面一个更好一些，因为求解出的答案不需要经过处理就可以直接输出）。

#include<bits/stdc++.h>
using namespace std;
const int N = 110;
const double eps = 1e-7;
int n;
double a[N][N], s[N][N], ans[N];
bool solve() {
   for (int i = 1; i <= n; ++i) {
       int r = i;
       for (int j = i + 1; j <= n; ++j)
           if (abs(a[r][i]) < abs(a[j][i])) r = j;
       if (abs(a[r][i]) < eps) return false;
       //C++11之后，可以这样直接交换两行
       //不过它复杂度不是O(n)而是O(N)，所以多组数据的时候还得手写
       if (i != r) swap(a[i], a[r]);

       double div = a[i][i];
       for (int j = i; j <= n + 1; ++j) a[i][j] /= div;
       for (int j = i + 1; j <= n; ++j) {
           div = a[j][i];
           for (int k = i; k <= n + 1; ++k)
               a[j][k] -= a[i][k] * div;
       }
   }
   for (int i = n; i >= 1; i--) {
       ans[i] = a[i][n + 1];
       for (int j = i + 1; j <= n; ++j)
           ans[i] -= a[i][j] * ans[j];
   }
   return true;
}
int main()
{
   cin >> n;
   for (int i = 0; i <= n; ++i)
       for (int j = 1; j <= n; ++j)
           cin >> s[i][j];
   for (int i = 1; i <= n; ++i)
       for (int j = 1; j <= n; ++j) {
           a[i][j] = 2 * (s[i][j] - s[0][j]);
           a[i][n + 1] += s[i][j] * s[i][j] - s[0][j] * s[0][j];
       }
   solve();
   for (int i = 1; i <= n; ++i)
       printf("%.3f ", ans[i]);
   return 0;
}

J题可乐（邻接矩阵，矩阵快速幂）

给定一张 /(n/) 点 /(m/) 边的无向图，现在有一个机器人在点 /(1/) 处，每秒钟他都可以选择下面三条指令中的一个进行操作：

在原地不动

走向另外一个相邻城市

自爆（自爆后无法进行其他操作）

现在时刻为 /(0/)，那么机器人可以在时刻 /(1,2,/cdots,t/) 处各进行一次操作，问该机器人的总行动方案数是多少？（在一个点走向不同城市，虽然都是操作 2，但是属于不同的行动）

/(t/leq 10^6,n/leq 30,m/leq 100/)

我们假定没有操作 3，那么根据离散数学的知识，我们直接拿出邻接矩阵 /(E/)，那么 /(E^1/) 中的 /(i/) 行 /(j/) 列就代表了从 /(i/) 到 /(j/)，经过时刻 1 的总方案数，/(E^2/) 则代表经过时刻 2 后的方案数，以此类推。（不太好想到，但是不难证明）

考虑到自爆，我们额外建立一个点，所有点向其连边，但是这个点不向其他点连出边（自己除外）。

得到邻接矩阵后，我们直接矩阵快速幂求 /(E^t/) 即可，复杂度 /(O(n^3/log t)/)。

#include<bits/stdc++.h>
using namespace std;
#define LL long long
const int N = 40;
const LL mod = 2017;
struct Matrix {
    int n;
    LL a[N][N];
    void init(int _n, int opt) {
        memset(a, 0, sizeof(a));
            n = _n;
            if (opt == 1)
                for (int i = 1; i <= n; ++i)
                    a[i][i] = 1;
    }
    Matrix operator = (const Matrix &B) {
        this->n = B.n;
        for (int i = 1; i <= n; ++i)
            for (int j = 1; j <= n; ++j)
                this->a[i][j] = B.a[i][j];
        return *this;
    }
    Matrix operator * (const Matrix &B) const {
        Matrix res;
        res.init(n, 0);
        for (int i = 1; i <= n; ++i)
            for (int j = 1; j <= n; ++j)
                for (int k = 1; k <= n; ++k) {
                    res.a[i][j] += (this->a[i][k]) * B.a[k][j];
                    res.a[i][j] %= mod;
                }
        return res;
    }
};
Matrix quickpow(Matrix P, LL n) {
    if (n == 1) return P;
    Matrix res;
    res.init(P.n, 1);
    while (n) {
        if (n & 1) res = res * P;
        n >>= 1;
        P = P * P;
    }
    return res;
}
int main()
{
    Matrix A;
    int n, m, t;
    scanf("%d%d", &n, &m);
    A.init(n + 1, 1);
    for (int i = 1; i <= n; ++i)
        A.a[i][n + 1] = 1;
    for (int i = 1, u, v; i <= m; ++i) {
        scanf("%d%d", &u, &v);
        A.a[u][v] = A.a[v][u] = 1;
    }
    scanf("%d", &t);
    Matrix D = quickpow(A, t);
    LL ans = 0;
    for (int i = 1; i <= n + 1; ++i)
        ans = (ans + D.a[1][i]) % mod;
    printf("%lld", ans);
    return 0;
}

组合数学

E题排队（组合数学，高精度）

现在有 /(n/) 名男生，/(m/) 名女生和 2 名老师，要让他们站成一列，并且任意两名女生不能相邻，且两名老师也不能相邻，问一共有多少种排队方案？（任意两人都是不同的）

/(n,m/leq 2*10^3/)

我们直接先排老师和男生，后面在空位上站女生，考虑两种情况：

两个老师相邻，等着中间站一个女生

那么仅考虑老师和男生的情况下，一共有 /(A_2^2A_{n+1}^{n+1}/) 种排列方法，随后选一个女生放在中间，剩下来的空位填上剩下来的 /(m-1/) 个女生，总计 /(A_2^2A_{n+1}^{n+1}*m*A_{n+2}^{m-1}/) 种方案（注意这里必须要乘上 /(m/)，因为选女生有 /(m/) 种方案，样例 /(n=m=1/) 太弱了，没查出来，很烦）
两个老师中间隔着男生

排列的总方案是 /(A_{n+2}^{n+2}-A_2^2A_{n+1}^{n+1}=n*A_{n+1}^{n+1}/) 种，然后再乘一个 /(A_{n+3}^m/) 即可

综上相加，合并一下就可以得到答案

/[A_{n+1}^{n+1}(2mA_{n+2}^{m-1}+nA_{n+3}^{m})
/]

注意一些无解的情况，也就是当 /(A_x^y/) 中出现了 /(y<0/) 或者 /(y>x/) 的时候就返回 0。

这题是一个老题，所以必须得高精度写，我图省事，直接 python 莽过去了（python 递归层次有限，只能递推求阶乘了）

def A(n, m, f):
    if m > n or m < 0:
        return 0
    return f[n] // f[n - m]

f = [0] * 2010
f[0] = 1
for i in range(1, 2010):
    f[i] = f[i - 1] * i
n, m = map(int, input().split())
ans = f[n + 1] * (2 * m * A(n + 2, m - 1, f) + n * A(n + 3, m, f))
print(ans)

K题 combination（组合数，Lucas定理）

给定 /(T/) 组数据，每组数据中需要求出 /(C_{n}^m/bmod 10007/) 的值。

/(T/leq 200,1/leq m/leq n/leq 2*10^8/)

对于这种 /(n,m/) 极大的组合数，我们直接使用 Lucas 定理即可：

对于给定数字 /(m=sp+q,n=tp+r(0/leq q,r <p)/)，那么有

/[C_n^m/bmod p=C_t^sC_r^q/bmod p
/]

（对于递归流程中出现的 /(q>r/) 的情况，返回 0 即可，这说明整体的计算结果是 /(p/) 的倍数）

那么，我们现在只需要考虑 /(n,m/leq 10^4/) 情况下的求解即可。

我们知道，/(C_n^m=/frac{n!}{m!(n-m)!}/)，我们可以预处理后 /(O(1)/) 得到 /(n!/)，那么我们考虑怎么 /(O(1)/) 得到 /(n!/) 的阶乘（当然，时限不紧张的时候，直接每次都快速幂+费马小定理即可）。

我们知道 /(/frac{1}{(k-1)!}=/frac{k}{k!}/)，所以我们可以先快速幂求出 /({n!}/) 的逆元，随后就可以一直逆推得到 /((n-1)!,(n-2)!,/cdots,2!,1!/) 的逆元，预处理完后直接 /(O(1)/) 查询即可。

本题预处理时间为 /(O(10007)/)，随后每组数据的求解都是接近 /(O(1)/) 的（严格来说是 /(O(/log_{10007}^{n})/)）。

#include<bits/stdc++.h>
using namespace std;
#define LL long long
const int N = 10010;
const LL P = 10007;
LL f[N];
LL power(LL a, LL b) {
    LL res = 1;
    while (b) {
        if (b & 1) res = res * a % P;
        b >>= 1;
        a = a * a % P;
    }
    return res;
}
LL inv(LL x) { return power(x, P - 2); }
LL C(LL n, LL m) {
    if (m > n) return 0;
    return f[n] * inv(f[m]) % P * inv(f[n - m]) % P;
}
LL Lucas(LL n, LL m) {
    if (!m) return 1;
    return Lucas(n / P, m / P) * C(n % P, m % P) % P;
}
int main()
{
    f[0] = 1;
    for (int i = 1; i < N; ++i)
        f[i] = f[i - 1] * i % P;
    int T;
    cin >> T;
    while (T--) {
        LL n, m;
        cin >> n >> m;
        cout << Lucas(n, m) << endl;
    }
    return 0;
}

数论

F题模积和（取模操作的转换，整数分块）

计算 /(/sum/limits_{i=1}^n/sum/limits_{j=1}^m(n /bmod i)(m/bmod j)(i/not=j)/) 的值。

/(n,m/leq 10^9/)

我们假定 /(n<m/)，那显然有

/[/sum/limits_{i=1}^n/sum/limits_{j=1}^m(n /bmod i)(m/bmod j)(i/not=j)
//
=/sum/limits_{i=1}^n/sum/limits_{j=1}^m(n /bmod i)(m/bmod j)-/sum/limits_{i=1}^n(n/bmod i)(m/bmod i)
//
=(/sum/limits_{i=1}^n(n/bmod i))(/sum/limits_{i=1}^m(m/bmod i))-/sum/limits_{i=1}^n(n/bmod i)(m/bmod i)
/]

我看到这个式子的时候一眼想到的是整数分块，但是没想起来具体写法，看到题解后才想起来咋搞：

/[n/bmod i=n – /lfloor/frac{n}{i}/rfloor*i
/]

我们来看一下前一个式子：

/[/sum/limits_{i=1}^n(n/bmod i))=/sum/limits_{i=1}^n(n – /lfloor/frac{n}{i}/rfloor*i)=n^2-/sum/limits_{i=1}^n i*/lfloor/frac{n}{i}/rfloor
/]

对于后面那个整除式子的处理，我们可以套一个数论分块来处理。

数论分块主要用于解决 /(/sum/limits_{i=1}^nf(i)g(/lfloor/frac{n}{i}/rfloor)/) 的求解问题，其本质在于通过 /(/lfloor/frac{n}{i}/rfloor/) 的本质数量的数量级是 /(O(/sqrt{n})/) 的性质来减少运算量。若 /(/lfloor/frac{n}{l}/rfloor=x/)，那么最大的满足 /(/lfloor/frac{n}{r}/rfloor=x/) 的 /(r/) 的值就是 /(/lfloor/frac{n}{x}/rfloor/)。根据此，我们得到了这些值隶属的不同区间，而前面的部分就直接用公式推导或者前缀和的形式来处理，从而 /(O(/sqrt{n})/) 的求出表达式的值。（有时候分块的时候不会数学推导，没关系，可以直接二分这个 /(r/)）。

我们再来看一下后一个式子：

/[/sum/limits_{i=1}^n(n/bmod i)(m/bmod i)=/sum/limits_{i=1}^n(n – /lfloor/frac{n}{i}/rfloor*i)(m – /lfloor/frac{m}{i}/rfloor*i)
//=n^2m+/sum/limits_{i=1}^ni^2/lfloor/frac{n}{i}/rfloor/lfloor/frac{m}{i}/rfloor-m/sum/limits_{i=1}^n/lfloor/frac{n}{i}/rfloor*i-n/sum/limits_{i=1}^n/lfloor/frac{m}{i}/rfloor*i
/]

后面两个减式不需要额外计算，我们只考虑中间的这个带 /(i^2/) 的式子咋处理。额，实际上，这题我们也可以直接用类似上面整数分块的方式来做，复杂度大概是在 /(O(n^/frac{2}{3})/) 这样。

那么，我们记 /(A(x,y)=/sum/limits_{i=1}^x i*/lfloor/dfrac{y}{i}/rfloor,B(x,y)=/sum/limits_{i=1}^x i^2/lfloor/dfrac{x}{i}/rfloor/lfloor/dfrac{y}{i}/rfloor/)，那么本题答案即为

/[(n^2-A(n,n))*(m^2-A(m,m))+mA(n,n)+nA(n,m)-B(n,m)-n^2m
/]

#include<bits/stdc++.h>
using namespace std;
#define LL long long
const LL mod = 19940417;
//手算2和3的逆元，省的写快速幂
LL f(LL n) { return n * (n + 1) % mod * 9970209 % mod; }
LL g(LL n) { return f(n) * (2 * n + 1) % mod * 6646806 % mod; }
LL A(LL n, LL m) {
    LL res = 0;
    for (int l = 1, r; l <= n; l = r + 1) {
        r = min(m / (m / l), n);
        res += ((f(r) - f(l - 1) + mod) % mod) * (m / l) % mod;
        res %= mod;
    }
    return res;
}
LL B(LL n, LL m) {
    LL res = 0;
    for (int l = 1, r; l <= n; l = r + 1) {
        r = min(n, min(n / (n / l), m / (m / l)));
        res += ((g(r) - g(l - 1) + mod) % mod) * (n / l) % mod * (m / l) % mod;
        res %= mod;
    }
    return res;
}
int main()
{
    LL n, m;
    cin >> n >> m;
    if (n > m) swap(n, m);
    //勤取模
    LL ans = (n * n % mod - A(n, n) + mod) % mod;
    ans = ans * ((m * m % mod - A(m, m) + mod) % mod) % mod;
    ans = (ans + m * A(n, n)) % mod;
    ans = (ans + n * A(n, m)) % mod;
    ans = (ans - (n * n % mod * m % mod) + mod) % mod;
    ans = (ans - B(n, m) % mod + mod) % mod;
    cout << ans << endl;
    return 0;
}

G题古代猪文（欧拉降幂，Lucas定理，中国剩余定理）

计算下面这个表达式的值，保证 /(1/leq n,g/leq 10^9/)。

/[g^s/bmod 999911659(s=/sum/limits_{d|n}C_n^d)
/]

/(s/) 的值不可能用 long long 存的下，只能强制欧拉降幂搞一波了。

考虑到 999911659 是一个质数，所以只要不是 /(g=999911659/)（这种情况直接输出 0 就完事了），那模数和 /(g/) 就是互质的，根据欧拉定理（当 /(/gcd(a,p)=1/) 时，/(a^{/phi(p)}/equiv1(/bmod p)/)），我们只需要求出 /(s/bmod /phi(999911659)/) 即可。

当 /(p/) 是质数时，/(/phi(p)=p-1/)（此时欧拉定理就变成了费马小定理），而 /(999911658/) 这个数对于 Lucas 来说也显得过于大了，需要考虑如何优化求解。

进行质因数分解，我们发现 /(999911658=2*3*4679*35617/)，那我们可以用 Lucas 来依次算出 /(C_n^d/) 对这四个数取模的结果，得到解 /(a_1,a_2,a_3,a_4/)，那么就有

/[/begin{cases}
x/equiv a_1/pmod 2
//
x/equiv a_2/pmod 3
//
x/equiv a_3/pmod {4679}
//
x/equiv a_4/pmod {35617}
//
/end{cases}
/]

而解 /(x/) 的流程就是中国剩余定理了。

#include<bits/stdc++.h>
using namespace std;
#define LL long long
const LL mod = 999911659;
const LL P[4] = {2, 3, 4679, 35617};
LL power(LL a, LL b, LL p) {
    LL res = 1;
    while (b) {
        if (b & 1) res = res * a % p;
        a = a * a % p;
        b >>= 1;
    }
    return res;
}
LL inv(LL x, LL p) { return power(x, p - 2, p); }
const int N = 40010;
LL f[N][4];
LL C(LL n, LL m, int id) {
    if (m > n) return 0;
    return f[n][id] * inv(f[m][id], P[id]) % P[id] * inv(f[n - m][id], P[id]) % P[id];
}
LL Lucas(LL n, LL m, int id) {
    if (m == 0) return 1;
    return Lucas(n / P[id], m / P[id], id) * C(n % P[id], m % P[id], id) % P[id];
}
//
LL solve(LL n, LL d) {
    LL a[4], m[4], t[4], M = mod - 1, res = 0;
    for (int i = 0; i < 4; ++i) {
        a[i] = Lucas(n, d, i), m[i] = M / P[i], t[i] = inv(m[i], P[i]);
        res += a[i] * t[i] * m[i];
    }
    return res % M;
}
int main()
{
    //init
    for (int i = 0; i < 4; ++i)
        f[0][i] = 1;
    for (int i = 1; i < N; ++i)
        for (int j = 0; j < 4; ++j)
            f[i][j] = f[i - 1][j] * i % P[j];
    //read
    LL n, g;
    cin >> n >> g;
    //special
    g %= mod;
    if (g == 0) {
        cout << 0 << endl;
        return 0;
    }
    //solve
    vector<LL> ds;
    for (LL d = 1; d * d <= n; ++d)
        if (n % d == 0) {
            ds.push_back(d);
            if (d * d != n) ds.push_back(n / d);
        }
    LL s = 0;
    for (LL d : ds)
        s = (s + solve(n, d)) % (mod - 1);
    cout << power(g, s, mod) << endl;
    return 0;
}

I题计算器（BSGS算法）

给定 /(T/) 次询问，每次询问给定 /(y,z,p/)，需要进行下面三种操作之一（本题比较特殊，对于同一组数据，里面所有询问的操作类型是相同的，即操作类型在发起询问前给出，所以代码的运行流程中只会进行一种计算）。

计算 /(y^z/bmod p/)

求出满足 /(xy/equiv z(/bmod p)/) 的最小非负整数 /(x/)

求出满足 /(y^x/equiv z(/bmod p)/) 的最小非负整数 /(x/)

/(1/leq y,z,p/leq 10^9/)，/(p/) 是质数，/(1/leq T/leq 10/)，无解时输出 Orz, I cannot find x!

操作 1 就是快速幂，写一下就行了。

操作 2 是求逆元，最保险的方式是扩展欧几里得求出通解后慢慢调，不过因为 /(p/) 是质数，所以我们可以偷懒一下，直接费马小定理搞起来，但是需要注意一点细节：当 /(y/) 是 /(p/) 的倍数时，逆元并不是很存在（肉眼也能看出来端倪吧），如果 /(z/) 也是 /(p/) 的倍数就输出 0，否则输出无解。

对于操作 3，这就是 BSGS 算法了，大家可以去洛谷的模板题上学一下（当然，这题也有亿些小细节：当 /(y/) 是 /(p/) 倍数时，若 /(z/) 是 /(p/) 的倍数则输出 1（注意了，是 1，因为 /(y^0=1/not=0/)），不是倍数的话则需要检查一下对 /(p/) 取模是不是 1，是的话直接返回 0 即可）。

~~讲道理，这题 debug 的时间要比我学 BSGS 的时间要长。~~

#include<bits/stdc++.h>
using namespace std;
#define LL long long
LL power(LL a, LL b, LL p) {
    LL res = 1;
    while (b) {
        if (b & 1) res = res * a % p;
        a = a * a % p;
        b >>= 1;
    }
    return res;
}
LL BSGS(LL a, LL b, LL p) {
    if (b == 1) return 0;
    LL t = sqrt(p) + 1;
    map<LL, LL> vis;
    //因为x = At-B，所以尽量保留B大的
    for (int i = 0; i < t; ++i)
        vis[b * power(a, i, p) % p] = i;
    LL aa = power(a, t, p);
    for (int A = 1; A <= t; ++A) {
        LL x = power(aa, A, p);
        if (vis.find(x) != vis.end()) {
            LL B = vis[x];
            return A * t - B;
        }
    }
    //遍历了整个解空间都找不到解，返回-1
    return -1;
}
int main()
{
    LL T, K;
    cin >> T >> K;
    while (T--) {
        LL a, b, p;
        cin >> a >> b >> p;
        if (K == 1) cout << power(a, b, p) << endl;
        else if (K == 2) {
            b %= p;
            if (a % p == 0) {
                if (b)
                    cout << "Orz, I cannot find x!" << endl;
                else
                    cout << 0 << endl;
            }
            else cout << b * power(a, p - 2, p) % p << endl;
        }
        else {
            b %= p;
            if (a % p == 0) {
                if (b == 0)
                    cout << 1 << endl;
                else if (b == 1)
                    cout << 0 << endl;
                else
                    cout << "Orz, I cannot find x!" << endl;
            }
            else {
                LL res = BSGS(a, b, p);
                if (res == -1) cout << "Orz, I cannot find x!" << endl;
                else cout << res << endl;
            }
        }
    }
}

原创文章，作者：ItWorker，如若转载，请注明出处：https://blog.ytso.com/281980.html

第五周专题（8.8-8.14）：数学（8/15）

第五周专题（8.8-8.14）：数学

线性代数

A题轮状病毒（递推，DP，矩阵树定理）

C题球形空间产生器（高斯消元）

J题可乐（邻接矩阵，矩阵快速幂）

组合数学

E题排队（组合数学，高精度）

K题 combination（组合数，Lucas定理）

数论

F题模积和（取模操作的转换，整数分块）

G题古代猪文（欧拉降幂，Lucas定理，中国剩余定理）

I题计算器（BSGS算法）

发表回复

第五周专题（8.8-8.14）：数学（8/15）

第五周专题（8.8-8.14）：数学

线性代数

A题 轮状病毒（递推，DP，矩阵树定理）

C题 球形空间产生器（高斯消元）

J题 可乐（邻接矩阵，矩阵快速幂）

组合数学

E题 排队（组合数学，高精度）

K题 combination（组合数，Lucas定理）

数论

F题 模积和（取模操作的转换，整数分块）

G题 古代猪文（欧拉降幂，Lucas定理，中国剩余定理）

I题 计算器（BSGS算法）

相关推荐

发表回复

A题轮状病毒（递推，DP，矩阵树定理）

C题球形空间产生器（高斯消元）

J题可乐（邻接矩阵，矩阵快速幂）

E题排队（组合数学，高精度）

F题模积和（取模操作的转换，整数分块）

G题古代猪文（欧拉降幂，Lucas定理，中国剩余定理）

I题计算器（BSGS算法）