当前位置：首页 > article >正文

【矩阵快速幂】P6601 「EZEC-2」机器|普及+

article 2026/2/24 16:14:11

本文涉及知识点

【矩阵快速幂】封装类及测试用例及样例

P6601 「EZEC-2」机器

题目背景

tlx 喜欢科幻小说。

小宇宙中只剩下漂流瓶和生态球。漂流瓶隐没于黑暗里,在一千米见方的宇宙中,只有生态球里的小太阳发出一点光芒。在这个小小的生命世界中,几只清澈的水球在零重力环境中静静地飘浮着,有一条小鱼从一只水球中蹦出,跃入另一只水球,轻盈地穿游于绿藻之间。在一小块陆地上的草丛中,有一滴露珠从一片草叶上脱离,旋转着飘起,向太空中折射出一缕晶莹的阳光。
$\qquad \qquad \qquad \qquad \qquad \qquad \qquad\qquad\qquad\qquad\qquad\qquad\qquad --$ 《三体》

在另一个宇宙，将是另一番奇景吧。

在那里，重力似乎变得微不足道了，引力机器成了司空见惯的东西。

引力机器装置内并没有重力，即若有物体在机器上运动，运动过程中只受机器给予的引力，这个力有一定几率使物体向施力物体快速移动，达到一定动力时就可以实现瞬移。

题目描述

一个引力机器由一个光滑圆轨道和 $2 n$ 个小孔组成（小孔按逆时针从 $1$ 到 $2 n$ 编号，每两个相邻的小孔所夹的劣弧度数为 $\dfrac{\pi}{n}$ ），每个小孔与和其夹角为 $\pi$ 的另一个小孔有通道相连，比如当 $n = 2$ 时， $1$ 号孔和 $3$ 号孔相连。

当 $n = 2$ 时，这个装置的构造大概是这样的：

现在我们在 $1$ 号孔处放一个小球，使它一直沿逆时针方向做匀速圆周运动，在不瞬移的情况下，每一秒恰好能从一个小孔运动至下一个小孔。

由于未来实验室构造奇特（内部的引力提供装置太神了！），每经过一个小孔时，有 $p$ 的概率立刻瞬移（即不花费时间）到通道对面的小孔并继续沿逆时针方向做匀速圆周运动，也就是有 $1 - p$ 的概率继续沿圆周向下一个小孔运动。

值得注意的是，每一单位时刻，小球只能瞬移一次。

简单地说，若某一时刻小球在小孔 $i$ ，则下一时刻它可能运动到小孔 $\bmod 2n + 1$ 或 $\bmod 2n + 1$ ，概率分别为 $1 - p$ 和 $p$ 。

现在 tlx 有两个一模一样的引力机器，两个小球同时从 $1$ 号孔开始运动。他会随机（所有可能选择的概率相同）选择一个二元组 $1\leqslant i\leqslant 2n,0\leqslant j\leqslant t,i,j\in \mathbb Z$ ) 分别代表小孔编号和时间，你需要求出时间为 $j$ 时两个引力机器的小孔 $i$ 同时有小球停留（运动经过小孔但瞬移到对面了不算停留） 的概率。

注意：小球刚开始运动时也可能瞬移到对面的小孔。

为方便计算，我们规定：所有概率都是在模 $10^9+7$ 意义下的。

输入格式

输入数据共一行，三个整数 $n, p, t$ ，分别代表引力装置的小孔数的一半，瞬移的概率对 $10^9+7$ 取模的结果，选择的时间的范围的上界。

输出格式

共一行，一个整数，代表两个小球同时经过所选位置的概率对 $10^9+7$ 取模的结果。

输入输出样例 #1

输入 #1

2 500000004 1

输出 #1

125000001

输入输出样例 #2

输入 #2

6 114514 11

输出 #2

756497239

说明/提示

【数据范围与约定】

本题采用捆绑测试。

具体计分方式如下：

