HDU 5895 Mathematician QSC（逆元应用+矩阵快速幂+数论知识）——2016 ACM/ICPC Asia Regional Shenyang Online

传送门

Mathematician QSC

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

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\geq 2)$Then the definition of the QSC sequence is $g(n)=\sum_{i=0}^{n}f(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
  
  

题目大意：

首先定义了一个 $f()$ 函数，$f(n)=2*f(n-1)+f(n-2),f(1)=1,f(0)=0$， 然后又定义了一个 $g()$ 函数 , $g(n)=\sum_{i=0}^nf(i)^2$,

最后让你求的是 $x^{g(n*y)}\% (s+1)$， 其中 $x,y,n,s$ 都是输入的。

解题思路：
这个题目其实不是很难，仔细分析一下还是挺简单的，首先我们一看到这个 $f()$ 函数的表达式，一定会联想到矩阵快速幂，那么我们首先将 $f(n)$

能够通过矩阵快速幂得到，具体怎么得到的呢，我现在来说一下：

首先，我们设一个矩阵 $A$ ，这就是那个转移的矩阵，也就是说需要构造一个矩阵 $A$ 使得 $(f(n-2),f(n-1))*A=(f(n-1),f(n))$ 成立，那么现在 $A$

一定是一个 $2*2$ 的矩阵，然后在根据 $f(n)$ 的递推式得到 $A$ 矩阵中的具体值，算出来是：
$\begin{matrix} 0 & 1 \\ 1 & 2 \\ \end{matrix}$

然后再来分析我们要求的式子： $x^{g(n*y)}\% (s+1)$，如果满足 $g(n*y)\ge Phi(s+1)$ 的时候我们可以将其转化为

$x^{g(x*y)\% Phi(s+1)+Phi(s+1)}\% (s+1)——(1)$ ，否则只能暴力计算，然后我们发现只有当 $n*y\ge 11$ 的时候才暴力计算也很简单了，直接求，

然后快速幂就ok了，关键是解决上面 $(1)$ 式的那个问题，首先我们看一下 $g(n)$ 能不能写成关于 $f()$ 的比较简单的式子，然后我们将其展开发现

与 $f(i)*f(i+1)$ 有关系，那么 $f(i)*f(i+1)$ 等于什么呢：

$f(n)*f(n+1)=f(n)*(2*f(n)-f(n-1))$

进而得到了：

$f(n)^2 = \frac {f(n)*f(n+1)-f(n)*f(n-1)}2$

$f(n-1)^2 = \frac {f(n-1)*f(n)-f(n-1)*f(n-2)}2$

$\ \ \ \ \ \ \ \ \ ...$

$f(1)^2 = \frac {f(1)*f(2)-f(1)*f(0)}2$

通过累加得到 $\sum_{i=1}^n=\frac{f(n)*f(n+1)}2$ 又因为当 $i==0$ 的时候 $f(0)^2=0$ 所以$g(n)=\sum_{i=0}^n=\frac{f(n)*f(n+1)}2$ ，

一看到 $g(n)$ 有除法我们就想到求逆元，可是有一个至关重要的问题是 这个式子可能没有逆元，那么怎么做呢，我们就用到了一个逆元的应用，

$\frac a b \% c$ 可以写成 $\frac {a\% (b*c)}b$， 然后就可以解决这个问题了，最后矩阵快速幂一下就行了，在需要注意的问题就是 输入的

顺序的问题，是 先输入的 $n,y,x,s$ 绝对不是 $n,x,y,s$ 我就被这个地方坑了好长时间。。。

$My\ Code：$

/**
2016 - 09 - 19 晚上
Author: ITAK

Motto:

今日的我要超越昨日的我，明日的我要胜过今日的我，
以创作出更好的代码为目标，不断地超越自己。
**/
#include <iostream>
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cmath>
#include <vector>
#include <queue>
#include <algorithm>
#include <set>
#include <time.h>
using namespace std;
typedef long long LL;
typedef unsigned long long ULL;
const LL INF = 1e9+5;
const int MAXN = 2;
const double eps = 1e-7;
const double PI = acos(-1);
using namespace std;
LL Scan_LL()///输入外挂
{
    LL res=0,ch,flag=0;
    if((ch=getchar())=='-')
        flag=1;
    else if(ch>='0'&&ch<='9')
        res=ch-'0';
    while((ch=getchar())>='0'&&ch<='9')
        res=res*10+ch-'0';
    return flag?-res:res;
}
void Out(LL a)///输出外挂
{
    if(a>9)
        Out(a/10);
    putchar(a%10+'0');
}

struct Matrix
{
    LL mat[MAXN][MAXN];
};
const Matrix I = {1, 0,
            0, 1,
           };

const Matrix P = {0, 1,
            1, 2,
           };
Matrix Matrix_Mul(Matrix a, Matrix b, LL MOD)
{
    Matrix c;
    for(int i=0; i<MAXN; i++)
        for(int j=0; j<i; j++)
            swap(b.mat[i][j], b.mat[j][i]);
    for(int i=0; i<MAXN; i++)
    {
        for(int j=0; j<MAXN; j++)
        {
            c.mat[i][j] = 0;
            for(int k=0; k<MAXN; k++)
            {
                c.mat[i][j] += a.mat[i][k] * b.mat[j][k];
                c.mat[i][j] %= MOD;
            }
        }
    }
    return c;
}
Matrix Matrix_Quick_Mod(Matrix a, LL m, LL MOD)
{
    Matrix ans = I;
    while(m)
    {
        if(m & 1)
            ans = Matrix_Mul(ans, a, MOD);
        m>>=1;
        a = Matrix_Mul(a, a, MOD);
    }
    return ans;
}
LL Quick_Mod(LL a, LL b, LL MOD)
{
    LL ans = 1LL;
    while(b)
    {
        if(b & 1)
            ans = (ans*a) % MOD;
        b>>=1;
        a = (a*a) % MOD;
    }
    return ans;
}
LL Phi(LL x)
{
    LL ans = x;
    for(LL i=2LL; i*i<=x; i++)
    {
        if(x % i == 0)
        {
            ans -= ans/i;
            while(x % i == 0)
                x /= i;
        }
    }
    if(x > 1)
        ans -= ans/x;
    return ans;
}
LL a[13] = {0, 1, 2, 5, 12, 29, 70, 169, 408, 985, 2378, 5741, 13860};
LL get_g(LL n)
{
    LL sum = 0;
    for(LL i=0; i<=n; i++)
        sum += a[i] * a[i];
    return sum;
}
int main()
{
    int T;
    scanf("%d",&T);
    while(T--)
    {
        LL n, x, y, s;
        n = Scan_LL();
        y = Scan_LL();
        x = Scan_LL();
        s = Scan_LL();
        n *= y;
        s++;
        if(n < 12LL)
        {
            printf("%I64d\n",Quick_Mod(x, get_g(n), s));
            continue;
        }
        LL phi = Phi(s);
        LL MOD1 = 2LL*phi;
        Matrix tmp1 = Matrix_Quick_Mod(P, n, MOD1);
        LL ans1 = ( ( (tmp1.mat[1][1]*tmp1.mat[1][0]) % MOD1 ) >> 1LL );
        ans1 += phi;
        LL ans = Quick_Mod(x, ans1, s);
        printf("%I64d\n",ans);
    }
    return 0;
}