【Author : DS】Regionals 2010, Asia - Amritapuri - E 题解

Sales Prediction

题目地址：

http://acm.bnu.edu.cn/bnuoj/problem_show.php?pid=11994

今天1Y了E题，对着俩2B的样例足足调了2.5h，不过1Y的时候相当兴奋。索性写个解题报告吧。

题目：

n , r , k 分别是 需要考察的n项 ， "猥"Fibonacci 数列的阶数 ， 考察数字的下表差。s1~sr给前r个数，a1~ar “猥”Fibonacci数列的系数。求 S_{k} + S_{k * 2} + ... + S_{k * n}

规模：

n < 1e9 ,  r , k <=8

分析：

看到Fib 数列，首先想到的就是矩阵乘法。介于矩阵乘法只能用logN的时间算出第N项，就想到用 E + A + A^{2} + ... + A^{k} 的形式累加。所以只用算出来相邻k项之间的转移矩阵，然后用log(n)的时间计算 E + A + A^{2} + ... + A^{k} 即可。比赛的时候想麻烦了。。。

比赛时解法：

1、求 r , k 的LCM  = C 。 在每一段C中划分等价类 ， 将 Const * c + j * k 化为一个等价类。并且n *= k。 也就是算1~n中 i % k == 0 的 Si 之和

2、预处理 1 ~ C 的所有值 ， 暴力O(C * C )。

3、枚举 1 ~ C 中，遇到 i % k == 0 时说明找到了一类， 如果i >= r 那么做4 ， 否则 ans = ans + a[i] , i += C 然后做4 。 最后累加即可 O(100)

4、3中，已经找到了某一个等价类的开始 i （i % k == 0） 并且 i >= C 。那么首先用Fib的方法建立1~r -> 2 - r+1 的转移矩阵 Meng ， [ O(r * r) ]  , 然后C次方求出对于一个等价类的转移矩阵（即从 x * k  转移到 x * k + C的转移矩阵） O(log(C) * r * r * r) 。 并将 k - r ~ k  放到向量 Meng2 中

5、将矩阵放入如下图的分块矩阵ZHU中，并且 n / c + 1 次方 得到 E + A + A^{2} + ... + A^{k}  矩阵 O(log(n) * 2 * 2 * 2 * r * r * r)  

6、将E + A + A^{2} + ... + A^{k} 矩阵乘以Meng2 那么得到的Meng2的（0，0）号位置就是此等价类的和 O(1 * r * r)

7、累加，输出答案。

所以最终的复杂度就是 O(log(n) * r ^ 3 * Const)

Code：

