【HDU5895【公式转化矩阵快速幂欧拉定义】Mathematician QSC 递推数列前n平方项和（转）

最新推荐文章于 2020-10-18 20:41:24 发布

ilblue

最新推荐文章于 2020-10-18 20:41:24 发布

阅读量352

点赞数

分类专栏：矩阵快速幂题解

题解同时被 2 个专栏收录

23 篇文章 0 订阅

订阅专栏

矩阵快速幂

5 篇文章 0 订阅

订阅专栏

给了我一种构造矩阵的新思路，值得学习哦。

Mathematician QSC

Time Limit: 2000/1000 MS (Java/Others) Memory Limit: 131072/131072 K (Java/Others)
Total Submission(s): 372 Accepted Submission(s): 196

Problem Description

QSC dream of becoming a mathematician, he believes that everything in this world has a mathematical law.

Through unremitting efforts, one day he finally found the QSC sequence, it is a very magical sequence, can be calculated by a series of calculations to predict the results of a course of a semester of a student.

This sequence is such like that, first of all,

f(0)=0,f(1)=1,f(n)=f(n−2)+2∗f(n−1)(n≥2) Then the definition of the QSC sequence is

g(n)=∑ni=0f(i)2 . If we know the birthday of the student is n, the year at the beginning of the semester is y, the course number x and the course total score s, then the forecast mark is

xg(n∗y)%(s+1) .
QSC sequence published caused a sensation, after a number of students to find out the results of the prediction is very accurate, the shortcoming is the complex calculation. As clever as you are, can you write a program to predict the mark?

Input

First line is an integer T(1≤T≤1000).

The next T lines were given n, y, x, s, respectively.

n、x is 8 bits decimal integer, for example, 00001234.

y is 4 bits decimal integer, for example, 1234.
n、x、y are not negetive.

1≤s≤100000000

Output

For each test case the output is only one integer number ans in a line.

Sample Input

  
  
   
   2
20160830 2016 12345678 666
20101010 2014 03030303 333

Sample Output

Source

