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?
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
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
#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
*/