ABC129F Takahashi’s Basics in Education and Learning

题目大意

有一个等差数列，有LLL个元素：s0,s1,s2,…,sL−1s_0,s_1,s_2,\dots,s_{L-1}s0​,s1​,s2​,…,sL−1​.初始元素为AAA，公差为BBB，si=A+B×is_i=A+B\times isi​=A+B×i.将这些元素从左到右拼接为一个整数，如3,7,11,15,193,7,11,15,193,7,11,15,19拼接起来得到371115193711151937111519这个整数。输出拼接后的整数模MMM的结果，数据保证等差数列中的前LLL项都不超过101810^{18}1018.

题解

当存在一对整数(l,r)(l,r)(l,r)，满足∀i∈[l,r],10d−1≤si<10d\forall i\in[l,r],10^{d-1}\leq s_i< 10^d∀i∈[l,r],10d−1≤si​<10d，设拼接到第i−1i-1i−1个整数时这个数为xxx，则拼接到第iii个整数时这个数为x+six+s_ix+si​。我们构造一个向量(x,si)(x,s_i)(x,si​)和一个矩阵CCC，使得(x,si)×C=(x×10d,si+B)(x,s_i)\times C=(x\times10^d,s_i+B)(x,si​)×C=(x×10d,si​+B)，其中si+Bs_i+Bsi​+B就是si+1s_{i+1}si+1​。

但在构造矩阵的时候，我们发现难以将sis_isi​变为si+Bs_i+Bsi​+B，所以我们将向量变为(x,si,B)(x,s_i,B)(x,si​,B)，那么：

(x,si,B)×[10d00110011]=(x×10d+si,si+B,B)(x,s_i,B)\times \left[\begin{matrix}10^d & 0 & 0 \\1 & 1 & 0 \\0 & 1 & 1\end{matrix}\right]=(x\times10^d+s_i,s_i+B,B)(x,si​,B)×⎣⎡​10d10​011​001​⎦⎤​=(x×10d+si​,si​+B,B)

那么，用矩阵快速幂就能求出对于每个ddd的答案了。

那么，怎么求对于每个ddd，矩阵要乘的次数呢？

设当前这一段每个数都有ddd位的范围为(ld,rd)(l_d,r_d)(ld​,rd​)，对于每个ddd，显然i=ldi=l_di=ld​是满足A+B×i≥10dA+B\times i\geq10^dA+B×i≥10d的最小值。变形之后就变为i≥10d−ABi\geq\dfrac{10^d-A}{B}i≥B10d−A​。那么ld=⌈10d−AB⌉l_d=\lceil\dfrac{10^d-A}{B}\rceilld​=⌈B10d−A​⌉。当然，有时候⌈10d−AB⌉\lceil\dfrac{10^d-A}{B}\rceil⌈B10d−A​⌉为负，所以得出的ldl_dld​要与000取max⁡\maxmax。显然，此时的rd=ld+1−1r_d=l_{d+1}-1rd​=ld+1​−1。

时间复杂度为O(18×27log⁡L)O(18\times27\log L)O(18×27logL)。

code

#include
using namespace std;
long long l,a,b,p,ans,mi[25],h[25];
struct matrix{long long a[4][4];matrix operator*(const matrix ax)const{matrix cx;memset(cx.a,0,sizeof(cx.a));for(int i=1;i<=3;i++){for(int j=1;j<=3;j++){for(int k=1;k<=3;k++){cx.a[i][j]=(cx.a[i][j]%p+a[i][k]%p*(ax.a[k][j]%p)%p)%p;}}}return cx;}
}re,w;
matrix ksm(matrix t,long long v){matrix c;memset(c.a,0,sizeof(c.a));if(v==0){c.a[1][1]=c.a[2][2]=c.a[3][3]=1;return c;}c=ksm(t,v/2);c=c*c;if(v%2) c=c*t;return c;
}
int main()
{scanf("%lld%lld%lld%lld",&l,&a,&b,&p);memset(re.a,0,sizeof(re.a));memset(w.a,0,sizeof(w.a));mi[0]=1;for(int i=1;i<=18;i++) mi[i]=mi[i-1]*10;for(int i=1;i<=18;i++){h[i]=(mi[i]-a)/b;if((mi[i]-a)%b>0) ++h[i];h[i]=max(h[i],0ll); }w.a[2][1]=w.a[2][2]=w.a[3][3]=1;w.a[3][2]=1;re.a[1][2]=a;re.a[1][3]=b;for(int i=0;i<=17;i++){w.a[1][1]=mi[i+1];re=re*ksm(w,min(h[i+1]-h[i],l));l-=h[i+1]-h[i];if(l<=0) break;}printf("%lld",re.a[1][1]);return 0;
}