Subtask $1$ ( $7$ points)：满足 $p\in \{0,1\}$ ；
Subtask $2$ ( $13$ points)：满足 $t\leqslant 20,n\leqslant50$ ；
Subtask $3$ ( $20$ points)：满足 $t\leqslant 10^3,n\leqslant50$ ；
Subtask $4$ ( $10$ points)：满足 $t\leqslant 10^3$ ；
Subtask $5$ ( $10$ points)：满足 $t\leqslant 10^6$ ；
Subtask $6$ ( $15$ points)：满足 $n\leqslant50$ ；
Subtask $7$ ( $25$ points)：无特殊限制。

对于 $100\%$ 的数据，满足 $2\leqslant n\leqslant 500$ ， $0\leqslant p\leqslant 10^9+6$ ， $0\leqslant t \leqslant 10^9$ 。

注意：不做说明的数据范围即为极限数据范围。

【样例解释 #1】

$500000004$ 是模 $10^9+7$ 意义下的 $\dfrac{1}{2}$ 。

下面为了方便，记 $P (i, j)$ 为选择的二元组为 $(i, j)$ 时的概率。

所有概率不为 $0$ 的二元组有：
$P(1,0)=\dfrac{1}{4},P(3,0)=\dfrac{1}{4},P(2,1)=\dfrac{1}{4},P(4,1)=\dfrac{1}{4}$ 。

所有可以选择的二元组有：
$(1, 0), (1, 1), (2, 0), (2, 1), (3, 0), (3, 1), (4, 0), (4, 1)$ ，共 $8$ 种。

所以总的概率：
$P=\dfrac{1}{8}×\dfrac{1}{4}×4+\dfrac{1}{8}×0×4=\dfrac{1}{8}$

在模 $10^9+7$ 意义下为 $125000001$ ，即为输出的答案。

【其他提示】

如果你不了解分数取模，可以查看这里。
如果你不明白题目中角度的表示方法，可以查看弧度制。

矩阵快速幂逆元

