2017 ACM-ICPC 亚洲区（西安赛区）网络赛 B. Coin(矩阵快速幂）

传送门
题意：一个硬币正面朝上的概率为q/p，问测试k次后，正面朝上次数为偶数的概率是多少。
题解：设 xi/yi为i次测试后正面朝上为偶次的概率，则x(i+1)/y(i+1)=(1-xi/yi)* q/p+xi/yi*(1-q/p)。x1/y1=1-q/p;但这复杂度是o(k)的，超时。
再拆一下 x(i+1)=q*yi+p*xi-2*q*xi,y(i+1)=p^(i+1).yi是o(1)的，但xi还是o(k)的。再考虑一下矩阵快速幂(矩阵在代码后面）。
代码:

#include<iostream>
#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<vector>
#include<queue>
#include<stack>
#include<set>
#include<algorithm>
#include<map>
#include<time.h>
#include<math.h>
//#define pb push_back
//#define mp make_pair
#define INF 0x3f3f3f3f
using namespace std;
typedef long long int ll  ;
typedef pair<int,int>pp;
const int N=1e7+100;
const ll   mod=1e9+7;
int read()
{
    int x=0;
    char ch = getchar();
    while('0'>ch||ch>'9')ch=getchar();
    while('0'<=ch&&ch<='9')
    {
        x=(x<<3)+(x<<1)+ch-'0';
        ch=getchar();
    }
    return x;
}
/***********************************************************/
ll xx[2];
ll  p,q,k;
inline ll fast_multi(ll m, ll n, ll mod)
{
    ll ans = 0;
    while (n)
    {
        if (n & 1)
            ans += m;
        m = (m + m) % mod;
        m %= mod;
        ans %= mod;
        n >>= 1;
    }
    return ans;
}
inline ll fast_pow(ll  a, ll  n, ll  mod)//快速幂
{
    ll  ans = 1;
    while (n)
    {
        if (n & 1)
            ans = fast_multi(ans, a, mod);//不能直接乘
        a = fast_multi(a, a, mod);
        ans %= mod;
        a %= mod;
        n >>= 1;
    }
    return ans;
}
inline void extgcd(ll   a,ll   b,ll  & d,ll  & x,ll  & y)
{
    if(!b)
    {
        d=a;
        x=1;
        y=0;
    }
    else
    {
        extgcd(b,a%b,d,y,x);
        y-=x*(a/b);
    }
}
inline ll inverse(ll   a,ll   n)
{
    ll d,x,y;
    extgcd(a,n,d,x,y);
    return d==1?(x+n)%n:-1;
}
typedef struct Matrix
{
    ll mat[2][2];
}matrix;
matrix A,B;
Matrix matrix_mul(matrix a,matrix b)
{
    matrix c;
    memset(c.mat,0,sizeof(c.mat));
    ll i,j,k;
    for(ll i=0;i<2;i++)
    {
        for(ll j=0;j<2;j++)
        {
            for(ll k=0;k<2;k++)
            {
                c.mat[i][j]+=a.mat[i][k]*b.mat[k][j];
                c.mat[i][j]%=mod;
            }
        }
    }
    return c;
}
Matrix matrix_quick_power(matrix a,ll k)//矩阵快速幂0.0
{
    matrix b;
    memset(b.mat,0,sizeof(b.mat));
    b.mat[0][0]=p-q;
    b.mat[0][1]=q;
    while(k)
    {
        if(k%2==1)
        {
            b=matrix_mul(b,a);
            k-=1;
        }
        else
        {
            a=matrix_mul(a,a);
            k/=2;
        }
    }
    return b;
}

int main()
{
    int t;
    scanf("%d",&t);
    while(t--)
    {
        scanf("%lld%lld%lld",&p,&q,&k);
        A.mat[0][0]=p-2*q;
        A.mat[0][1]=0;
        A.mat[1][0]=p;
        A.mat[1][1]=p;
        B=matrix_quick_power(A,k-1);
        printf("%lld\n",fast_multi(B.mat[0][0],inverse(fast_pow(p,k,mod),mod),mod)%mod);
    }
    return 0;
}
//x(i-1)  q *  (p-2*q)  0 = xi   q
//              p      p