hdu4686Arc of Dream 矩阵快速幂

//a[n] = ax*a[n-1] + ay//b[n] = bx*b[n-1] + by//给出n求segma(a[i]*b[i])(0<=i<=n-1)//f[n] = a[n]*b[n] = ax*bx*f[n-1] + ax*by*a[n-1] + ay*bx*b[n-1] + ay*by ;//s[n] = segma(a[i]*b[i]) = s[n-1] + f[n]//可以构造矩阵// |f[n]  | |ax*bx ax*by ay*by ay*bx  0| |f[n-1]|// |a[n]  | |  1     0     0     0    0| |a[n-1]|// |b[n]  | |  0    ax     0     ay   0| |b[n-1]|// |1     | |  0     0     0     1    0| | 1    |// |s[n]  | |  1     0     0     0    1| |s[n-1]|//然后矩阵快速幂就可求解#include<cstdio>#include<cstring>#include<iostream>using namespace std ;typedef long long ll ;const ll mod = 1e9+7 ;const int maxn = 10 ;struct node{    ll p[maxn][maxn] ;};void debug(node a){    for(int i = 1;i <= 5;i++)      for(int j = 1;j <= 5;j++)      printf("%lld%c" , a.p[i][j] , j ==  5 ?'\n':' ') ;    puts("") ;}node mul(node a , node b){    node c ;    memset(c.p , 0 , sizeof(c.p)) ;    for(int i = 1;i <= 5;i++)      for(int j = 1;j <= 5;j++)        for(int k = 1;k <= 5;k++)        c.p[i][j] = (c.p[i][j] + a.p[i][k]*b.p[k][j])%mod ;    return c ;}node pow(node a , ll k){    node c ;    memset(c.p , 0 , sizeof(c.p)) ;    for(int i = 1;i <= 5;i++)    c.p[i][i] = 1 ;    while(k)    {        if(k&1)          c = mul(c , a) ;        a = mul(a , a) ;        k >>= 1 ;    }    return c ;}int main(){    //freopen("in.txt" , "r" , stdin) ;    ll ax , ay , a0 ;    ll bx , by , b0 ;    ll n ;    while(~scanf("%lld" , &n))    {        scanf("%lld%lld%lld" , &a0 , &ax , &ay ) ;        scanf("%lld%lld%lld" , &b0 , &bx , &by ) ;        node a ;        memset(a.p , 0 , sizeof(a.p)) ;        a.p[1][1] = (ax*bx)%mod ;        a.p[1][2] = (ax*by)%mod ;        a.p[1][3] = (ay*bx)%mod ;        a.p[1][4] = (ay*by)%mod ;        a.p[2][2] = ax%mod ;        a.p[2][4] = ay%mod ;        a.p[3][3] = bx%mod ;        a.p[3][4] = by%mod ;        a.p[4][4] = 1 ;        a.p[5][1] = 1 ;        a.p[5][5] = 1 ;        if(n == 0)        {            puts("0") ;            continue ;        }        node c = pow(a , n-1) ;        ll a1 =  (a0*ax%mod + ay)%mod ;        ll b1 =  (b0*bx%mod + by)%mod ;        ll s1 =  a0*b0%mod ;        ll f1 =  a1*b1%mod ;        node x ;        memset(x.p , 0 , sizeof(x.p)) ;        x.p[1][1] = f1 ;        x.p[2][1] = a1 ;        x.p[3][1] = b1 ;        x.p[4][1] = 1 ;        x.p[5][1] = s1 ;        node ans = mul(c , x) ;        printf("%lld\n" , ans.p[5][1]%mod) ;    }}