2016 ACM/ICPC Asia Regional Shenyang Online

    #include<stdio.h>  
    #include<iostream>  
    #include<string.h>  
    #include<string>  
    #include<ctype.h>  
    #include<math.h>  
    #include<set>  
    #include<map>  
    #include<vector>  
    #include<queue>  
    #include<bitset>  
    #include<algorithm>  
    #include<time.h>  
    using namespace std;  
    void fre() { freopen("c://test//input.in", "r", stdin); freopen("c://test//output.out", "w", stdout); }  
    #define MS(x,y) memset(x,y,sizeof(x))  
    #define MC(x,y) memcpy(x,y,sizeof(x))  
    #define ls o<<1  
    #define rs o<<1|1  
    typedef long long LL;  
    typedef unsigned long long UL;  
    typedef unsigned int UI;  
    template <class T1, class T2>inline void gmax(T1 &a, T2 b) { if (b>a)a = b; }  
    template <class T1, class T2>inline void gmin(T1 &a, T2 b) { if (b<a)a = b; }  
    const int N = 0, M = 0, inf = 0x3f3f3f3f;  
    int casenum, casei;  
    int n, y, x, s;  
    /* 
    也可以维护！ 
    于是原始矩阵为1*4矩阵，4项分别为F[i-2] F[i-1] f(i-2)f(i-1) g(i-1) 
    于是转移矩阵为4*4矩阵， 
    0   1   0   1 
    1   4   2   4 
    0   4   1   4 
    0   0   0   1 
    */  
    const int G = 4;  
    int Z;  
    struct MX  
    {  
        int v[G][G];  
        void O() { MS(v, 0); }  
        void E() { MS(v, 0); for (int i = 0; i<G; i++)v[i][i] = 1; }  
        MX operator * (const MX &b) const  
        {  
            MX c; c.O();  
            for (int i = 0; i<G; i++)  
            {  
                for (int j = 0; j<G; j++)  
                {  
                    for (int k = 0; k<G; k++)  
                    {  
                        c.v[i][j] = (c.v[i][j] + (LL)v[i][k] * b.v[k][j]) % Z;  
                    }  
                }  
            }  
            return c;  
        }  
        MX operator ^ (LL p) const  
        {  
            MX y; y.E();  
            MX x; MC(x.v, v);  
            while (p)  
            {  
                if (p & 1)y = y*x;  
                x = x*x;  
                p >>= 1;  
            }  
            return y;  
        }  
    }a, b;  
    int euler(int x)  
    {  
        int ret = x;  
        for (int i = 2; i*i <= x; ++i)if (x % i == 0)  
        {  
            ret -= ret / i;  
            while (x % i == 0)x /= i;  
        }  
        if (x > 1)ret -= ret / x;  
        return ret;  
    }  
    int qpow(LL x, int p)  
    {  
        LL y = 1;  
        while (p)  
        {  
            if (p & 1)y = y*x%Z;  
            x = x*x%Z;  
            p >>= 1;  
        }  
        return y;  
    }  
    int main()  
    {  
        scanf("%d", &casenum);  
        for (casei = 1; casei <= casenum; ++casei)  
        {  
            scanf("%d%d%d%d", &n, &y, &x, &s); ++s;  
            a.v[0][0] = 0; a.v[0][1] = 1; a.v[0][2] = 0; a.v[0][3] = 1;  
            b.v[0][0] = 0; b.v[0][1] = 1; b.v[0][2] = 0; b.v[0][3] = 1;  
            b.v[1][0] = 1; b.v[1][1] = 4; b.v[1][2] = 2; b.v[1][3] = 4;  
            b.v[2][0] = 0; b.v[2][1] = 4; b.v[2][2] = 1; b.v[2][3] = 4;  
            b.v[3][0] = 0; b.v[3][1] = 0; b.v[3][2] = 0; b.v[3][3] = 1;  
            LL p = (LL)n * y - 1;  
            int ans = 0;  
            if (p >= 0)  
            {  
                Z = euler(s);  
                a = a*(b ^ p);  
                p = a.v[0][3] + Z;  
                Z = s;  
                ans = qpow(x, p);  
            }  
            printf("%d\n", ans);  
        }  
        return 0;  
    }  
    /* 
    【trick&&吐槽】 
     
     
    【题意】 
    我们定义f[0]=1,f[1]=1,f[n]=f[n-2]+2f[n-1] 
    定义F[i]=f[i]*f[i] 
    定义g[i]=∑F[0~i] 
    给定（x，y） 
    给定n,y,x,s 
    让你输出x^g(n*y)%(s+1) 
     
    【类型】 
    矩阵快速幂 + 欧拉定理 
     
    【分析】 
    x∈[0,1e8]，显然x并不重要 
    s∈[1,1e8]，使得s+1可能不是素数 
    n∈[0,1e8]，y∈[0,1e4]。 
    n*y会是一个[0,1e12]范围的数，这使得我们要学会快速求g[] 
     
     
    问题的核心就是求g。 
    这里肯定要使用矩阵快速幂。 
     
    我们的目的是求g(i) 
    有f(i) = f(i-2) + 2f(i-1) 
    有f(i) * f(i) = f(i-2) * f(i-2) + 4f(i-1) * f(i-1) + 4f(i-2)f(i-1) 
     
    即我们可以搞到F[i]，F[i]=F[i-2]+4F[i-1]+4f(i-2)f(i-1) 
    我们发现，维护F[i]，还需要f(i-2)f(i-1)这个东西 
     
    于是我们就考虑维护f(x-1)f(x) 
    如果我们知道f(i-2)f(i-1)，能否推出f(i-1)f(i)呢？ 
     
    考虑公式化简—— 
    f(i-1)f(i)=f(i-1)( f(i-2)+2f(i-1) ) 
    =f(i-1)f(i-2)+2F[i-1] 
     
    也可以维护！ 
    于是原始矩阵为1*4矩阵，4项分别为F[i-2] F[i-1] f(i-2)f(i-1) g(i-1) 
    于是转移矩阵为4*4矩阵， 
    0   1   0   1 
    1   4   2   4 
    0   4   1   4 
    0   0   0   1 
     
    因为这个作为指数，底数mod(s+1)，所以这里对应要mod(euler(s+1))（最后再对指数+euler(s+1)） 
    然后做快速幂，就能得到答案。 
     
    【时间复杂度&&优化】 
    O(4^3 * log(1e12) + sqrt(s+1)) 
     
    【数据】 
    n y x s 
    x^g(n*y)%(s+1) 
     
    1 0 2 1000000006 g = 0, ans = 1 
    1 1 2 1000000006 g = 1, ans = 2 
    1 2 2 1000000006 g = 5, ans = 32 
     
    */

snowy_smile原创博客，地址：http://blog.csdn.net/snowy_smile/article/details/52619292