本文涉及知识点
【矩阵快速幂】封装类及测试用例及样例
P6601 「EZEC-2」机器
题目背景
tlx 喜欢科幻小说。
小宇宙中只剩下漂流瓶和生态球。漂流瓶隐没于黑暗里,在一千米见方的宇宙中,只有生态球里的小太阳发出一点光芒。在这个小小的生命世界中,几只清澈的水球在零重力环境中静静地飘浮着,有一条小鱼从一只水球中蹦出,跃入另一只水球,轻盈地穿游于绿藻之间。在一小块陆地上的草丛中,有一滴露珠从一片草叶上脱离,旋转着飘起,向太空中折射出一缕晶莹的阳光。
− − \qquad \qquad \qquad \qquad \qquad \qquad \qquad\qquad\qquad\qquad\qquad\qquad\qquad -- −−《三体》
在另一个宇宙,将是另一番奇景吧。
在那里,重力似乎变得微不足道了,引力机器成了司空见惯的东西。
引力机器装置内并没有重力,即若有物体在机器上运动,运动过程中只受机器给予的引力,这个力有一定几率使物体向施力物体快速移动,达到一定动力时就可以实现瞬移。
题目描述
一个引力机器由一个光滑圆轨道和 2 n 2n 2n 个小孔组成(小孔按逆时针从 1 1 1 到 2 n 2n 2n 编号,每两个相邻的小孔所夹的劣弧度数为 π n \dfrac{\pi}{n} nπ ),每个小孔与和其夹角为 π \pi π 的另一个小孔有通道相连,比如当 n = 2 n=2 n=2 时, 1 1 1 号孔和 3 3 3 号孔相连。
当 n = 2 n=2 n=2 时,这个装置的构造大概是这样的:
现在我们在 1 1 1 号孔处放一个小球,使它一直沿逆时针方向做匀速圆周运动,在不瞬移的情况下,每一秒恰好能从一个小孔运动至下一个小孔。
由于未来实验室构造奇特(内部的引力提供装置太神了!),每经过一个小孔时,有 p p p 的概率立刻瞬移(即不花费时间)到通道对面的小孔并继续沿逆时针方向做匀速圆周运动,也就是有 1 − p 1-p 1−p 的概率继续沿圆周向下一个小孔运动。
值得注意的是,每一单位时刻,小球只能瞬移一次。
简单地说,若某一时刻小球在小孔 i i i,则下一时刻它可能运动到小孔 i m o d 2 n + 1 i \bmod 2n + 1 imod2n+1 或 ( i + n ) m o d 2 n + 1 (i + n) \bmod 2n + 1 (i+n)mod2n+1,概率分别为 1 − p 1-p 1−p 和 p p p。
现在 tlx 有两个一模一样的引力机器,两个小球同时从 1 1 1 号孔开始运动。他会随机(所有可能选择的概率相同)选择一个二元组 ( i , j ) ( 1 ⩽ i ⩽ 2 n , 0 ⩽ j ⩽ t , i , j ∈ Z (i,j)( 1\leqslant i\leqslant 2n,0\leqslant j\leqslant t,i,j\in \mathbb Z (i,j)(1⩽i⩽2n,0⩽j⩽t,i,j∈Z ) 分别代表小孔编号和时间,你需要求出时间为 j j j 时两个引力机器的小孔 i i i 同时有小球停留(运动经过小孔但瞬移到对面了不算停留) 的概率。
注意:小球刚开始运动时也可能瞬移到对面的小孔。
为方便计算,我们规定:所有概率都是在模 1 0 9 + 7 10^9+7 109+7 意义下的。
输入格式
输入数据共一行,三个整数 n , p , t n,p,t n,p,t,分别代表引力装置的小孔数的一半,瞬移的概率对 1 0 9 + 7 10^9+7 109+7 取模的结果,选择的时间的范围的上界。
输出格式
共一行,一个整数,代表两个小球同时经过所选位置的概率对 1 0 9 + 7 10^9+7 109+7 取模的结果。
输入输出样例 #1
输入 #1
2 500000004 1
输出 #1
125000001
输入输出样例 #2
输入 #2
6 114514 11
输出 #2
756497239
说明/提示
【数据范围与约定】
本题采用捆绑测试。
具体计分方式如下:
- Subtask 1 1 1 ( 7 7 7 points):满足 p ∈ { 0 , 1 } p\in \{0,1\} p∈{0,1};
- Subtask 2 2 2 ( 13 13 13 points):满足 t ⩽ 20 , n ⩽ 50 t\leqslant 20,n\leqslant50 t⩽20,n⩽50;
- Subtask 3 3 3 ( 20 20 20 points):满足 t ⩽ 1 0 3 , n ⩽ 50 t\leqslant 10^3,n\leqslant50 t⩽103,n⩽50;
- Subtask 4 4 4 ( 10 10 10 points):满足 t ⩽ 1 0 3 t\leqslant 10^3 t⩽103;
- Subtask 5 5 5 ( 10 10 10 points):满足 t ⩽ 1 0 6 t\leqslant 10^6 t⩽106;
- Subtask 6 6 6 ( 15 15 15 points):满足 n ⩽ 50 n\leqslant50 n⩽50;
- Subtask 7 7 7 ( 25 25 25 points):无特殊限制。
对于 100 % 100\% 100% 的数据,满足 2 ⩽ n ⩽ 500 2\leqslant n\leqslant 500 2⩽n⩽500, 0 ⩽ p ⩽ 1 0 9 + 6 0\leqslant p\leqslant 10^9+6 0⩽p⩽109+6, 0 ⩽ t ⩽ 1 0 9 0\leqslant t \leqslant 10^9 0⩽t⩽109。
注意:不做说明的数据范围即为极限数据范围。
【样例解释 #1】
500000004 500000004 500000004 是模 1 0 9 + 7 10^9+7 109+7 意义下的 1 2 \dfrac{1}{2} 21。
下面为了方便,记 P ( i , j ) P(i,j) P(i,j) 为选择的二元组为 ( i , j ) (i,j) (i,j) 时的概率。
所有概率不为
0
0
0 的二元组有:
P
(
1
,
0
)
=
1
4
,
P
(
3
,
0
)
=
1
4
,
P
(
2
,
1
)
=
1
4
,
P
(
4
,
1
)
=
1
4
P(1,0)=\dfrac{1}{4},P(3,0)=\dfrac{1}{4},P(2,1)=\dfrac{1}{4},P(4,1)=\dfrac{1}{4}
P(1,0)=41,P(3,0)=41,P(2,1)=41,P(4,1)=41。
所有可以选择的二元组有:
(
1
,
0
)
,
(
1
,
1
)
,
(
2
,
0
)
,
(
2
,
1
)
,
(
3
,
0
)
,
(
3
,
1
)
,
(
4
,
0
)
,
(
4
,
1
)
(1,0),(1,1),(2,0),(2,1),(3,0),(3,1),(4,0),(4,1)
(1,0),(1,1),(2,0),(2,1),(3,0),(3,1),(4,0),(4,1),共
8
8
8 种。
所以总的概率:
P
=
1
8
×
1
4
×
4
+
1
8
×
0
×
4
=
1
8
P=\dfrac{1}{8}×\dfrac{1}{4}×4+\dfrac{1}{8}×0×4=\dfrac{1}{8}
P=81×41×4+81×0×4=81
在模 1 0 9 + 7 10^9+7 109+7 意义下为 125000001 125000001 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)
则结果为:
∑
(
g
0
+
g
1
)
\sum (g0+g1)
∑(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*3
static 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 _DEBUG
freopen("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++**实现。