#include <iostream>
#include <stdio.h>
#include <string.h>
using namespace std;
const int MATN  = 10;
#define LL long long
const LL modo =1000000007LL;
struct Mat{		//矩阵类
    int Matn , Matm;
    LL a[MATN][MATN];
    Mat(){
        Matn = 0;
        Matm = 0;
        memset(a , 0 , sizeof(a));
    }
    void output(){
        cout << "OUTPUT" << endl;
        for (int i = 0 ; i < Matn ; ++i){
            for (int j = 0 ; j < Matm ; ++j){
                printf("%d ",a[i][j]);
            }
            printf("\n");
        }
    }
    void init(){
        Matn = 0;Matm = 0;
        memset(a , 0 , sizeof(a));
    }
    Mat operator + (const Mat & A) const{
        Mat ret = A;
        for (int i = 0 ; i < Matn ; ++i){
            for (int j = 0 ; j < Matm ; ++j)
                ret.a[i][j] = (ret.a[i][j] + a[i][j]) % modo;
        }
        return ret;

    }
    Mat operator * (const Mat & A) const{
        Mat c ;
        c.Matn = Matn;
        c.Matm = A.Matm;
        for (int i = 0 ; i < Matn ; ++i)
            for (int j = 0 ; j < A.Matm ; ++j)
                for (int k = 0 ; k < Matm ; ++k)
                    c.a[i][j] = (c.a[i][j] + (a[i][k] * A.a[k][j])%modo)%modo;
        return c;
    }
    void initI(){
        memset(a, 0 , sizeof(a));
        for (int i = 0 ; i < Matn ; ++i)
            a[i][i] = 1;
    }
    Mat power(LL k){
        Mat c = *this  , b;
        b.init();
        b.Matn = Matn ; b.Matm = Matm;
        b.initI();
        while(k){
            if (k & 1LL)
                b = b * c;
            c = c * c;
            k >>= 1LL;
        }
        return b;
    }
};
LL gcd(LL a, LL b){
   return b ? gcd(b , a % b) : a;
}
LL GCD[10][10];
void init(){			//预处理GCD
    for (int i = 1 ; i <= 8 ; ++i)
        for(int j = 1; j <= 8 ; ++j)
            GCD[i][j] = gcd(i , j);
}
const int N = 200;
LL s[N] , a[N] ;
LL n , k , r;
LL c;
Mat Meng , res , Meng2;
struct ZHU{	//用来计算 E + A + ...
    Mat a[2][2];
    void init(Mat p){
        for (int i = 0 ; i < 2 ; ++i)
            for (int j=  0 ; j < 2 ;++j){
                a[i][j].init();
                a[i][j].Matn = a[i][j].Matm = r;
            }
        a[0][0] = p;
        a[0][1].initI();
        a[1][1].initI();
    }
    void initI(){
        for (int i = 0 ; i < 2 ; ++i)
            for (int j=  0 ; j < 2 ;++j){
                a[i][j].Matn = a[i][j].Matm = r;
                a[i][i].initI();
            }
    }
    ZHU operator * (const ZHU & A) const{
        ZHU res;
        for (int i = 0 ; i < 2 ; ++i)
            for (int j= 0 ; j < 2 ;++j){
                res.a[i][j].init();
                res.a[i][j].Matn = res.a[i][j].Matm = r;
            }
        for (int i = 0 ; i < 2 ; ++i)
            for (int j = 0 ; j < 2 ; ++j)
                for (int k = 0 ; k < 2 ; ++k)
                    res.a[i][j] = (res.a[i][j] + a[i][k] * A.a[k][j]);
        return res;
    }
    Mat power(LL k){
        ZHU c = *this  , b;
        bool f = false;
        while(k){
            if (k & 1LL){
                if (!f) b = c;
                else b = b * c;
                f = true;
            }
            c = c * c;
            k >>= 1LL;
        }
        if (!f){
            b = c;
            b.a[0][1].initI();
        }
        return b.a[0][1];
    }
};
ZHU zhu;
void solve(){
    scanf("%lld%lld%lld",&n , &r , &k);
    for (int i = 1; i <= r ; ++i)
        scanf("%lld", &s[i]);
    for (int i = 1 ; i <= r ; ++i)
        scanf("%lld",&a[i]);
    c = k * r / gcd(k , r); 	//计算循环节





    for (int i = r + 1; i <= 2 * c ; ++i){	//预处理循环内S的值
        s[i] = 0;
        for (int j = 1 ; j <= r ; ++j)
            s[i] = (s[i] + (a[j] * s[i - j]) % modo) %modo;
    }
    LL ans = 0;
    n *= k;
    for (int i = 1; i <=  c && i <= n ; ++i){	//在循环节内枚举等价类的开始并且计算
        if (i % k != 0)continue;
        //printf("%d   %d  DS\n",i , (n - i) / c + 1);

        Meng.init();
        Meng.Matm = Meng.Matn = r;
        for (int j = 0 ; j < r ; ++j)
            Meng.a[j][0] = a[j + 1];
        for (int j = 0 ; j < r ; ++j)
            Meng.a[j][j + 1] = 1LL;
        //Meng.output();
        Meng= Meng.power(c);
        int ii;
        if (i <= r){		// 如果 i <= r ，那么从下一个k开始
            ans = (ans + s[i]) % modo;
            ii = c + i;
        }
        else ii = i;
        if (ii > n) continue;
        //Meng.output();
        Meng2.init();		//将前r项放到向量中
        Meng2.Matn = 1;
        Meng2.Matm = r;
        for(int j = 0 ; j < r  ; ++j)
            Meng2.a[0][j]  = s[ii - j];
        zhu.init(Meng);
        res = zhu.power( (n - ii) / (c) + 1);	// 计算 E + A...
        //res.output();
        Meng2 = Meng2 * res;	//	计算此等价类的累加和
        //cout << "DS:" << Meng2.a[0][0] << endl;
        ans = (ans + Meng2.a[0][0]) % modo;	//放入ans中。。
    }
    printf("%lld\n",ans);
}
int main(){
    int _;
    scanf("%d", &_);
    while (_--) solve();
}
/*
2
3 2 3
1 1
1 1
*/

总结：

这个题目一开始想的方法方向对了， 很自豪。不过错误的地方太多，分块矩阵乘法也没有写过，比较生疏。不过是一个好题，巩固了矩阵乘法，也练了一把Coding能力。好题啊Orz。传说现场赛只有仨人过- -爽翻了！

屌丝和孟神在长春要加油。BUAA && CUGB要给力~

DSLoveMz