根据题意:最后一步是寻求f(b) + f(k + b) + f(2 * k + b) + …+ f((n-1) * k + b)

清除f(b) = A^b

间A =

1 1

1 0

所以sum(n - 1) = A^b(E + A^ k + A ^(2 * k) + … + A ^((n - 1) * k)

设D = A^k

sum(n-1) = A^b(E + D + D ^ 2 + … + D ^(n - 1))

括号中的部分就能够二分递归求出来了

而单个矩阵就能够用矩阵高速幂求出来

/*************************************************************************

    > File Name: hdu1588.cpp

    > Author: ALex

    > Mail: zchao1995@gmail.com

    > Created Time: 2015年03月12日 星期四 18时25分07秒

 ************************************************************************/

#include <map>

#include <set>

#include <queue>

#include <stack>

#include <vector>

#include <cmath>

#include <cstdio>

#include <cstdlib>

#include <cstring>

#include <iostream>

#include <algorithm>

using namespace std;

const double pi = acos(-1.0);

const int inf = 0x3f3f3f3f;

const double eps = 1e-15;

typedef long long LL;

typedef pair <int, int> PLL;

LL mod, k, b;

class MARTIX

{

    public:

        LL mat[3][3];

        MARTIX();

        MARTIX operator * (const MARTIX &b)const;

        MARTIX operator + (const MARTIX &b)const;

        MARTIX& operator = (const MARTIX &b);

}A, E, D;

MARTIX :: MARTIX()

{

    memset (mat, 0, sizeof(mat));

}

MARTIX MARTIX :: operator * (const MARTIX &b)const

{

    MARTIX ret;

    for (int i = 0; i < 2; ++i)

    {

        for (int j = 0; j < 2; ++j)

        {

            for (int k = 0; k < 2; ++k)

            {

                ret.mat[i][j] += this -> mat[i][k] * b.mat[k][j];

                ret.mat[i][j] %= mod;

            }

        }

    }

    return ret;

}

MARTIX MARTIX :: operator + (const MARTIX &b)const

{

    MARTIX ret;

    for (int i = 0; i < 2; ++i)

    {

        for (int j = 0; j < 2; ++j)

        {

            ret.mat[i][j] = this -> mat[i][j] + b.mat[i][j];

            ret.mat[i][j] %= mod;

        }

    }

    return ret;

}

MARTIX& MARTIX :: operator = (const MARTIX &b)

{

    for (int i = 0; i < 2; ++i)

    {

        for (int j = 0; j < 2; ++j)

        {

            this -> mat[i][j] = b.mat[i][j];

        }

    }

    return *this;

}

MARTIX fastpow(MARTIX ret, LL n)

{

    MARTIX ans;

    ans.mat[0][0] = ans.mat[1][1] = 1;

    while (n)

    {

        if (n & 1)

        {

            ans = ans * ret;

        }

        n >>= 1;

        ret = ret * ret;

    }

    return ans;

}

void Debug(MARTIX A)

{

    for (int i = 0; i < 2; ++i)

    {

        for (int j = 0; j < 2; ++j)

        {

            printf("%lld ", A.mat[i][j]);

        }

        printf("\n");

    }

}

MARTIX binseach(LL n)

{

    if (n == 1)

    {

        return D;

    }

    MARTIX nxt = binseach(n >> 1);

    MARTIX B = fastpow(D, n / 2);

    B = B + E;

    nxt = nxt * B;

    if (n & 1)

    {

        MARTIX C = fastpow(D, n);

        nxt = nxt + C;

    }

    return nxt;

}

int main()

{

    LL n;

    E.mat[0][0] = E.mat[1][1] = 1;

    A.mat[0][0] = A.mat[0][1] = A.mat[1][0] = 1;

//  Debug(A);

    while (~scanf("%lld%lld%lld%lld", &k, &b, &n, &mod))

    {

        if (n == 1)

        {

            MARTIX x = fastpow(A, b);

            printf("%lld\n", x.mat[0][1]);

            continue;

        }

        D = fastpow(A, k);

        MARTIX ans = binseach(n - 1);

        ans = ans + E;

        MARTIX base = fastpow(A, b);

        ans = base * ans;

//      Debug(ans);

        printf("%lld\n", ans.mat[0][1]);

    }

    return 0;

}

版权声明：本文博主原创文章，博客，未经同意不得转载。