HDU 4686(矩阵快速幂)

最新推荐文章于 2022-07-29 11:41:21 发布

原创最新推荐文章于 2022-07-29 11:41:21 发布 · 1.2k 阅读

1 ·

CC 4.0 BY-SA版权

------矩阵快速幂专栏收录该内容

3 篇文章

订阅专栏

本文介绍了一种利用矩阵快速幂优化递推式的方法，以解决弧梦曲线问题。通过构建特定的矩阵A并运用快速幂运算，可以在大规模数据集下高效计算出指定项的值。

An Arc of Dream is a curve defined by following function:

where
a ₀ = A0
a _i = a _i-1*AX+AY
b ₀ = B0
b _i = b _i-1*BX+BY
What is the value of AoD(N) modulo 1,000,000,007?

Input

There are multiple test cases. Process to the End of File.
Each test case contains 7 nonnegative integers as follows:
N
A0 AX AY
B0 BX BY
N is no more than 10 ¹⁸, and all the other integers are no more than 2×10⁹.

Output

For each test case, output AoD(N) modulo 1,000,000,007

转自：http://www.cnblogs.com/frog112111/archive/2013/08/21/3273660.html

分析：ai*bi = (ai-1 *ax + ay) * (bi-1 * bx + by)

=(ai-1 * bi-1 * ax * bx) + (ai-1 * ax * by) + (bi-1 * bx * ay) + (ay * by)

设p = ax*bx, q = ax*by, r = ay*bx, s = ay*by

所以ai * bi= p·(ai-1 * bi-1) + q·(ai-1) + r·(bi-1) + s

虽然可以用递推来求出每一项，但是n太大了，直接求绝对会超时的。

设 f(n) = an*bn, a(n)=an, b(n)=bn

s(n) = sum(ai * bi) , i=0,1,...n

则f(i)= p * f(i-1) + q * a(i-1) + r * b(i-1) + s

这是一个递推式，对于任何一个递推式，我们都可以用矩阵法来优化，加快速度求出第n项或前n项和。

我们可以构造一个5*5的矩阵A，使得

【f(n-1), a(n-1), b(n-1), 1, s(n-2)】* A =【f(n) ,a(n) ,b(n) , 1, s(n-1)】

=【p * f(n-1) + q * a(n-1) + r*b(n-1) + s, a(n-1) * ax + ay, b(n-1) * bx + by, 1, s(n-2) + f(n-1)】

所以我们容易得出矩阵A：【 axbx 0 0 0 1

axby ax 0 0 0

aybx 0 bx 0 0

ayay ay by 1 0

0 0 0 0 1 】

由【f(1), a(1) , b(1), 1, s(0) 】* A = 【f(2), a(2), b(2), 1, s(1)】

以此类推得，【f(1), a(1) , b(1), 1, s(0)】* A^(n-1) = 【f(n), a(n), b(n), 1, s(n-1)】

这样就可以快速的求出s(n-1)了，

其中 f(1) = a1 * b1, a(1) = a0 * ax + ay,

b(1) = b0 * bx + by, s(0) = a0 * b0

接下来就是矩阵快速幂了。

注意：n==0时，直接输出0，不然会死循环TLE的，还有就是要用long long,也要记得mod

/***********************************************
 * Author: fisty
 * Created Time: 2015-08-26 23:32:51
 * File Name   : A.cpp
 *********************************************** */
#include <iostream>
#include <cstring>
#include <deque>
#include <cmath>
#include <queue>
#include <stack>
#include <list>
#include <map>
#include <set>
#include <string>
#include <vector>
#include <cstdio>
#include <bitset>
#include <algorithm>
using namespace std;
#define Debug(x) cout << #x << " " << x <<endl
#define Memset(x, a) memset(x, a, sizeof(x))
const int INF = 0x3f3f3f3f;
typedef long long LL;
typedef pair<int, int> P;
#define FOR(i, a, b) for(int i = a;i < b; i++)
#define lson l, m, k<<1
#define rson m+1, r, k<<1|1
#define MOD 1000000007
struct Matrix{
    LL a[6][6];   
}ans, tem, A, o, _ans;
int n;
Matrix mul(Matrix x, Matrix y){
    Memset(tem.a, 0);
    for(int i = 1;i <= n; i++){
        for(int j = 1;j <= n; j++){
            for(int k = 1;k <= n; k++){
                tem.a[i][j] += (x.a[i][k] * y.a[k][j]) % MOD;
                tem.a[i][j] %= MOD;
            }
        }
    }
    return tem;
}
    
void quickpow(LL k){
    Memset(ans.a, 0);
    for(int i = 1;i <= n; i++){
        ans.a[i][i] = 1;
    }
    while(k){
        if(k & 1){
            ans = mul(ans, A);
        }
        A = mul(A, A);
        k >>= 1;
    }
}
LL m, a0, ax, ay, b0, bx, by;
int main() {
    //freopen("in.cpp", "r", stdin);
    //cin.tie(0);
    //ios::sync_with_stdio(false);
    while(~scanf("%lld%lld%lld%lld%lld%lld%lld", &m, &a0, &ax, &ay, &b0, &bx, &by)){
        if(!m) {    
            printf("0\n");
            continue;
        }
        n = 5;
        LL a1 = (a0*ax + ay) % MOD;LL b1 = (b0 * bx + by) % MOD;
        LL f1 = (a1 * b1) % MOD;
        LL s0 = (a0 * b0) % MOD;
        Memset(o.a, 0);
        o.a[1][1] = f1; o.a[1][2] = a1; o.a[1][3] = b1; o.a[1][4] = 1; o.a[1][5] = s0;
        Memset(A.a, 0);
        A.a[1][1] = (ax * bx) % MOD; A.a[1][5] = 1;
        A.a[2][1] = (ax * by) % MOD; A.a[2][2] = ax % MOD;
        A.a[3][1] = (ay * bx) % MOD; A.a[3][3] = bx % MOD;
        A.a[4][1] = (ay * by) % MOD; A.a[4][2] = ay % MOD;
        A.a[4][3] = by % MOD;        A.a[4][4] = 1;
        A.a[5][5] = 1;
        quickpow(m-1);
        _ans = mul(o, ans);
        printf("%lld\n", _ans.a[1][5]);
    }
    return 0;
}