[原创]HDU 5895&&2016 ACM/ICPC Asia Regional Shenyang Online1004 Mathematician QSC [矩阵加速+欧拉降幂]【数论】

2016-09-19 20:02:35 Tabris_ 阅读数：454

博客爬取于2020-06-14 22:43:15
以下为正文

版权声明：本文为Tabris原创文章，未经博主允许不得私自转载。
https://blog.csdn.net/qq_33184171/article/details/52588817

题目链接：http://acm.hdu.edu.cn/showproblem.php?pid=5895
------------------------.
Mathematician QSC

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

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 x^g(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
1
317

Source
2016 ACM/ICPC Asia Regional Shenyang Online

------------------------.

题目大意：
就是给你四个数n,y,x,s,
让你求x^g(n∗y)%(s+1).
其中g(n)=∑(i->n)f(i)^2;
f(0)=0,f(1)=1,f(n)=f(n−2)+2∗f(n−1)(n≥2)

解题思路：
对于x的指数g(n) 是一个很大的数 所以需要想办法把它改成我们能计算的 就是欧拉降幂

然后只要求解g(n)就行了
很容易的想到g(n)=f(n)*f(n+1)/2 其中f(n)只用一个矩阵加速就能很快地求解
但是随后发现这样并不行 因为f(n)已经是对Phi(s+1)取模之后的数了 在/2之后之就会出错 然后想到求逆元的办法解决 但是随后发现虽然2是一个质数但是并不能满足gcd(2,s+1)==1 因为s+1%2可能等于0 于是这个思路就GG了。。。

最后还是看了ICPCCamp的题解 才知道解法
最终还是一个矩阵快速幂
（f[n]^2,f[n+1]^2,f[n]*f[n+1],g[n]）
↑这是左矩阵
[0,1,0,0]
[1,4,2,1]
[0,4,1,0]
[0,0,0,1] ←这是右矩阵

是这么解释的
这里写图片描述
相信你已经看懂了

//
就是思维太局限了 首先欧拉降幂不知道
再后来矩阵不会构造 导致这场的GG。。
//

附本题代码
--------------------------------------.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
# include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int maxn = 505;
const double  Pi =  acos(-1);
# define pb push_back
# define lalal puts("****");
const int M = 4;
int MOD ;
struct Matrix
{
    LL m[M][M];
    void clearO()
    {
        for(int i=0; i<M; i++) //初始化矩阵
            for(int j=0; j<M; j++)
                m[i][j]= 0;
    }
    void clearE()
    {
        for(int i=0; i<M; i++) //初始化矩阵
            for(int j=0; j<M; j++)
                m[i][j]= (i==j);
    }
    void display()
    {
        for(int i=0; i<M; i++)
            {
                for(int j=0; j<M; j++)
                printf("%d ",m[i][j]);
                puts("");
            }
    }
};

Matrix operator * (Matrix a,Matrix b)
{
    Matrix c;
    c.clearO();

    for(int k=0; k<M; k++)
        for(int i=0; i<M; i++) //实现矩阵乘法
        {
            if(a.m[i][k] <= 0)  continue;
            for(int j=0; j<M; j++)
            {
                if(b.m[k][j] <= 0)    continue;
                c.m[i][j]=(c.m[i][j]+a.m[i][k]*b.m[k][j]+MOD)%MOD;
            }
        }
    return c;
}

Matrix operator ^ (Matrix a,LL b)
{
    Matrix c;
    c.clearE();
    while(b)
    {
        if(b&1) c= c * a ;
        b >>= 1;
        a = a * a ;
    }
    return c;
}

int Is_or[101001];
int prime[13000],kpri;
void Prime()
{
    int n=100001;
    kpri=0;
    memset(Is_or,1,sizeof(Is_or));
    Is_or[0]=Is_or[1]=0;
    for(int i=2; i<n; i++)
    {
        if(Is_or[i])
        {
            prime[kpri++]=i;
            for(int j=i+i; j<n; j+=i)
            {
                Is_or[j]=0;
            }
        }
    }
    //prime[kpri]=10007;
    //printf("%d\n",kpri);
    return ;
}

LL Phi(LL n)
{
    LL rea=n;
    for(int i=0; prime[i]*prime[i]<=n; i++)
    {
        if(n%prime[i]==0)
        {
            rea=rea-rea/prime[i];
            while(n%prime[i]==0)
                n/=prime[i];
        }
    }
    if(n>1)  rea=rea-rea/n;
    return rea;
}

LL qmod(LL a,LL b)
{
    LL res= 1;
    while(b)
    {
        if(b&1) res=(res*a)%MOD;
        b>>=1;
        a=(a*a)%MOD;
    }
    return res;
}

int main()
{
    Prime();
    int _;
    while(~scanf("%d",&_))
    {
        while(_--)
        {
            LL n,y,x,s;
            scanf("%I64d%I64d%I64d%I64d",&n,&y,&x,&s);

            MOD = Phi(s+1);

            Matrix a,b;
            a.clearO(),b.clearO();
            a.m[0][1]=1;

            b.m[0][1]=1;
            b.m[1][0]=1,b.m[1][1]=4,b.m[1][2]=2,b.m[1][3]=1;
            b.m[2][1]=4,b.m[2][2]=1;
            b.m[3][3]=1;

            b=b^(n*y);
            a=a*b;

            LL zhi = a.m[0][3]%MOD+MOD;

            MOD=s+1;
            printf("%I64d\n",qmod(x%MOD,zhi));
        }
    }
    return 0;
}