SPOJ AMR 10E Stocks Prediction 矩阵快速幂 + 二分求和

最新推荐文章于 2019-04-08 19:43:29 发布

Gatevin

最新推荐文章于 2019-04-08 19:43:29 发布

阅读量984

点赞数

分类专栏： ACM_Dichotomy SPOJ ACM_Matrix_quick_pow 文章标签： SPOJ AMR 10E Stocks Prediction 矩阵快速幂

本文链接：https://blog.csdn.net/gatevin/article/details/38305145

版权

SPOJ 同时被 3 个专栏收录

17 篇文章 0 订阅

订阅专栏

ACM_Matrix_quick_pow

12 篇文章 0 订阅

订阅专栏

ACM_Dichotomy

4 篇文章 0 订阅

订阅专栏

题目大意：

每个月的销量满足递推式：

S(n) = a1*S(n - 1) + a2*S(n - 2) + a3*S(n - 3) + a4*S(n - 4) + ... + aR*S(n - R)

在知道 T, a1~aR， R, N, K 和 S1~SR的情况下，要求求出 sigma(S( i * K)) MOD 1000000007 ( i 从 1 到 N)

数据范围：

T <= 40

1 <= N <= 1000000000

1 <= R <= 8

1 <= K <= 8

0 <= All other input values < 1000000007

大致思路：

首先从数据范围和这个递推式来看是矩阵快速幂的题

不难发现构造这样的一个R*R的矩阵满足关系：

我们用A表示上面的递推式中的R*R的那个矩阵，那么对于前面那个向量，每次乘上A^k之后都会变成(S(n + k)...)

那么对于初始的向量( S(R) S(R - 1) ... S(1) ) 如果这个向量当中包括 S(k) 我们可以直接对于每次要算的 S( i * k) 求和

也就是说这个向量乘上( I + A^k + (A^k)^2 + (A^k)^3 + ... + (A^k)^(N - 1))之后对应的 S(k) 所在的那个位置就变成了要求的和

而对于那个矩阵型的等比数列求和可以直接用二分求和(常用的技巧)，这样就可以在限制的时间内完成计算了

当然还可以在算出A^k 之后将得到的军阵用来构造一个新的(R +1)*(R + 1) 的矩阵，这个学长的做法,,,得到的新的转换矩阵每乘上一次就对应后k项，同时也完成了求和操作，，

这个想法的确挺巧妙地吧，复杂度比用二分法求等比数列的和也要低一些

两种方法的代码都贴一下吧，，

方法一，等比数列求和的方法：

Result : Accepted Memory : 3379 KB Time : 2960 ms

/*
 * Author: Gatevin
 * Created Time:  2014/7/30 13:20:23
 * File Name: 123.cpp
 */
#include<iostream>
#include<sstream>
#include<fstream>
#include<vector>
#include<list>
#include<deque>
#include<queue>
#include<stack>
#include<map>
#include<set>
#include<bitset>
#include<algorithm>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cctype>
#include<cmath>
#include<ctime>
#include<iomanip>
using namespace std;
const double eps(1e-8);
typedef long long lint;

int t,r,k,siz;
lint n;
lint s[10];
lint h[10];
const lint mod = 1000000007LL;
struct Matrix
{
    lint a[11][11];
    Matrix()
    {
        memset(a, 0, sizeof(a));
        for(int i = 1; i <= 10; i++)
        {
            a[i][i] = 1;
        }
    }
};

Matrix operator * (const Matrix & m1, const Matrix & m2)
{
    Matrix m;
    for(int i = 1; i <= siz; i++)
    {
        for(int j = 1; j <= siz; j++)
        {
            m.a[i][j] = 0;
            for(int k = 1; k <= siz; k++)
            {
                m.a[i][j] = (m.a[i][j] + (m1.a[i][k]* m2.a[k][j]) % mod) % mod;
            }
        }
    }
    return m;
}

Matrix operator + (const Matrix & m1, const Matrix & m2)
{
    Matrix m;
    memset(m.a, 0, sizeof(m.a));
    for(int i = 1; i <= siz; i++)
    {
        for(int j = 1; j <= siz; j++)
        {
            m.a[i][j] = (m1.a[i][j] + m2.a[i][j]) % mod;
        }
    }
    return m;
}

Matrix quick_pow(Matrix base, lint pow)
{
    Matrix I;
    while(pow)
    {
        if(pow & 1)
        {
            I = I * base;
        }
        base = base * base;
        pow >>= 1;
    }
    return I;
}

Matrix Geo(Matrix base, lint pow)
{
    Matrix I;
    if(pow == 1) return I;
    if(pow & 1) return Geo(base, pow - 1) + quick_pow(base, pow - 1);
    else return Geo(base, pow >> 1)*(I + quick_pow(base, pow >> 1));
}

