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

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

Ansi={(Ansi−1+Ansi−2)%P(Ansi−1+Ansi−2−1)%P(Ansi−1+Ansi−2)%K≠1(Ansi−1+Ansi−2)%K=1

那么，令两个3×3矩阵A,B为：

A=⎛⎝⎜110100001⎞⎠⎟,B=⎛⎝⎜ 1 0−1010001⎞⎠⎟

可以看出，只需在行向量(0 1 1)上不断右乘矩阵A和B，就可以完成递推和减1的操作。然而递推过程中我们还要对(Ansi−1+Ansi−2)modK进行查找，怎么样才能对这一过程进行加速呢？
我们把原来的斐波那契数列的第i项称为Fi。手动计算AnsimodK这一数列的前几项值，可以发现一个规律：每一次出现0，后面都紧接着两个同样的数，这个数等于这个0前面的那个数。我们把这个数称为a，那么我们发现可以将这一段写为：...,a,0,F1amodK,F2amodK,...,Fl−1amodK,0,b,...，可以发现这一段的值是由a唯一确定的，要求这一段的长度（也就是上面写的l），实际上就是求关于x的同余方程Fxa≡1(modK)的最小的解，而这个解实际上就是a关于模K的乘法逆元在数列FimodK中最小的出现位置。注意，为了使算出来的l有意义，计算1的出现位置时是不算F1,F2两项的。那么对于所有0≤a<K运行上面过程，并记录len(a)=l,next(a)=b=Fl−1amodK，若逆元不存在，就令len(a)=inf,next(a)=−1（即看成是一个无限延伸的段），这样我们就把数列AnsimodK划分成了若干段。另外，在计算数列FimodK时，只要算出一个循环节即可，因为根据某个定理，该数列的循环节长度不会超过6K，具体我不会证……
将数列AnsimodK分段后我们发现可以每次计算一个段，若这个段长度为l，就要在当前算出的矩阵上右乘Al×B，计算过程中只需要记录当前来到的点以及这一段开头的那个a即可，每计算完一段，令d=d+len(a),a=next(a)，这样我们就可以把原来的按点计算转变为按段计算，有了一定的加速。然而在面对一些情况时还是超时，注意到next(a)是唯一的，那么根据抽屉原理，数列AnsimodK一定会从某处开始循环，而且每个循环节长度不超过K段，那么我们就可以算出每计算一个循环节所需要乘的矩阵。那么整个的数列就被我们分成：不循环部分+若干个循环节+若干个段+若干个单点，每个部分中使用矩阵快速幂优化时间，这样就可以通过这道题了。
以下是本人（又臭又长的）代码：

#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;
}