Interesting Fibonacci HDU 2814 （斐波那契数列模m性质+欧拉函数降幂（指数循环节）详解）

Interesting Fibonacci

In mathematics, the Fibonacci numbers are a sequence of numbers named after Leonardo of Pisa, known as Fibonacci (a contraction of filius Bonaccio, “son of Bonaccio”). Fibonacci’s 1202 book Liber Abaci introduced the sequence to Western European mathematics, although the sequence had been previously described in Indian mathematics.
The first number of the sequence is 0, the second number is 1, and each subsequent number is equal to the sum of the previous two numbers of the sequence itself, yielding the sequence 0, 1, 1, 2, 3, 5, 8, etc. In mathematical terms, it is defined by the following recurrence relation:

That is, after two starting values, each number is the sum of the two preceding numbers. The first Fibonacci numbers (sequence A000045 in OEIS), also denoted as F[n];
F[n] can be calculate exactly by the following two expressions:

A Fibonacci spiral created by drawing arcs connecting the opposite corners of squares in the Fibonacci tiling; this one uses squares of sizes 1, 1, 2, 3, 5, 8, 13, 21, and 34;

So you can see how interesting the Fibonacci number is.
Now AekdyCoin denote a function G(n)

Now your task is quite easy, just help AekdyCoin to calculate the value of G (n) mod C 

