【Bzoj3240】【Luogu1397】矩阵游戏

Source and Judge

Noi2013
Bzoj3240
Luogu1397

Problem

【Brief description】
$F[1][1]=1$
$F[i][j]=a\times F[i][j-1]+b (j!=1)$
$F[i][1]=c\times F[i-1][m]+d (i!=1)$
输出F[n][m]除以1,000,000,007的余数。
【Input】
六个整数n,m,a,b,c,d
【Output】
F[n][m]除以1,000,000,007的余数
【Limited conditions】

【Sample input】
3 4 1 3 2 6
【Sample output】
85
【Sample explanation】
样例中的矩阵为:
1 4 7 10
26 29 32 35
76 79 82 85

Record

2h
呃我的矩阵是3乘3的,时间可能是别人的两倍以上?反正问题不大

Analysis

请先思考后再展开

首先,这个n、m的大小这么奇葩,那就是有一定暗示的了
(其实有的时候我就在想为啥非要搞特大、小数据范围来提示?有检验正确性能力即可)
(UP:好像有人直接写十进制快速幂+常数优化搞过去了……)
那么我们引用伟大的费马小定理:【OI之路】11更高级数论-1定理杂烩
这个神奇玩意居然对矩阵也有效!【flag,见后文】
在mod MOD(也就是1e9+7)下,设操作矩阵为A,
那么$A^{1e9+6}=1 (\mod MOD)$ ,相当于单位矩阵
所以对于同一行下,第m个=第m%(MOD-1)个
而如果把【一行的转移+m到下一行第一个的转移】看作一个操作矩阵B,
同理第n行=第n%(MOD-1)行
综上所述,$f[n][m]=f[n\%(MOD-1)][m\%(MOD-1)]$

然后我到这里就懵逼了:这么大的数字怎么存储?还要写高精度取模???
如果你也是这么想,嗯英雄所见略同
事实上,%(MOD-1)=%( (MOD-1)*10^k ),意思就是说,
完全可以放到字符串里面,取一个模一下
真是让人涨见识的骚操作QAQ
剩下的就是推一个sb矩阵了

UP:
发现自己wa了两个点【90分,其实也该满足了】?别着急
我也是自信满满地提交,然后看别人题解才看到:a=1和c=1的特殊情况
网上清一色“要特判”,但都没讲理由?

首先,实验证明,我的两个操作矩阵$A^{MOD-1}\neq 1 (\mod MOD)$,事实上循环结长度是MOD
为什么会出现这种情况呢?【开始解决flag】
我们回顾一下,费马小定理的条件之一:gcd(A,MOD)=1
但这里A是一个矩阵呀,怎么会有gcd?咳咳,别着急。
前面我们直接啥也不管,以为网上题解说能就直接用了,所以才会有现在的状况

那么我先是问了下师兄,然后两人一起捣鼓半天,大概搞了个解释:
注意,所有公式都是在模意义下进行
我们回归到最初的公式【以行内转移举例,忽略行数,反正也就是二维】
$$
f[i]=a\times f[i-1]+b (\mod MOD)——①
$$

我们希望把它变成一个等比数列来搞出一个通项公式
【没学过高中数学没关系,我也没学】
构造一个b’使满足这个柿子:
$$
f[i]+b’=a\times (f[i-1]+b’) (\mod MOD)——②
$$

现在通过①和②推导出这道题的$b’=\frac{b}{a-1} (a\neq1)$
辣么现在就能搞出一个通项公式啦
$$
f[i]+b’=a^{i-1}\times (f[1]+b’) (\mod MOD)
$$

既然是等比数列,那就搞上费马小定理吧~
【gcd(a,MOD)=1并没有影响,因为题目条件里面限制了a的范围】
$$
f[i]+\frac{b}{a-1}=a^{ (i-1)\%(MOD-1) }\times (f[1]+\frac{b}{a-1}) (\mod MOD) (a\neq1)
$$

好了,我们搞这么多有什么用呢?
师兄的说法:用来给你十进制快速幂呀!啊那行之间怎么转移呢?不会呀
这么说原来根本就不是同一个做法好吗。
那我干嘛要写在博客上呀 仅仅因为那个$a\neq1$
没错我们现在说这么多就是为了解决我们丢了10分的问题【终于回到正题了】
不扯了
那a=1的时候,其实就变成了一个等差数列,通项公式:
$$
f[i]=f[1]+b\times (i-1) (\mod MOD)
$$

