HDU 3802 【二阶递推】

最新推荐文章于 2017-10-12 19:54:00 发布

beginendzrq

最新推荐文章于 2017-10-12 19:54:00 发布

阅读量615

点赞数

CC 4.0 BY-SA版权

分类专栏：数学其他文章标签：数列

本文链接：https://blog.youkuaiyun.com/beginendzrq/article/details/51693964

数学其他专栏收录该内容

12 篇文章

订阅专栏

本文介绍了一种求解特定类斐波拉契数列循环节的方法，并给出了具体的数学推导过程及C++实现代码。通过特征值、二次剩余等概念确定循环节长度，并利用矩阵快速幂技巧高效计算数列项。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

类斐波拉契数列的循环节
设转移矩阵

A = (a 1 b 0)

$A = \begin{pmatrix}a & b\\1 & 0\end{pmatrix}$
A 的特征值为

|λE−A| $|\lambda E - A|$ =

λ2−aλ−b $\lambda^2 - a\lambda - b$ = 0
所以 A 的特征值为

λ 1, 2 = a \pm a 2 + 4 b - - - - - - \sqrt 2

$\lambda_{1,2} = \frac {a \pm \sqrt {a^2 + 4b}}2$
设

c=a2+4b $c = a^2 + 4b$
(1) c 是模 p 的二次剩余，n = p - 1
(2) c 是模 p 的非二次剩余，n = (p + 1)(p - 1)

所以枚举 n 的因子就可以找到最小循环节了

题目大意：求 $G(a,b,n,p) = (a^{\frac {p-1}2}+1)(b^{\frac{p-1}2}+1)[(\sqrt a + \sqrt b)^{2F_n} + (\sqrt a - \sqrt b)^{2F_n}] (mod \; p)$
原式 = $(a^{\frac {p-1}2}+1)(b^{\frac{p-1}2}+1)[(a + b + \sqrt {ab})^{F_n} + (a + b -\sqrt {ab})^{F_n}] (mod \; p)$
前一部分可以快速幂计算，若不为零则可以知道是二次剩余，p - 1为循环节,所以 $F_n$ 可以对 $p - 1$ 取模
后面的式子可以化简得到递推式

G n = 2 (a + b) G n - 1 - (a - b) 2 G n - 2

$G_n = 2(a+b)G_{n-1} - (a - b)^2G_{n-2}$

#include<iostream>
#include<algorithm>
#include<cstdlib>
#include<cstdio>
#include<string>
#include<cstring>
#include<ctime>
#include<cmath>
#define LL long long
using namespace std;

LL T,a,b,n,p;
LL ret[2][2],A[2][2],t[2][2];

LL ksm(LL a,LL b)
{
    LL ret = 1;
    for (;b;b >>= 1,a = a * a % p)
        if (b & 1) ret = ret * a % p;
    return ret;
}

LL mul(LL a[][2],LL b[][2],LL p)
{
    for (int i = 0;i < 2;i ++)
        for (int j = 0;j < 2;j ++)
        {
            t[i][j] = 0;
            for (int k = 0;k < 2;k ++)
                (t[i][j] += a[i][k] * b[k][j]) %= p;
        }
    for (int i = 0;i < 2;i ++)
        for (int j = 0;j < 2;j ++)
            a[i][j] = t[i][j];
}

LL fib(LL n)
{
    n --;
    ret[0][0] = ret[1][1] = A[1][0] = A[0][0] = A[0][1] = 1;
    ret[0][1] = ret[1][0] = A[1][1] = 0;
    for (;n;n >>= 1,mul(A,A,p-1))
        if (n & 1) mul(ret,A,p-1);
    return (ret[0][0] + ret[0][1]) % (p - 1);
}

LL mat_ksm(LL n)
{
    if (!n) n = p - 1;n --;
    ret[0][0] = ret[1][1] = A[1][0] = 1;
    ret[0][1] = ret[1][0] = A[1][1] = 0;
    A[0][0] = (a+b)*2%p,A[0][1] = -(a - b) * (a - b) % p;
    for (;n;n >>= 1,mul(A,A,p))
        if (n & 1) mul(ret,A,p);
    return ((a+b)*2*ret[0][0]+ret[0][1]*2)%p;
}

int main()
{
    cin >> T;
    while (T --)
    {
        cin >> a >> b >> n >> p;
        LL ans = (ksm(a,p >> 1) + 1) * (ksm(b,p >> 1) + 1) % p;
        if (ans == 0) {puts("0");continue;}
        if (n < 2) {printf("%d\n",(a+b)*2%p*ans%p);continue;}
        (ans *= mat_ksm(fib(n))) %= p;
        if (ans < 0) ans += p;
        cout << ans << endl;
    }

    return 0;
}