HDU 5674 Function(斐波那契模数列循环节)

Problem Description
There is a function
f(n)=(a+√b)^n+(a−√b)^n. a and b are integers
(1≤a,b≤1,000,000).
Maybe the function looks complex but it is actually an integer.The question is to calculate
f(x^y).The answer can bevery large,so just output the answer mod 1,000,000,007.

Input
There are multiple test cases. The first line of input contains an integer
T(1≤T≤200)
indicating the number of test cases. For each test case:

One line contains four integers
a,b (1≤a,b≤1,000,000), x(1≤x≤50),y(1≤y≤10^18).

Output
For each test case, output one integer.

Sample Input
3
3 5 1 1
3 5 2 1
1 1 49 99999

Sample Output
6
28
160106184

斐波那契数列循环节的定理就不在这里证明了,结论就是:假设MOD为模数,则在f[2]!=MOD时,(MOD+1)(MOD-1)一定为其一个循环节.这道题f[n]明显是裴波那契数列通项的形式,所以可以利用裴波那契数列的特征方程来构造裴波那契数列 f[n]=2a*f[n-1]+(b-a^2)*f[n-2]; 之后就可以利用循环节+快速幂求出答案了.

#include <stdio.h>
#include <iostream>
#include <algorithm>
#include <string>
#include <string.h>
using namespace std;
#define ll long long 
#define F(x,a,b) for (int x=a;x<=b;x++)
#define _fast F(i,1,2)F(j,1,2)F(k,1,2) b[i][j]=((b[i][j]+(a[i][k]%MOD)*(a[k][j]%MOD))%MOD+MOD)%MOD;
#define _reset F(i,1,2)F(j,1,2) {a[i][j]=b[i][j];b[i][j]=0;}
#define _c F(i,1,2)F(j,1,2)F(k,1,2) b[i][j]=((b[i][j]+(a[i][k]%MOD)*(c[k][j]%MOD))%MOD+MOD)%MOD;
ll a[4][4],b[4][4],c[4][4];
ll f[3];
const ll MOD = 1e9+7;
void Init(ll aa,ll bb)
{
    a[1][1]=aa;a[1][2]=bb;
    a[2][1]=1;a[2][2]=0;
    for (int i=1;i<=2;i++)
        for (int j=1;j<=2;j++) c[i][j]=a[i][j];
}
ll multi(ll x,ll y,ll mod)
{
    ll ans=0;
    while (y)
    {
        if (y&1) ans=(ans+x+mod)%mod;
        x=((x<<1)+mod)%mod;
        y=y>>1;
    }
    return ans;
}
void fastmat(ll x)
{
    if (x<=1) return;
    if (x&1) {fastmat(x-1);_c _reset }
    else 
    {
        fastmat(x/2); _fast _reset
    }
}
ll _f(ll x,ll k, ll mod)
{
    if (!k) return 1;
    if (k&1) return multi(_f(x,k-1,mod),x,mod);
    else
    {
        ll t=_f(x,k/2,mod);
        return multi(t,t,mod);
    }
}
void out()
{
    for (int i=1;i<=2;i++)
    {
        for (int j=1;j<=2;j++)
            printf("%I64d ",a[i][j]);
        printf("\n");
    }
}
int main()
{
     int T;
     scanf("%d",&T);
     while (T--)
     {

         ll aa,bb,x,y;
         scanf("%I64d%I64d%I64d%I64d",&aa,&bb,&x,&y);
         f[1]=aa*2;f[2]=(aa*aa+bb)*2;
         ll a1=2*aa,b1=bb-aa*aa;
         Init(a1,b1);
         if (b1==0) 
         {
            ll tt=_f(x,y,MOD-1);
            printf("%I64d\n",_f(f[1],tt,MOD));
         }else 
         {
            ll mmod=(MOD-1)*(MOD+1);
            ll tt=_f(x,y,mmod);
            if (!tt)
            {
                fastmat(mmod-2);
                printf("%I64d\n",(a[1][1]%MOD*f[2]%MOD+a[1][2]%MOD*f[1]%MOD)%MOD);
            }else 
            {
                if (tt<=2){printf("%I64d\n",f[tt]%MOD);}
                else 
                {
                  fastmat(tt-2);
                  printf("%I64d\n",((a[1][1]%MOD)*(f[2]%MOD)+(a[1][2]%MOD)*(f[1]%MOD))%MOD);
                  //out();  
                }
            }
         }


     }


}