小数（分数）取模 a/b = a*(b的逆元）
下标改成从0开始,故初始在0，目标是j-1。
性质一：瞬移两次和瞬移0次是一样的。故瞬移x次相当于瞬移x%2次。
结论：瞬移0次，处于t%(2n)；瞬移一次处于(t+n)%2n。
x记录瞬移偶数次的概率，初始1-p。y记录瞬移奇数次的概率，初始p。我们通过二元一次方程组迭代，即求x1=ax+by,y1=cx+dy的a,b,c,d。
x：有1-p的概率不瞬移，故a=1-p；有p的概率瞬移，故c=p。
y：有1-p的概率不瞬移，d=1-p；b=p。

t秒时，瞬移偶数次的概率是f0,奇数次的概率是f1。两台机器同时瞬移偶数次的概率为g0=f0f0,奇数次的概率为g1=f1f1。
我们假定已经选中(i,t),i是未知，t是已知，则:
如果i是0次瞬移的位置，则 g0的概率同时有球。2n分之一的几率选中i，故几率为:g0/2n。
如果i是1次瞬移的位置，则g1的概率同时有球，2n分之一的几率选中i，故几率为:g1/2n。
i是其它位置，有球的概率为0。
故选中(?,t)的情况下：有 (g0+g1)/2n,令h(t) = (g1+j2)
选中(?,t1)的概率为：1/(1+T)
则结果为： $\sum (g0+g1)$ /(T+1)/(2n)
我通过{a,b,c,d}可以求出求xx,xy,yy的转移矩阵，我们要求其和，详见下文。故需要4阶，最后一个元素记录和。
pre = {xx,xy,yy,xx+yy}
4*4矩阵，初始为求{xx,xy,yy}，然后：
mat[3][3]=1 h(0…t-1)之和。
mat[0…2][3] += mat[0…2][0] xx之和。
mat[0…2][3] += mat[0…2][2] yy之和。
pre.back()/(2n)/(T+1)便是答案。
时间复杂度：O(logt)。

代码

核心代码

#include <iostream>
#include <sstream>
#include <vector>
#include<map>
#include<unordered_map>
#include<set>
#include<unordered_set>
#include<string>
#include<algorithm>
#include<functional>
#include<queue>
#include <stack>
#include<iomanip>
#include<numeric>
#include <math.h>
#include <climits>
#include<assert.h>
#include<cstring>
#include<list>#include <bitset>
using namespace std;template<class T1, class T2>
std::istream& operator >> (std::istream& in, pair<T1, T2>& pr) {in >> pr.first >> pr.second;return in;
}template<class T1, class T2, class T3 >
std::istream& operator >> (std::istream& in, tuple<T1, T2, T3>& t) {in >> get<0>(t) >> get<1>(t) >> get<2>(t);return in;
}template<class T1, class T2, class T3, class T4 >
std::istream& operator >> (std::istream& in, tuple<T1, T2, T3, T4>& t) {in >> get<0>(t) >> get<1>(t) >> get<2>(t) >> get<3>(t);return in;
}template<class T = int>
vector<T> Read() {int n;scanf("%d", &n);vector<T> ret(n);for (int i = 0; i < n; i++) {cin >> ret[i];}return ret;
}template<class T = int>
vector<T> Read(int n) {vector<T> ret(n);for (int i = 0; i < n; i++) {cin >> ret[i];}return ret;
}template<int N = 1'000'000>
class COutBuff
{
public:COutBuff() {m_p = puffer;}template<class T>void write(T x) {int num[28], sp = 0;if (x < 0)*m_p++ = '-', x = -x;if (!x)*m_p++ = 48;while (x)num[++sp] = x % 10, x /= 10;while (sp)*m_p++ = num[sp--] + 48;AuotToFile();}void writestr(const char* sz) {strcpy(m_p, sz);m_p += strlen(sz);AuotToFile();}inline void write(char ch){*m_p++ = ch;AuotToFile();}inline void ToFile() {fwrite(puffer, 1, m_p - puffer, stdout);m_p = puffer;}~COutBuff() {ToFile();}
private:inline void AuotToFile() {if (m_p - puffer > N - 100) {ToFile();}}char  puffer[N], * m_p;
};template<int N = 1'000'000>
class CInBuff
{
public:inline CInBuff() {}inline CInBuff<N>& operator>>(char& ch) {FileToBuf();ch = *S++;return *this;}inline CInBuff<N>& operator>>(int& val) {FileToBuf();int x(0), f(0);while (!isdigit(*S))f |= (*S++ == '-');while (isdigit(*S))x = (x << 1) + (x << 3) + (*S++ ^ 48);val = f ? -x : x; S++;//忽略空格换行		return *this;}inline CInBuff& operator>>(long long& val) {FileToBuf();long long x(0); int f(0);while (!isdigit(*S))f |= (*S++ == '-');while (isdigit(*S))x = (x << 1) + (x << 3) + (*S++ ^ 48);val = f ? -x : x; S++;//忽略空格换行return *this;}template<class T1, class T2>inline CInBuff& operator>>(pair<T1, T2>& val) {*this >> val.first >> val.second;return *this;}template<class T1, class T2, class T3>inline CInBuff& operator>>(tuple<T1, T2, T3>& val) {*this >> get<0>(val) >> get<1>(val) >> get<2>(val);return *this;}template<class T1, class T2, class T3, class T4>inline CInBuff& operator>>(tuple<T1, T2, T3, T4>& val) {*this >> get<0>(val) >> get<1>(val) >> get<2>(val) >> get<3>(val);return *this;}template<class T = int>inline CInBuff& operator>>(vector<T>& val) {int n;*this >> n;val.resize(n);for (int i = 0; i < n; i++) {*this >> val[i];}return *this;}template<class T = int>vector<T> Read(int n) {vector<T> ret(n);for (int i = 0; i < n; i++) {*this >> ret[i];}return ret;}template<class T = int>vector<T> Read() {vector<T> ret;*this >> ret;return ret;}
private:inline void FileToBuf() {const int canRead = m_iWritePos - (S - buffer);if (canRead >= 100) { return; }if (m_bFinish) { return; }for (int i = 0; i < canRead; i++){buffer[i] = S[i];//memcpy出错			}m_iWritePos = canRead;buffer[m_iWritePos] = 0;S = buffer;int readCnt = fread(buffer + m_iWritePos, 1, N - m_iWritePos, stdin);if (readCnt <= 0) { m_bFinish = true; return; }m_iWritePos += readCnt;buffer[m_iWritePos] = 0;S = buffer;}int m_iWritePos = 0; bool m_bFinish = false;char buffer[N + 10], * S = buffer;
};template<int MOD = 1000000007>
class C1097Int
{
public:C1097Int(int iData = 0) :m_iData(iData% MOD){}C1097Int(long long llData) :m_iData(llData% MOD) {}C1097Int  operator+(const C1097Int& o)const{return C1097Int((m_iData + o.m_iData) % MOD);}C1097Int& operator+=(const C1097Int& o){m_iData = (m_iData + o.m_iData) % MOD;return *this;}C1097Int& operator-=(const C1097Int& o){m_iData = (m_iData - o.m_iData) % MOD;return *this;}C1097Int  operator-(const C1097Int& o){return C1097Int((m_iData - o.m_iData) % MOD);}C1097Int  operator*(const C1097Int& o)const{return((long long)m_iData * o.m_iData) % MOD;}C1097Int& operator*=(const C1097Int& o){m_iData = ((long long)m_iData * o.m_iData) % MOD;return *this;}C1097Int  operator/(const C1097Int& o)const{return *this * o.PowNegative1();}C1097Int& operator/=(const C1097Int& o){*this /= o.PowNegative1();return *this;}bool operator==(const C1097Int& o)const{return m_iData == o.m_iData;}bool operator<(const C1097Int& o)const{return m_iData < o.m_iData;}C1097Int pow(long long n)const{C1097Int iRet = 1, iCur = *this;while (n){if (n & 1){iRet *= iCur;}iCur *= iCur;n >>= 1;}return iRet;}C1097Int PowNegative1()const{return pow(MOD - 2);}int ToInt()const{return (m_iData + MOD) % MOD;}
private:int m_iData = 0;;
};class CCreateMat {
public://mat = { {a,c},{b,d} }是x1 = ax + by, y1 = cx + dy的转移矩阵。即{ x,y }*mat的结果是{ x1,y1 }。//{ xx,xy,y2 }到{ x1x1,x1y1,x2y2 }的转移矩阵mat={{aa,ac,cc},{2ab,ad+bc,2cd},{bb,bd,dd}}。// mat至少3*3，pre至少1*3static void X2XYY2(vector<vector<long long>>& mat, vector<vector<long long>>& pre, unsigned int x, unsigned int y, unsigned int a, unsigned int b, unsigned int c, unsigned int d, int iMod) {auto Mod = [&](long long val) { return val % iMod; };pre[0][0] = Mod(1LL * x * x); pre[0][1] = Mod(1LL * x * y); pre[0][2] = Mod(1LL * y * y);mat[0][0] = Mod(1LL * a * a); mat[0][1] = Mod(1LL * a * c); mat[0][2] = Mod(1LL * c * c);mat[1][0] = Mod(2LL * a * b); mat[1][1] = Mod(1LL * a * d + 1LL * b * c); mat[1][2] = Mod(2LL * c * d);mat[2][0] = Mod(1LL * b * b); mat[2][1] = Mod(1LL * b * d); mat[2][2] = Mod(1LL * d * d);}
};class CMatMul
{
public:CMatMul(long long llMod = 1e9 + 7) :m_llMod(llMod) {}// 矩阵乘法vector<vector<long long>> multiply(const vector<vector<long long>>& a, const vector<vector<long long>>& b) {const int r = a.size(), c = b.front().size(), iK = a.front().size();assert(iK == b.size());vector<vector<long long>> ret(r, vector<long long>(c));for (int i = 0; i < r; i++){for (int j = 0; j < c; j++){for (int k = 0; k < iK; k++){ret[i][j] = (ret[i][j] + a[i][k] * b[k][j]) % m_llMod;}}}return ret;}// 矩阵快速幂vector<vector<long long>> pow(const vector<vector<long long>>& a, vector<vector<long long>> b, long long n) {vector<vector<long long>> res = a;for (; n; n /= 2) {if (n % 2) {res = multiply(res, b);}b = multiply(b, b);}return res;}vector<vector<long long>> TotalRow(const vector<vector<long long>>& a){vector<vector<long long>> b(a.front().size(), vector<long long>(1, 1));return multiply(a, b);}
protected:const  long long m_llMod;
};class Solution {
public:int Ans(const int N, int p, int T) {const int a = (C1097Int<>(1) - C1097Int<>(p)).ToInt();const int b = C1097Int<>(p).ToInt();const int c = b;const int d = a;const int x = a, y = b;vector<vector<long long>> pre(1, vector<long long>(4)), mat(4, vector<long long>(4));CCreateMat::X2XYY2(mat, pre, x, y, a, b, c, d, 1e9 + 7);pre[0][3] = C1097Int<>(pre[0][0] + pre[0][2]).ToInt();mat[3][3] = 1;for (int i = 0; i < 3; i++){mat[i][3] = C1097Int<>(mat[i][0] + mat[i][2]).ToInt();}CMatMul matMul;C1097Int<> ans = matMul.pow(pre, mat, T)[0].back();ans *= C1097Int<>(2 * N).PowNegative1();ans *= C1097Int<>(T + 1).PowNegative1();return ans.ToInt();}
};int main() {
#ifdef _DEBUGfreopen("a.in", "r", stdin);
#endif // DEBUG	long long N;int p,t;cin >> N >> p >> t;auto res = Solution().Ans(N, p, t);cout << res;
#ifdef _DEBUG		//printf("N=%lld,", N);//Out(a, "a=");//Out(b, ",b=");/*Out(edge, "edge=");Out(que, "que=");*/
#endif // DEBUG		return 0;
}

单元测试

long long N;int p,t;TEST_METHOD(TestMethod1){auto res = Solution().Ans(2, 500000004 ,1);AssertEx(125000001, res);}TEST_METHOD(TestMethod2){auto res = Solution().Ans(2, 500000004, 0);AssertEx(C1097Int<>(8).PowNegative1().ToInt(), res);}TEST_METHOD(TestMethod3){auto res = Solution().Ans(2 ,C1097Int<>(3).PowNegative1().ToInt(), 0);auto acu = C1097Int<>(5) * C1097Int<>(36).PowNegative1();AssertEx(acu.ToInt(), res);}TEST_METHOD(TestMethod4){auto res = Solution().Ans(6, 114514, 11);AssertEx(756497239, res);}

扩展阅读

我想对大家说的话
工作中遇到的问题，可以按类别查阅鄙人的算法文章，请点击《算法与数据汇总》。
学习算法：按章节学习《喜缺全书算法册》，大量的题目和测试用例，打包下载。重视操作
有效学习：明确的目标及时的反馈拉伸区（难度合适）专注
闻缺陷则喜(喜缺)是一个美好的愿望，早发现问题，早修改问题，给老板节约钱。
子墨子言之：事无终始，无务多业。也就是我们常说的专业的人做专业的事。
如果程序是一条龙，那算法就是他的是睛
失败+反思=成功成功+反思=成功

视频课程

先学简单的课程，请移步CSDN学院，听白银讲师（也就是鄙人）的讲解。
https://edu.csdn.net/course/detail/38771
如何你想快速形成战斗了，为老板分忧，请学习C#入职培训、C++入职培训等课程
https://edu.csdn.net/lecturer/6176

测试环境

操作系统：win7 开发环境： VS2019 C++17
或者操作系统：win10 开发环境： VS2022 C++17
如无特殊说明，本算法用**C++**实现。