Input
The input consists of T test cases. The number of test cases (T is given in the first line of the input. Each test case begins with a line containing A, B, N, C (10<=A, B<2^64, 2<=N<2^64, 1<=C<=300)
Output
For each test case, print a line containing the test case number( beginning with 1) followed by a integer which is the value of G(N) mod C
Sample Input

1
17 18446744073709551615 1998 139

Sample Output

Case 1: 120

题意：

$G(1) = F(a^b)$

$G(n) = G(n-1)^{F(a^b)} \ (n \ge 2)$

求$G(n)\ mod \ c$

分析：

首先我们观察所给的函数

$G(1) = F(a^b) \ \ \ \% \ c$

$G(2) = G(1)^{F(a^b)} = F(a^b)^{F(a^b)} \ \ \ \% \ c$

$G(3) = G(2)^{F(a^b)} = (F(a^b)^{F(a^b)})^{F(a^b)} = F(a^b)^{F(a^b)^2} \ \ \ \% \ c$
.
.
.
$G(n) = F(a^b)^{F(a^b)^{n-1}} \ \ \ \% \ c$

所以我们找到了最终需要求的公式

$G(n) = F(a^b)^{F(a^b)^{n-1}} \ \ \ \% \ c$

但是我们发现每个数可能都会很大，不要着急我们将问题分解一步一步看

---

一、首先将$指数 F(a^b)^{n-1}$看成一个整体，我们设

$e = F(a^b)^{n-1}$

因此原式变成了

$G(n) = F(a^b)^e \ \ \ \% \ c$

我们先看底数部分$F(a^b)\ \ \%\ \ c$

按照一般思路必然要先求出$a^b$,然后再去求斐波那契数列，可以递推，可以公式，可以矩阵快速幂，但是毫无疑问，$a^b$太大了，先不看计算斐波那契数在时间上的花费问题，因为我们要求第$a^b$个斐波那契数，因此$a^b$不能取模，所以必然会long long，所以直接求是不可能的。

斐波那契序列模m
如果对斐波那契数列中的数模m简化会怎样呢？
先举几个简单的例子：
$F_n\ (mod\ 2)$ 1,1,0, 1,1,0,1,1,0,1,1,0…
$F_n\ (mod\ 3)$ 1,1,2,0,2,2,1,0, 1,1,2,0,2,2,1,0,1,1,2,0,2,2,1,0…
$F_n\ (mod\ 4)$ 1,1,2,3,1,0, 1,1,2,3,1,0,1,1,2,3,1,0,1,1,2,3,1,0…
$F_n\ (mod\ 5)$ 1,1,2,3,0,3,3,1,4,0,4,4,3,2,0,2,2,4,1,0,1,1,2,3,0,3,3…
于是我们发现：
对于计算斐波那契数列模m时，存在整数$N \ge 1$使得
$F_{n+N} \equiv F_n \ \ (mod \ m)$ 对所有的$n = 1,2...$

根据上面我们发现的斐波那契序列模m的性质（注意只有在取模时才有这种性质）

回到我们的问题上来，现在想求啥来？

$F(a^b)\ \ \%\ \ c$

所以我们完全可以先通过暴力找出斐波那契序列模c的循环节，及其长度len1，这样我们只需要算出$a^b$在第一个循环节中对应的位置即可

即求处$a^b \ \% \ c$快速幂取模很好解决吧

这样我们求出了斐波那契数$F(a^b)在模c情况下的值，我们记为f_1$

在来看看现在变成求什么了

求$G(n) = (f_1)^e \ \ \ \% \ c$

---

二、这时我们关注指数部分$e$怎么求

在上一部分我们设

$e = F(a^b)^{n-1}$

因为指数中含有$F(a^b)^{n-1}$，如果不取模的话，很显然，我们依然算不出来，因为他太大了，必然会超long long.

这时需要用到第二个知识点

指数循环节
使用欧拉函数可对幂模运算进行降幂
它实际上是费马小定理降幂的推广

公式：

$a^b \ \%\ c = a^{b\ \%\ \phi(c) + \phi(c)} \ \%\ c$

因此

求$G(n) = (f_1)^e \ \ \ \% \ c$

可以转化为

$G(n) = (f_1)^{e\ \% \ \phi(c) + \phi(c)} \ \ \ \% \ c$

所以我们得求

$e\ \% \ \phi(c) + \phi(c) = F(a^b)^{n-1}\ \% \ \phi(c) + \phi(c)$

发现什么没有，是不是和我们求的第一步的步骤很类似

没错类比上面我们需要先求出模$\phi(c)$下的斐波那契序列的循环节及其长度len2，

这样我们可以得到$a^b$在其第一个循环节中所处的位置

即求出$a^b \ \% \ \phi(c)$同样快速幂取模可求

这样就可以求得了模$\phi(c)$下的斐波那契数$F(a^b)$，我们记为$f_2$

---

三、求：

$G(n) = (f_1)^{e\ \% \ \phi(c) + \phi(c)} \ \ \ \% \ c$

变成了求

$G(n) = (f_1)^{(f_2)^{n-1}\ \% \ \phi(c) + \phi(c)} \ \ \ \% \ c$

现在只需要求出

$e = (f_2)^{n-1}\ \% \ \phi(c) + \phi(c)$

直接快速幂取模求解没问题了吧

这样也求出了e

我们就可以求最终答案$G(n) = (f_1)^e \ \ \ \% \ c$了

---

注意问题（坑）：

1）数据范围是$2^{64}-1$，所以要用unsigned long long

2） c=1的时候，求斐波那契序列的循环节是无效的，因为任何数模1都是0，因此如果c=1，则说明斐波那契数列全是0，直接输出0即可

3）当$\phi(c) = 1$的时候，有模$\phi(c)$的斐波那契数列全是0，即指数部分$F(a^b) \ \%\ \phi(c)$为零，所以$e = F(a^b)^{n-1} \ \%\ \phi(c) + \phi(c) = 1$
故答案为第一步中求出的$f_1$

4）求快速幂的时候，因为a可能很大，第一次乘就有可能超了long long，因此在运算之前先对a取模，然后再进行快速幂运算。

code：

#include <iostream>
#include <cstdio>
#include <cstring>
using namespace std;
typedef unsigned long long ll;
ll F[3000];
int get_phi(int n){
    int ans = n;
    for(int i = 2; i * i <= n; i++){
        if(n % i == 0){
            ans = ans / i * (i - 1);
            while(n % i == 0) n /= i;
        }
    }
    if(n != 1) ans = ans / n * (n - 1);
    return ans;
}
ll q_pow(ll a,ll b,ll mod){
    ll ans = 1;
    a %= mod;
    while(b){
        if(b & 1)
            ans = ans * a % mod;
        b >>= 1;
        a = a * a % mod;
    }
    return ans;
}
int get_loop(int m){
    F[0] = 0;
    F[1] = F[2] = 1;
    for(int i = 3; ; i++){
        F[i] = (F[i-1] % m + F[i-2] % m) % m;
        if(F[i] == 1 && F[i-1] == 0) return i - 1;
    }
}
int main(){
    ll a,b,n,c;
    int T,cas = 0;
    scanf("%d",&T);
    while(T--){
        scanf("%llu%llu%llu%llu",&a,&b,&n,&c);
        if(c == 1){
            printf("Case %d: 0\n",++cas);
            continue;
        }
        ll phic = get_phi(c);
        int len1 = get_loop(c);
        ll f1 = q_pow(a,b,len1);
        f1 = F[f1];
        if(phic == 1){
            printf("Case %d: %llu\n",++cas,f1);
            continue;
        }
        int len2 = get_loop(phic);
        ll f2 = q_pow(a,b,len2);
        f2 = F[f2];
        ll e = q_pow(f2,n-1,phic) + phic;
        ll ans = q_pow(f1,e,c);
        printf("Case %d: %llu\n",++cas,ans);
    }
    return 0;
}