【HDU5895 2016 ACM ICPC Asia Regional Shenyang Online D】【公式转化 矩阵快速幂 欧拉定义】Mathematician QSC 递推数列前n平方项和

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

  
  

   
   1
317
  
  

 

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

*/