int main()
{
    scanf("%d",&t);
    while(t--)
    {
        scanf("%I64d%d%d",&n,&r,&k);
        siz = r;
        for(int i = 1; i <= r; i++)
        {
            scanf("%I64d",&s[i]);
        }
        for(int i = 1; i <= r; i++)
        {
            scanf("%I64d",&h[i]);
        }
        Matrix A;
        memset(A.a, 0, sizeof(A.a));
        for(int i = 1; i <= r; i++)
        {
            A.a[i][1] = h[i];
            if(i != r)
            A.a[i][i + 1] = 1;
        }
        Matrix last = Geo(quick_pow(A, k), n);
        lint answer = 0;
        if(k <= r)
        {
            for(int i = r; i >= 1; i--)
            {
                answer = (answer + (s[i]*last.a[r - i + 1][r - k + 1]) % mod) % mod;
            }
            printf("%I64d\n",answer);
        }
        else
        {
            for(int i = r + 1; i <= k; i++)
            {
                s[i] = 0;
                for(int j = 1; j <= r; j++)
                {
                    s[i] = (s[i] + (s[i - j]*h[j]) % mod) % mod;
                }
            }
            for(int i = 0; i <= r - 1; i++)
            {
                answer = (answer + (s[k - i]*last.a[i + 1][1]) % mod) % mod;
            }
            printf("%I64d\n",answer);
        }
    }
    return 0;
}

方法二，利用得到的A^k构造新的矩阵的方法：

/*
 * Author:  Eyelids
 * Created Time:  2014/7/30 16:41:47
 * File Name: E.cpp
 */
#include<iostream>
#include<sstream>
#include<fstream>
#include<vector>
#include<list>
#include<deque>
#include<queue>
#include<stack>
#include<map>
#include<set>
#include<bitset>
#include<algorithm>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cctype>
#include<cmath>
#include<ctime>
using namespace std;
const double eps(1e-8);
typedef long long lint;
#define clr(x) memset( x , 0 , sizeof(x) )
#define sz(v) ((int)(v).size())
#define rep(i, n) for (int i = 0; i < (n); ++i)
#define repf(i, a, b) for (int i = (a); i <= (b); ++i)
#define repd(i, a, b) for (int i = (a); i >= (b); --i)
#define clrs( x , y ) memset( x , y , sizeof(x) )
typedef long long LL;
const LL mod = 1000000007LL;
int n, k, r;
LL s[110], a[110], A[10][10], B[10][10];

void Mul( LL A[][10], LL B[][10] ) {
    LL ret[10][10];
    for ( int i = 1; i < 10; i ++ )
        for ( int j = 1; j < 10; j ++ ) {
            ret[i][j] = 0;
            for ( int k = 1; k < 10; k ++ )
                ret[i][j] = ( ret[i][j] + A[i][k] * B[k][j] ) % mod;
        }
    memcpy( A, ret, sizeof(ret) );
}
void Pow( LL A[][10], int k ) {
    LL ret[10][10];
    memset( ret, 0, sizeof(ret) );
    for ( int j = 1; j < 10; j ++ ) ret[j][j] = 1;
    while ( k ) {
        if ( k & 1 ) Mul( ret, A );
        Mul( A, A );
        k >>= 1;
    }
    memcpy( A, ret, sizeof(ret) );
}

int main() {
    int T;
    cin >>T;
    while ( T -- ) {
        cin >>n>>r>>k;
        memset( s, 0, sizeof(s) );
        memset( B, 0, sizeof(B) );
        memset( a, 0, sizeof(a) );
        for ( int i = 1; i <= r; i ++ ) cin >>s[i];
        for ( int i = 1; i <= r; i ++ ) cin >>a[i];
        
        memset( A, 0, sizeof(A) );
        for ( int i = r + 1; i <= 100; i ++ )
            for ( int j = 1; j <= r; j ++ ) {
                s[i] += s[i - j] * a[j];
                s[i] %= mod;
            }

        for ( int i = 1, j = k + r - 1; i <= r; i ++, j -- ) 
            A[1][i] = s[j];
        A[1][r + 1] = 0;
        /*for ( int i = 1; i <= r + 1; i ++ )
            cout <<A[1][i]<<" ";
        cout <<endl;*/
            
        for ( int i = 1; i <= r; i ++ )     
            B[i][1] = a[i];
        for ( int i = 2; i <= r; i ++ ) 
            B[i - 1][i] = 1;
       
        /*for ( int i = 1; i <= r; i ++ ) {
            for ( int j = 1; j <= r; j ++ )
                cout <<B[i][j]<<" ";
            cout <<endl;
        }  */
        Pow( B, k );
        
        B[r][r + 1] = B[r + 1][r + 1] = 1;
        Pow( B, n );
        /*for ( int i = 1; i <= r + 1; i ++ ) {
            for ( int j = 1; j <= r + 1; j ++ )
                cout <<B[i][j]<<" ";
            cout <<endl;
        } */

        Mul( A, B );
        
        cout <<A[1][r + 1]<<endl;
    }
    
    return 0;
}