Gauss Fibonacci hdu1588

/* 很明显要用矩阵乘法 求的是F(b) + F(b + k) + ... + F((n - 1) * k + b) 的和 A ^ (k + m) = A ^ k * A ^ m A = 0 1 1 1 K = F(k) = I * A ^ k B = F(b) = I * A ^ b B * |K I| ^ n |0 I| /*最早的做的实话是超时，用了很没有头脑的矩阵二分。终于AC了一开始我构造的矩阵是A = （1 1 0 1） 没有考虑好b = 0 的情况改成 A = (0 1 1 1)测试数据过了，但是wrong 在把int 改 long long int 就OK了 首先，我们将问题整理一下，就是对等差数列 ai=k*i+b，求所有的f(ai)之和除以M的余数 当0<=i<N 大家有没有想到，因为ai是等差数列，倘若f(ai)也是个等什么序列，那说不定就有公式求了 f(ai)显然不是等差数列，直接看上去也不是等比数列 但是如果把f(ai)换成我们刚才所说的Fibonacci矩阵呢？ 是的，可是我们对矩阵是直接求幂即可，由于矩阵加法的性质，我们要求A^ai的右上角元素之和，只要求A^ai之和的右上角元素 就矩阵这个东西来说，完全可以看作一个等比数列， 首项是：A^b 公比是：A^k 项数是：N 因为矩阵的加法对乘法也符合分配律，我们提出一个A^b来，形成这样的式子： A^b*( I + A^k + (A^k)^2 + .... + (A^k)^(N-1) ) A^b 和 A^k 显然都可以用我们之前说过的方法计算出来，这剩下一部分累加怎么解决呢 简单起见，设A^k=B 要求 G(N)=I + ... + B^(N-1)，设i=N/2 若N为偶数，G(N)=G(i)+G(i)*B^i 若N为奇数，G(N)=I+ G(i)*B + G(i) * (B^(i+1)) 我们来设置这样一个矩阵 B I O I 其中O是零矩阵，I是单位矩阵 将它乘方，得到 B^2 I+B O I 乘三方，得到 B^3 I+B+B^2 O I 乘四方，得到 B^4 I+B+B^2+B^3 O I */ */ #include<iostream> #include<string> using namespace std; int m; struct Mat { long long int b[4][4]; Mat() { for(int i = 0; i < 4; i++) for(int j = 0; j < 4; j++) b[i][j] = 0; } }; Mat util,intl; Mat mul(Mat aa, Mat bb, int n) { Mat r; int i, j, k; for(i = 0; i < n; i++) for(j = 0; j < n; j++) { r.b[i][j] = 0; for(k = 0 ;k < n; k++) r.b[i][j] += aa.b[i][k] * bb.b[k][j]; if(r.b[i][j] >= m) r.b[i][j] %= m; } return r; } Mat mal(int n, int nn) { Mat p, q; p = util; q = intl; while(n > 1) { if(n % 2 == 0) { p = mul(p, p, nn); n /= 2; } else { q = mul(q, p, nn); n--; } } return mul(q, p, nn); } int main() { int n, k, b; while( cin >> k >> b >> n >> m ) { intl.b[0][0] = 1; intl.b[1][1] = 1; intl.b[0][1] = 0; intl.b[1][0] = 0; util.b[0][0] = 0; util.b[1][1] = 1; util.b[0][1] = 1; util.b[1][0] = 1; Mat ansb, ansk; if(b == 0) { ansb = intl; } else ansb = mal(b, 2); ansk = mal(k, 2); int i, j; for(i = 0; i< 4; i++) for(j = 0; j < 4; j++) { util.b[i][j] = 0; if(i == j) intl.b[i][j] = 1; } for(i = 0; i < 2; i++) for(j = 0; j < 2; j++) util.b[i][j] = ansk.b[i][j]; util.b[0][2] = 1; util.b[1][3] = 1; util.b[2][2] = 1; util.b[3][3] = 1; Mat anssum; anssum = mal(n, 4); Mat ans; ans.b[0][0] = anssum.b[0][2] % m; ans.b[0][1] = anssum.b[0][3] % m; ans.b[1][0] = anssum.b[1][2] % m; ans.b[1][1] = anssum.b[1][3] % m; Mat sum; sum = mul(ansb, ans, 4); printf("%I64d/n", sum.b[0][1] % m); } } code:

