题解 P3216 【[HNOI2011]数学作业】
阔睡王子
2020-09-10 18:30:36
### 简单的矩阵递推练手题
前言：听完了矩阵递推的网课后我打开了团队练习题单，发现自己并不会[矩阵游戏](https://www.luogu.com.cn/problem/P1397)中矩阵满足费马小定理的证明，于是只好来切这道题，但事实上，抛开细节对代码构筑能力的考验，这道题比矩阵游戏简单，甚至可以说在你明白这是道矩阵递推的题后，简直就是水题。\
一开始先不想矩阵递推，不难推出公式：
```
long long ans=0;
for(long long i=1;i<=n;i++)
{
//ws[i]是i的位数
ans=ans*ws[i]+i;//不模的简洁版
ans=(ans*ws[i]%mod+i)%mod;//取模的细节版
}
```
明显会爆longlong， $n$ 可是小于等于 $10^{18}$ 的啊。\
思考矩阵递推，发现在矩阵里面搞位数这玩意儿实属不会，看看 $n$ 的范围，想到按位数来搞矩阵快速幂，总共做 $18$ 次，这样就不难处理位数了，而且代码跑的比香港记者还快。
然后想出矩阵如下（以下以位数为 $1$ 时为例子）：
$\begin{vmatrix}
10 & 1 & 0 \\
0 & 1 & 1 \\
0 & 0 & 1
\end{vmatrix}$
$\times$
$\begin{vmatrix}
sum[i-1] \\
i \\
1
\end{vmatrix}$
$=$
$\begin{vmatrix}
sum[i]\times10+i \\
i+1 \\
1
\end{vmatrix}$
$=$
$\begin{vmatrix}
sum[i] \\
i+1 \\
1
\end{vmatrix}$
又比如位数为3时：
$\begin{vmatrix}
1000 & 1 & 0 \\
0 & 1 & 1 \\
0 & 0 & 1
\end{vmatrix}$
$\times$
$\begin{vmatrix}
sum[i-1] \\
i \\
1
\end{vmatrix}$
$=$
$\begin{vmatrix}sum[i]\times 1000+i \\
i+1 \\
1
\end{vmatrix}$
$=$
$\begin{vmatrix}
sum[i] \\
i+1 \\
1
\end{vmatrix}$
举了两个例子后思路就清晰了，我们可以先求1到10，再求10到100，然后是100到1000等等。
```
struct node
{
int x,y;//行列数
int a[5][5];
void empty()//清空哦
{
x=y=0;
for(int i=0;i<=4;i++)
for(int j=0;j<=4;j++)
{
a[i][j]=0;
}
}
void check()//验错用的
{
for(int i=1;i<=x;i++)
{
for(int j=1;j<=y;j++)
printf("%lld ",a[i][j]);
printf("\n");
}
}
}q[20],p,p2;//q表示对应位数下使用的3X1初始矩阵,p是用来维护递推关系的柿子，p是单位矩阵0.
/*
q[1]：
0
1
1
q[2]：
123456789
10
1
*/
void pre()//初始化
{
int temp=1;
for(int i=1;i<=18;i++)
{
q[i].x=3,q[i].y=1;
q[i].a[2][1]=temp;
temp=temp*10;
q[i].a[3][1]=1;
}
p.x=p.y=p2.x=p2.y=3;
p.a[1][1]=10;
p.a[2][2]=p.a[1][2]=p.a[2][3]=p.a[3][3]=1;
p2.a[1][1]=p2.a[2][2]=p2.a[3][3]=1;
}
void work(int now,int pre,int to)//现在的位数,上一次计算到的位置，这一次该计算到的位置
{
//q[now]是递推出来的3X1矩阵，p是左边用于递推的3X3矩阵
int temp=q[now-1].a[1][1];//继承
q[now].a[1][1]=temp;
if(now==ws_ans)//如果当前位数和答案一致,可以合法的到达答案了
{
node b=quickpow(p,n-pre);//矩阵快速幂
q[now]=b*q[now];//计算出最终的矩阵
printf("%lld",q[now].a[1][1]);
return;
}
else q[now]=quickpow(p,to-pre)*q[now];//直接求解下一个位数的第一个矩阵
p.a[1][1]=p.a[1][1]*10%mod;//此处记得取模
work(now+1,to,(to+1)*10-1);//递归下一位数
}
```
然后还要注意取模的问题，不要爆longlong，基本上自己手测几个样例之后就能过了。\
送一组大样例，输入19191145141919 415411 ，输出246506。\
原代码如下：（有些变量名与上面不一样，注意区分）
```
#include<cstdio>
#define int unsigned long long
using namespace std;
int n,m,ws;
struct node
{
int x,y;
int a[5][5];
void empty()
{
x=y=0;
for(int i=0;i<=4;i++)
for(int j=0;j<=4;j++)
{
a[i][j]=0;
}
}
void check()
{
for(int i=1;i<=x;i++)
{
for(int j=1;j<=y;j++)
printf("%lld ",a[i][j]);
printf("\n");
}
}
}q[20],p,p2;
void pre()
{
int temp=1;
for(int i=1;i<=18;i++)
{
q[i].x=3,q[i].y=1;
q[i].a[2][1]=temp;
temp=temp*10;
q[i].a[3][1]=1;
}
p.x=p.y=p2.x=p2.y=3;
p.a[1][1]=10;
p.a[2][2]=p.a[1][2]=p.a[2][3]=p.a[3][3]=1;
p2.a[1][1]=p2.a[2][2]=p2.a[3][3]=1;
}
node operator*(node a,node b)
{
node o;
o.empty();
int x=a.x;
int y=a.y;
int z=b.y;
for(int i=1;i<=x;i++)
for(int j=1;j<=y;j++)
for(int k=1;k<=z;k++)
o.a[i][k]=(o.a[i][k]+(a.a[i][j]%m)*(b.a[j][k]%m))%m;
o.x=x;
o.y=z;
return o;
}
node quickpow(node a,int b)
{
node ans=p2,temp=a;
while(b)
{
if(b&1)ans=ans*temp;
temp=temp*temp;
b/=2;
}
return ans;
}
void work(int now,int pre,int to)
{
int vra=q[now-1].a[1][1];
q[now].a[1][1]=vra;
if(now==ws)
{
node b=quickpow(p,n-pre);
q[now]=b*q[now];
printf("%lld",q[now].a[1][1]);
return;
}
else q[now]=quickpow(p,to-pre)*q[now];
p.a[1][1]=p.a[1][1]*10%m;
work(now+1,to,(to+1)*10-1);
}
signed main()
{
scanf("%lld%lld",&n,&m);
pre();
int temp=n;
while(temp)
{
temp/=10;
ws++;
}
work(1,0,9);
}
```
撰写不易，求通过。