【NOI2011T1】兔农-矩阵快速幂+乘法逆元

测试地址：兔农
做法：这题太神了orz，在跪拜了大佬的题解，再自己调了很久之后终于A掉了……这真的是第一题吗？！
这题需要用到矩阵快速幂以及乘法逆元。
把要求的数列的第$i$项称为$Ans_i$，考虑计算$Ans_i \mod P$。可以得到递推式：

$$Ans_i=\left\{ \begin{array}{} (Ans_{i-1}+Ans_{i-2})\%P & (Ans_{i-1}+Ans_{i-2})\%K\ne 1\\ (Ans_{i-1}+Ans_{i-2}-1)\%P & (Ans_{i-1}+Ans_{i-2})\%K=1\\ \end{array} \right.$$
那么，令两个 $3\times3$矩阵 $A,B$为：
$$A=\left( \begin{array}{} 1 & 1 & 0\\ 1 & 0 & 0\\ 0 & 0 & 1 \end{array} \right), B=\left( \begin{array}{} \ 1 & 0 & 0\\ \ 0 & 1 & 0\\ -1 & 0 & 1 \end{array} \right)$$
可以看出，只需在行向量 $(0 \ 1 \ 1)$上不断右乘矩阵 $A$和$B$，就可以完成递推和减1的操作。然而递推过程中我们还要对 $(Ans_{i-1}+Ans_{i-2})\mod K$进行查找，怎么样才能对这一过程进行加速呢？
我们把原来的斐波那契数列的第 $i$项称为$F_i$。手动计算 $Ans_i\mod K$这一数列的前几项值，可以发现一个规律：每一次出现 $0$，后面都紧接着两个同样的数，这个数等于这个$0$前面的那个数。我们把这个数称为 $a$，那么我们发现可以将这一段写为：$...,a,0,F_1a\mod K,F_2a\mod K,...,F_{l-1}a\mod K,0,b,...$，可以发现这一段的值是由 $a$唯一确定的，要求这一段的长度（也就是上面写的$l$），实际上就是求关于 $x$的同余方程$F_xa\equiv 1(\mod K)$的最小的解，而这个解实际上就是 $a$关于模$K$的乘法逆元在数列 $F_i\mod K$中最小的出现位置。注意，为了使算出来的 $l$有意义，计算$1$的出现位置时是不算 $F_1,F_2$两项的。那么对于所有 $0\le a<K$运行上面过程，并记录 $len(a)=l,next(a)=b=F_{l-1}a\mod K$，若逆元不存在，就令 $len(a)=inf,next(a)=-1$（即看成是一个无限延伸的段），这样我们就把数列 $Ans_i\mod K$划分成了若干段。另外，在计算数列 $F_i\mod K$时，只要算出一个循环节即可，因为根据某个定理，该数列的循环节长度不会超过 $6K$，具体我不会证……
将数列 $Ans_i\mod K$分段后我们发现可以每次计算一个段，若这个段长度为 $l$，就要在当前算出的矩阵上右乘$A^l\times B$，计算过程中只需要记录当前来到的点以及这一段开头的那个 $a$即可，每计算完一段，令$d=d+len(a),a=next(a)$，这样我们就可以把原来的按点计算转变为按段计算，有了一定的加速。然而在面对一些情况时还是超时，注意到 $next(a)$是唯一的，那么根据抽屉原理，数列 $Ans_i\mod K$一定会从某处开始循环，而且每个循环节长度不超过 K <script type="math/tex" id="MathJax-Element-633">K</script>段，那么我们就可以算出每计算一个循环节所需要乘的矩阵。那么整个的数列就被我们分成：不循环部分+若干个循环节+若干个段+若干个单点，每个部分中使用矩阵快速幂优化时间，这样就可以通过这道题了。
以下是本人（又臭又长的）代码：

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <iostream>
#include <algorithm>
#define ll long long
using namespace std;
ll N,K,P,fib[6000010],ap[1000010]={0},next[1000010];
int len[1000010],vis[1000010]={0};
struct matrix {ll s[3][3];} now,M[70],A,B,Fst,Pr;

void init()
{
  scanf("%lld%lld%lld",&N,&K,&P);
  now.s[0][0]=1,now.s[0][1]=0,now.s[0][2]=0;
  now.s[1][0]=0,now.s[1][1]=1,now.s[1][2]=0;
  now.s[2][0]=0,now.s[2][1]=0,now.s[2][2]=1;
  Fst=now,Pr=now;
  A.s[0][0]=1,A.s[0][1]=1,A.s[0][2]=0;
  A.s[1][0]=1,A.s[1][1]=0,A.s[1][2]=0;
  A.s[2][0]=0,A.s[2][1]=0,A.s[2][2]=1;
  B.s[0][0]=1,B.s[0][1]=0,B.s[0][2]=0;
  B.s[1][0]=0,B.s[1][1]=1,B.s[1][2]=0;
  B.s[2][0]=-1,B.s[2][1]=0,B.s[2][2]=1;
}

ll exgcd(ll a,ll b,ll &x,ll &y)
{
  ll x0=1,y0=0,x1=0,y1=1;
  while(b)
  {
    ll tmp,q;
    q=a/b;
    tmp=x0,x0=x1,x1=tmp-q*x1;
    tmp=y0,y0=y1,y1=tmp-q*y1;
    tmp=a,a=b,b=tmp%b;
  }
  x=x0,y=y0;
  return a;
}

ll find_inverse(ll x)
{
  if (x==0) return 0;
  ll px,py,d=exgcd(x,K,px,py);
  if (d>1) return 0;
  else return (px%(K/d)+(K/d))%(K/d);
}

matrix mult(matrix a,matrix b)
{
  matrix S;
  memset(S.s,0,sizeof(S.s));
  for(int i=0;i<=2;i++)
    for(int j=0;j<=2;j++)
      for(int k=0;k<=2;k++)
        S.s[i][j]=(S.s[i][j]+a.s[i][k]*b.s[k][j])%P;
  return S;
}

matrix power(ll x)
{
  matrix S;
  memset(S.s,0,sizeof(S.s));
  for(int i=0;i<=2;i++) S.s[i][i]=1;
  ll p=0;
  while(x)
  {
    if (x&1) S=mult(S,M[p]);
    x>>=1;p++;
  }
  return S;
}

int main()
{
  init();

  ll p=2;fib[0]=0,fib[1]=fib[2]=1;
  do
  {
    fib[++p]=(fib[p-1]+fib[p-2])%K;
    ap[fib[p]]=ap[fib[p]]?ap[fib[p]]:p;
  }
  while(fib[p-1]!=0||fib[p]!=1);

  for(int i=0;i<K;i++)
  {
    ll inv=find_inverse(i);
    if (inv&&ap[inv])
    {
      len[i]=ap[inv];
      next[i]=(fib[len[i]-1]*i)%K;
    }
    else len[i]=-1,next[i]=-1;
  }

  ll d=0,a=1;bool flag=0;
  M[0]=A;vis[1]=d;
  for(int i=1;i<=69;i++) M[i]=mult(M[i-1],M[i-1]);
  while(a!=-1&&len[a]!=-1&&d+len[a]<=N)
  {
    now=mult(now,power(len[a]));
    now=mult(now,B);
    d+=len[a];a=next[a];
    if (vis[a]) {flag=1;break;}
    else vis[a]=d;
  }

  if (flag)
  {
    ll pr=d-vis[a],fst,aa=a;
    d=0;a=1;
    while(d!=vis[aa])
    {
      Fst=mult(Fst,power(len[a]));
      Fst=mult(Fst,B);
      d+=len[a];a=next[a];
    }
    fst=d;
    d=vis[aa];a=aa;
    while(d!=pr+vis[aa])
    {
      Pr=mult(Pr,power(len[a]));
      Pr=mult(Pr,B);
      d+=len[a];a=next[a];
    }
    now=Fst;
    M[0]=Pr;
    for(int i=1;i<=69;i++) M[i]=mult(M[i-1],M[i-1]);
    now=mult(now,power((N-vis[aa])/pr));
    M[0]=A;
    for(int i=1;i<=69;i++) M[i]=mult(M[i-1],M[i-1]);
    d=fst+pr*((N-vis[aa])/pr);a=aa;
    while(a!=-1&&len[a]!=-1&&d+len[a]<=N)
    {
      now=mult(now,power(len[a]));
      now=mult(now,B);
      d+=len[a];a=next[a];
    }
    now=mult(now,power(N-d));
  }
  else now=mult(now,power(N-d));

  printf("%lld",((now.s[1][0]+now.s[2][0])%P+P)%P);

  return 0;
}