/*最早的做的实话是超时，用了很没有头脑的矩阵二分。终于AC了一开始我构造的矩阵是A = （1 1 0 1） 没有考虑好b = 0 的情况改成 A = (0 1 1 1)测试数据过了，但是wrong
在把int 改 long long int 就OK了
首先，我们将问题整理一下，就是对等差数列 ai=k*i+b，求所有的f(ai)之和除以M的余数

当0<=i<N

大家有没有想到，因为ai是等差数列，倘若f(ai)也是个等什么序列，那说不定就有公式求了

f(ai)显然不是等差数列，直接看上去也不是等比数列

但是如果把f(ai)换成我们刚才所说的Fibonacci矩阵呢？

是的，可是我们对矩阵是直接求幂即可，由于矩阵加法的性质，我们要求A^ai的右上角元素之和，只要求A^ai之和的右上角元素

就矩阵这个东西来说，完全可以看作一个等比数列，
首项是：A^b
公比是：A^k
项数是：N

因为矩阵的加法对乘法也符合分配律，我们提出一个A^b来，形成这样的式子：
A^b*( I + A^k + (A^k)^2 + .... + (A^k)^(N-1) )

A^b 和 A^k 显然都可以用我们之前说过的方法计算出来，这剩下一部分累加怎么解决呢

简单起见，设A^k=B
要求 G(N)=I + ... + B^(N-1)，设i=N/2
若N为偶数，G(N)=G(i)+G(i)*B^i
若N为奇数，G(N)=I+ G(i)*B + G(i) * (B^(i+1))

 
我们来设置这样一个矩阵
B I
O I
其中O是零矩阵，I是单位矩阵

将它乘方，得到
B^2 I+B
O   I
乘三方，得到
B^3 I+B+B^2
O   I
乘四方，得到
B^4 I+B+B^2+B^3
O   I
*/
#include<iostream>//2311179 2010-04-08 11:55:54 Accepted 1588 15MS 296K 2460 B C++ 悔惜晟
#include<string>
using namespace std;
int m;

struct Mat
{
    long long int a[2][2];//用int 会超出范围
    long long int b[4][4];
}util,intl;

Mat mul2(Mat aa, Mat bb)
{
    Mat r;
    int i, j, k;
    for(i = 0; i < 2; i++)
      for(j = 0; j < 2; j++)
      {   
          r.a[i][j] = 0;
        for(k = 0 ;k <2; k++)
            r.a[i][j] += aa.a[i][k] * bb.a[k][j];
        if(r.a[i][j] >= m)
            r.a[i][j] %= m;
    }
    return r;
}

Mat mul4(Mat aa, Mat bb)
{
    Mat r;
    int i, j, k;
    for(i = 0; i < 4; i++)
      for(j = 0; j < 4; j++)
      {   
          r.b[i][j] = 0;
        for(k = 0 ;k < 4; k++)
            r.b[i][j] += aa.b[i][k] * bb.b[k][j];
        if(r.b[i][j] >= m)
            r.b[i][j] %= m;
    }
    return r;
}

Mat mal2(int n)
{
    Mat p, q;
    p = util;
    q = intl;
    while(n > 1)
    {
        if(n % 2 ==0)
        {
            p = mul2(p, p);
            n /= 2;
        }
        else
        {
            q = mul2(q, p);
            n--;
        }
    }
    return mul2(q, p);
}

Mat mal4(int n)
{
    Mat p, q;
    p = util;
    q = intl;
    while(n > 1)
    {
        if(n % 2 == 0)
        {
            p = mul4(p, p);
            n /= 2;
        }
        else
        {
            q = mul4(q, p);
            n--;
        }
    }
    return mul4(q, p);
}

int main()
{
    int n, k, b;
    while(cin>>k>>b>>n>>m)
    {
  //memset(util, 0 ,sizeof(util));
        intl.a[0][0] = 1;
        intl.a[1][1] = 1;
        intl.a[0][1] = 0;
        intl.a[1][0] = 0;

        util.a[0][0] = 0;
        util.a[1][1] = 1;
        util.a[0][1] = 1;
        util.a[1][0] = 1;
        Mat ansb, ansk;
 //ansb = mal2(b);
 if(b == 0)
 {
  ansb = intl;
 }
 else
  ansb = mal2( b );
 ansk = mal2(k);
 int i, j;
 for(i = 0;  i< 4; i++)
  for(j =0; j < 4; j++)
  {
   util.b[i][j] = 0;
   if(i == j)
   intl.b[i][j] = 1;
  }

 for(i = 0; i < 2; i++)
  for(j = 0; j < 2; j++)
   util.b[i][j] = ansk.a[i][j];

 util.b[0][2] = 1;
 util.b[1][3] = 1;
 util.b[2][2] = 1;
 util.b[3][3] = 1;

 Mat anssum;
 anssum = mal4(n);
 Mat ans;

 ans.a[0][0] = anssum.b[0][2] % m;
 ans.a[0][1] = anssum.b[0][3] % m;
 ans.a[1][0] = anssum.b[1][2] % m;
 ans.a[1][1] = anssum.b[1][3] % m;

 Mat sum;
 sum = mul2(ansb, ans);

 printf("%I64d/n", sum.a[0][1] % m);
    }
}