第五周专题(8.8-8.14):数学
比赛链接
线性代数
A题 轮状病毒(递推,DP,矩阵树定理)
这题是可以暴力打表找规律来求通项,或者硬推出 DP 方程,但是作为数学场的第一题,我们还是小小思考一下这题背后的数学性质:矩阵树定理。
定义 (G) 是一个 (n) 顶点的无向图,那么有度数矩阵 (D(G)) 为:
同时定义 (e(i,j)) 为点 (i,j) 之间的边数,并定义邻接矩阵 (A) 为:(附:不允许自环)
那么有 Laplace 矩阵(也被称为基尔霍夫矩阵)(L(G)=D(G)-A(G)),而这张图的生成树的个数就等于该矩阵中任何一个 (n-1) 阶主子式的行列式值((k) 阶主子式指从矩阵中选取 (k) 行 (k) 列组成的新的子矩阵,且选取的行列号必须相同,例如选了 1, 2, 7 行,那么列也必须选 1, 2, 7)。
对于本题,我们可以相对轻松的构造出基尔霍夫矩阵:
随后,我们仅需要求出该矩阵的一个 (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) 个点的坐标,试求出球心坐标。
(nleq 10,|x|leq 2*10^4),给定坐标的值均为六位浮点数,输出答案要求精确到三位
我们假定一个球面(球心坐标 ((a_1,a_2,cdots,a_n)),半径为 (R))的方程为
这个方程看起来是一个 (n) 元二次方程,无法正常线性求解,实际上
-
通项公式
我们可以直接类似圆的求出该球的通项公式:(sumlimits_{i=1}^nx_i^2+sumlimits_{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)),那么可得方程组
[2sumlimits_{i=1}^n(x_i-y_i)a_i=sumlimits_{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,但是属于不同的行动)
(tleq 10^6,nleq 30,mleq 100)
我们假定没有操作 3,那么根据离散数学的知识,我们直接拿出邻接矩阵 (E),那么 (E^1) 中的 (i) 行 (j) 列就代表了从 (i) 到 (j),经过时刻 1 的总方案数,(E^2) 则代表经过时刻 2 后的方案数,以此类推。(不太好想到,但是不难证明)
考虑到自爆,我们额外建立一个点,所有点向其连边,但是这个点不向其他点连出边(自己除外)。
得到邻接矩阵后,我们直接矩阵快速幂求 (E^t) 即可,复杂度 (O(n^3log 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,mleq 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_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}^mbmod 10007) 的值。
(Tleq 200,1leq mleq nleq 2*10^8)
对于这种 (n,m) 极大的组合数,我们直接使用 Lucas 定理即可:
对于给定数字 (m=sp+q,n=tp+r(0leq q,r <p)),那么有
[C_n^mbmod p=C_t^sC_r^qbmod p ]
(对于递归流程中出现的 (q>r) 的情况,返回 0 即可,这说明整体的计算结果是 (p) 的倍数)
那么,我们现在只需要考虑 (n,mleq 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题 模积和(取模操作的转换,整数分块)
计算 (sumlimits_{i=1}^nsumlimits_{j=1}^m(n bmod i)(mbmod j)(inot=j)) 的值。
(n,mleq 10^9)
我们假定 (n<m),那显然有
我看到这个式子的时候一眼想到的是整数分块,但是没想起来具体写法,看到题解后才想起来咋搞:
我们来看一下前一个式子:
对于后面那个整除式子的处理,我们可以套一个数论分块来处理。
数论分块主要用于解决 (sumlimits_{i=1}^nf(i)g(lfloorfrac{n}{i}rfloor)) 的求解问题,其本质在于通过 (lfloorfrac{n}{i}rfloor) 的本质数量的数量级是 (O(sqrt{n})) 的性质来减少运算量。若 (lfloorfrac{n}{l}rfloor=x),那么最大的满足 (lfloorfrac{n}{r}rfloor=x) 的 (r) 的值就是 (lfloorfrac{n}{x}rfloor)。根据此,我们得到了这些值隶属的不同区间,而前面的部分就直接用公式推导或者前缀和的形式来处理,从而 (O(sqrt{n})) 的求出表达式的值。(有时候分块的时候不会数学推导,没关系,可以直接二分这个 (r))。
我们再来看一下后一个式子:
后面两个减式不需要额外计算,我们只考虑中间的这个带 (i^2) 的式子咋处理。额,实际上,这题我们也可以直接用类似上面整数分块的方式来做,复杂度大概是在 (O(n^frac{2}{3})) 这样。
那么,我们记 (A(x,y)=sumlimits_{i=1}^x i*lfloordfrac{y}{i}rfloor,B(x,y)=sumlimits_{i=1}^x i^2lfloordfrac{x}{i}rfloorlfloordfrac{y}{i}rfloor),那么本题答案即为
#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定理,中国剩余定理)
计算下面这个表达式的值,保证 (1leq n,gleq 10^9)。
[g^sbmod 999911659(s=sumlimits_{d|n}C_n^d) ]
(s) 的值不可能用 long long 存的下,只能强制欧拉降幂搞一波了。
考虑到 999911659 是一个质数,所以只要不是 (g=999911659)(这种情况直接输出 0 就完事了),那模数和 (g) 就是互质的,根据欧拉定理(当 (gcd(a,p)=1) 时,(a^{phi(p)}equiv1(bmod p))),我们只需要求出 (sbmod 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),那么就有
而解 (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^zbmod p)
- 求出满足 (xyequiv z(bmod p)) 的最小非负整数 (x)
- 求出满足 (y^xequiv z(bmod p)) 的最小非负整数 (x)
(1leq y,z,pleq 10^9),(p) 是质数,(1leq Tleq 10),无解时输出
Orz, I cannot find x!
操作 1 就是快速幂,写一下就行了。
操作 2 是求逆元,最保险的方式是扩展欧几里得求出通解后慢慢调,不过因为 (p) 是质数,所以我们可以偷懒一下,直接费马小定理搞起来,但是需要注意一点细节:当 (y) 是 (p) 的倍数时,逆元并不是很存在(肉眼也能看出来端倪吧),如果 (z) 也是 (p) 的倍数就输出 0,否则输出无解。
对于操作 3,这就是 BSGS 算法了,大家可以去洛谷的模板题上学一下(当然,这题也有亿些小细节:当 (y) 是 (p) 倍数时,若 (z) 是 (p) 的倍数则输出 1(注意了,是 1,因为 (y^0=1not=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;
}
}
}
}