普通的快速幂

问题：[COGS 1433]圣庙里的汉诺塔

计算$a^b\mod m$，其中$0\leq a,b\leq10^{16}$。

我们可以写出如下的暴力求幂程序，复杂度$O(n)$:

typedef long long ll;
ll pow_mod(ll a,ll b,ll m)
{
  	ll ans=1;
  	while(b--)
      	ans=(ans*a)%m;
  	return ans;
}

不过注意到对于任意正整数k，我们都可以拆成至多$log k$个二进制下的1,所以我们可以计算出$a^{2^0},a^{2^1},a^{2^2},\cdots$（可以递推出来）,然后再用$\log b$次乘法将其合并起来，复杂度为$O(\log n)$，代码如下：

ll fast_pow_mod(ll a,ll b,ll m)
{
    ll ans=1;
  	while(b)
    {
      	if(b&1)
          	ans*=a;
      	a*=a;
      	b>>=1;
    }
  	return ans;
}

由快速幂魔改的快速乘

问题：[COGS 2037]Asm.def大点兵

由于$a\times b$可能溢出，因此不能直接乘起来再取膜，考虑到一个事实： $$ a^b= \begin{cases}

1 & b=0 \\
\underbrace { a \times \cdots \times a }_b & b \ge 1

\end{cases} $$ $$ a\times b= \begin{cases} 0&b=0
\underbrace{a+\cdots+a}_b&b\not=0 \end{cases} $$ 所以我们把上面的代码稍微改一改就好了：

ll fast_mul_mod(ll a,ll b,ll m)
{
    ll ans=0;
  	while(b)
    {
      	if(b&1)
          	ans+=a;
      	a+=a;
      	b>>=1;
    }
  	return ans;
}

矩阵快速幂

问题：[COGS 1493]递推关系

先用蓝书上的方法搞出来友矩阵，然后利用快速幂的思想瞎搞一波。代码如下：

#include<bits/stdc++.h>
using namespace std;
const int maxn=20;
typedef long long ll;
int n,m,dd;
ll a[maxn],f[maxn];
struct matrix
{
    int x,y;
    ll d[maxn][maxn];
    matrix()
    {
        x=y=0;
        memset(d,0,sizeof(d));
    }
    matrix(int s,int mx,int my)
    {
        x=mx,y=my;
        for(int i=0;i<x;++i)
            for(int j=0;j<y;++j)
                d[i][j]=0;
        if(s==1)
            for(int i=0;i<x;++i)
                d[i][i]=1;
        if(s==-1)
        {
            for(int i=1;i<x;++i)
                d[i-1][i]=1;
            for(int i=y-1,j=0;i>=0;--i,++j)
                d[x-1][i]=a[j];
        }
    }
    const matrix operator*(const matrix& n) const
    {
        matrix ans(0,x,n.y);
        for(int i=0;i<ans.x;++i)
            for(int k=0;k<ans.y;++k)
                for(int j=0;j<y;++j)
                    ans.d[i][k]=(ans.d[i][k]+d[i][j]*n.d[j][k])%m;
        return ans;
    }
};
const matrix pow_mod(matrix & a,int t)
{
    matrix ans(1,dd,dd);
    while(t)
    {
        if(t&1)
            ans=ans*a;
        a=a*a;t>>=1;
    }
    return ans;
}
void init()
{
    freopen("recurrences.in","r",stdin);
    freopen("recurrences.out","w",stdout);
}
void print_matrix(matrix& a)
{
    cout<<"---------------\n";
    for(int i=0;i<a.x;++i)
    {
        for(int j=0;j<a.y;++j)
        {
            cout<<a.d[i][j]<<' ';
        }
        cout<<endl;
    }
    cout<<"---------------\n\n";
}
int main()
{
    init();
    while(cin>>dd>>n>>m)
    {
        if(!(dd||m||n))
            return 0;
        for(int i=0;i<dd;++i)
            cin>>a[i];
        for(int i=0;i<dd;++i)
            cin>>f[i];
        matrix a(-1,dd,dd);
        //DEBUG_print_matrix(a);
        matrix mt(0,dd,1);
        for(int i=0;i<dd;++i)
            mt.d[i][0]=f[i];
        mt=pow_mod(a,n-dd)*mt;
        //print_matrix(mt);
        //for(int i=0;i<mt.x;++i)
        //    cout<<(mt.d[i][0])%m<<' ';
        //cout<<endl;
        cout<<mt.d[mt.x-1][0]%m<<endl;
    }
}