然后有个东西叫欧拉函数,$\phi(x)=在1到x的正整数中与x互质的数的个数$
显然在x为质数的时候,$\phi(x)=x-1$
那么在我们刚才矩阵乘法的时候,循环节可以记作$p=\phi(MOD)$
由此而知$f[p+1]=f[1]$
【当$a\neq1$的时候,p已经算出来了,现在算a=1的情况,p未知】
$$
f[p+1]=f[1]+b\times p (\mod MOD)
$$

那么$b\times p=0 (\mod MOD)$
而这道题中,b一不是MOD的倍数,二不是0,则$p=0 (\mod MOD)$
而p肯定又不是0,所以可以考虑把p看作MOD的倍数
哇那就是说
$$
p=MOD 【a=1】
$$

$$
p=\phi(MOD)=MOD-1 【Otherwise】
$$

哈哈搞定
~其实我自己觉得有种东扯扯西扯扯的感觉
但好歹也是有理有据的嘛
总结:
综上所述,我们在搞n和m的时候,分别根据a和c是否是1来决定MOD(详见代码)
然后经验就是,这种有特殊条件的定理,先不要引入矩阵乘法,而是把公式搞清楚,把各种细节考虑周到再去优化【虽然这道题也就10分】

Code

请先思考后再展开
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
//*******************主函数******************
#include<cstdio>
#include<cstring>
#include<cmath>
#include<queue>
#include<vector>
#include<algorithm>
#ifdef WIN32
#define BIGN "%I64d"
#else
#define BIGN "%lld"
#endif
using namespace std;
typedef long long ll;
//*******************全局常量****************
const ll MOD=1e9+7;
//*******************全局定义****************
struct martix
{
int row,col;
ll m[5][5];
martix()
{
memset(m,0,sizeof(m));
}
};
//*******************实现******************
martix cheng(martix a,martix b)
{
martix c;
int n=a.row,m=a.col,p=b.col;
c.row=n;c.col=p;
for(int i=1;i<=n;i++)
for(int j=1;j<=p;j++)
for(int k=1;k<=m;k++)
(c.m[i][j]+=a.m[i][k]*b.m[k][j])%=MOD;
return c;
}
martix pre()
{
martix c;c.row=c.col=3;
c.m[1][1]=c.m[2][2]=c.m[3][3]=1;
return c;
}
martix power(martix a,int e)
{
martix ans=pre();
while(e>0)
{
if(e&1) ans=cheng(ans,a);
a=cheng(a,a);e>>=1;
}
return ans;
}
//*******************主函数******************
char s1[1000010],s2[1000010];
int main()
{
ll mod1=MOD-1,mod2=MOD-1;
scanf("%s%s",s1+1,s2+1);
ll a,b,c,d;scanf("%lld%lld%lld%lld",&a,&b,&c,&d);
if(a==1) mod1++;if(c==1) mod2++;
ll n=0;int t1=strlen(s1+1);for(int i=1;i<=t1;i++) n=(n*10+s1[i]-'0')%mod1;
ll m=0;int t2=strlen(s2+1);for(int i=1;i<=t2;i++) m=(m*10+s2[i]-'0')%mod2;
if(n==0) n=mod1;if(m==0) m=mod2;
martix f1;f1.row=3;f1.col=1;f1.m[1][1]=1;f1.m[2][1]=b;f1.m[3][1]=d;
martix A;A.row=3;A.col=3;A.m[1][1]=a;A.m[1][2]=1;A.m[2][2]=1;A.m[3][3]=1;
martix B;B.row=3;B.col=3;B.m[1][1]=c;B.m[1][3]=1;B.m[2][2]=1;B.m[3][3]=1;
A=power(A,m-1);
if(n>1) A=cheng( A,power(cheng(B,A),n-1) );
printf("%lld",cheng(A,f1).m[1][1]);
}

本文基于 知识共享署名-相同方式共享 4.0 国际许可协议发布
本文地址:http://zory.ink/posts/e426.html
转载请注明出处,谢谢!