题解 P5175 数列

· · 题解

更好的阅读体验

CTime_Pup_314的博客

P5175 数列

求二阶线性常系数齐次递推的平方和

蒟蒻不会构造矩阵，所以用了一种很暴力的做法

由于这道题意要求，下文数列下标都从 1 开始

前置知识

线性常系数齐次递推关系

若递推关系满足

f_n\ =\ \sum_{i=1}^kc_if_{n-i}\ (n >k)

且对于 f_1,\ c_1,\ f_2,\ c_2\ ...\ f_k,\ c_k 都为常数

则称这种递推关系为线性常系数齐次递推关系，其中 k 为阶数

斐波那契数列就是常见的满足二阶线性常系数齐次递推关系的数列

特征方程

求解线性常系数齐次递推方法的有很多，特征方程是常见的一种

对于递推 f_n\ =\ 2f_{n-1}+3f_{n-2} 来说

我们构造等比数列 \langle\ 1,\ -1,\ 1\ ...\ \rangle 满足这种关系，其中 g_n\ =\ (-1)^{n-1}

又另一个构造等比数列 \langle\ 1,\ 3,\ 9\ ...\ \rangle 满足这种关系，其中 h_n\ =\ 3^{n-1}

仔细思考发现对于常数 \alpha 和 \beta 而言 \lbrace \alpha g_n\rbrace 和 \lbrace \beta h_n\rbrace 都满足

进一步思考后感觉对于 \lbrace \alpha g_n\ +\ \beta h_n\rbrace 也都满足，其正确性十分明显无需多言

也就是说我们可以通过一些等比数列作为基底，通过它们的线性组合，生成很多满足线性常系数齐次递推的数列

如何找到这些等比数列呢？拿二阶来说

我们不妨设等比数列 q_n 的公比为 x，我们发现递推关系可以写成

x^2q_n\ =\ axq_n+\ bq_n

两边同时消去 q_n，只剩下

x^2\ =\ ax\ + b

这个方程叫做特征方程，通过这个我们就可以找到所有满足条件的等比数列

现在用斐波那契数列举例，特征方程为

x^2\ =\ x\ +\ 1

解得 x_1\ =\ \frac{1+\sqrt{5}}{2}，x_2\ =\ \frac{1-\sqrt{5}}{2}

我们现在要解出 \alpha 和 \beta，则有 f_n= \alpha(x_1)^{n-1}+\beta(x_2)^{n-1}

带入 f_1=1 和 f_2=1 得

\begin{cases}\alpha+\beta=1\\\\\alpha x_1+\beta x_2=1\\\end{cases}\Rightarrow\begin{cases}\alpha=\frac{\sqrt{5}-1}{2\sqrt{5}}\\\\\beta=\frac{\sqrt{5}+1}{2\sqrt{5}}\\\end{cases}

所以 f_n=\frac{1}{\sqrt{5}}[(\frac{1+\sqrt{5}}{2})^n-\frac{1-\sqrt{5}}{2})^n]

模意义下的根式运算

假如我们在计算过程中出现了根式，但保证根号下的数不变且结果一定不包含根式，要求在模P意义下进行，保证计算过程中的任何数有逆元，如何解决呢？

我们可以把根式写作类似复数的形式 a+\sqrt{k}b 其中 k 是所有数的根号下面的部分，如 a+\sqrt{2}\ b 等

我们把 a 称作有理部，b 称作无理部，\sqrt{k} 称作基

下面给出几项重要的运算，基不会被取模

加法 (a+\sqrt{k}b)+(c+\sqrt{k}d)\ \equiv\ (a+b)+\sqrt{k}(c+d)\ (mod\ P)

减法 (a-\sqrt{k}b)+(c-\sqrt{k}d)\ \equiv\ (a-b)+\sqrt{k}(c-d)\ (mod\ P)

共轭 conj(a+\sqrt{k}b)\ \equiv\ a+\sqrt{k}(-b)\ (mod\ P)

乘法 (a+\sqrt{k}b)(c+\sqrt{k}d)\ \equiv\ (ac+bdk)+\sqrt{k}(ad+bc)\ (mod\ P)

除法比较特殊，我们需要先将分母有理化，即同时将分子和分母同时乘上分母的共轭，这时分母就只剩下有理部了，再把分子部分除以分母有理部的逆元即可

除法 \frac{a+\sqrt{k}b}{c+\sqrt{k}d}\ \equiv\ \frac{(a+\sqrt{k}b)(c+\sqrt{k}d)}{c^2-kd^2} (mod\ P)

写出来就是

struct Complex
{
    int rat, irr; // rat 有理部 irr 无理部 irr_base 基
    Complex(int rat = 0, int irr = 0):rat(rat), irr(irr) {};
    Complex conj() const { return Complex(rat, P-irr); }
    Complex operator - (const Complex &_) const { return Complex((rat-_.rat+P)%P, (irr-_.irr+P)%P); }
    Complex operator + (const Complex &_) const { return Complex((rat+_.rat)%P, (irr+_.irr)%P); }   
    Complex operator * (const Complex &_) const
    {
        Complex ret; 
        ret.rat = (rat*_.rat%P+irr*_.irr%P*irr_base%P)%P;
        ret.irr = (rat*_.irr%P+irr*_.rat%P)%P;
        return ret;
    }
    Complex operator / (const Complex &_) const
    {
        Complex a = Complex(rat, irr)*_.conj(); Complex b = _*_.conj();
        a = a*Complex(qpow(b.rat, P-2), 0); // 求逆
        return a;
    }
};

题目主体

题干给出的是一个二阶线性常系数齐次递推数列，而让求出前 n 项平方和

有前置知识得出这个数列的通项的是几个等比数列的线性组合，那么最后也无非是几个等比数列求和罢了

同时这道题保证了特征方程有两根和根式模运算的前提条件，所以下面做法极其暴力

设数列为 \lbrace f_n\rbrace 其他变量名含义与前置知识中的相同

则有

\begin{cases} x_1\ =\ \frac{a+\sqrt{a^2+4b}}{2}\\\\ x_2\ =\ \frac{a-\sqrt{a^2+4b}}{2}\\\\ \alpha\ =\ \frac{x_1f_1-f_2}{x_1-x_2}\\\\ \beta\ =\ \frac{f_2-x_2f_1}{x_1-x_2} \end{cases}

而 f_n= \alpha(x_1)^{n-1}+\beta(x_2)^{n-1} 那么 f_n^2= \alpha^2(x_1^2)^{n-1}+\beta^2(x_2^2)^{n-1}+2\alpha\beta(x_1x_2)^{n-1}

每一部分都是等比数列，用公式 S_n\ =\ \frac{x^n-1}{x-1} 即可

暴力代码如下，没有加任何常数优化居然跑得最快

#include <cctype>
#include <algorithm>
#include <cstdio>
#define int long long // 虽然习惯不好，但是在这题比较方便
using namespace std;
const int P = 1e9+7;
int n, t, a, b, irr_base, f1, f2;
int qpow(int a, int b)
{
    int ret = 1;
    for( ; b; b >>= 1, a = a*a%P) if(b&1) ret = ret*a%P;
    return ret;
}
struct Complex
{
    int rat, irr;
    Complex(int rat = 0, int irr = 0):rat(rat), irr(irr) {};
    Complex conj() const { return Complex(rat, P-irr); }
    Complex operator - (const Complex &_) const { return Complex((rat-_.rat+P)%P, (irr-_.irr+P)%P); }
    Complex operator + (const Complex &_) const { return Complex((rat+_.rat)%P, (irr+_.irr)%P); }   
    Complex operator * (const Complex &_) const
    {
        Complex ret; 
        ret.rat = (rat*_.rat%P+irr*_.irr%P*irr_base%P)%P;
        ret.irr = (rat*_.irr%P+irr*_.rat%P)%P;
        return ret;
    }
    Complex operator / (const Complex &_) const
    {
        Complex a = Complex(rat, irr)*_.conj(); Complex b = _*_.conj();
        a = a*Complex(qpow(b.rat, P-2), 0);
        return a;
    }
};
Complex qpow(Complex a, int b)
{
    Complex ret = Complex(1, 0);
    for( ; b; b >>= 1, a = a*a) if(b&1) ret = ret*a;
    return ret;
}
Complex alpha, beta, x1, x2, fn;
Complex sum1, sum2, sum3, ans;
signed main()
{
    scanf("%lld", &t);
    while(t--)
    {
        scanf("%lld%lld%lld%lld%lld", &n, &f1, &f2, &a, &b);
        irr_base = (a*a%P+4*b%P)%P;
        x1 = Complex(a, 1)/Complex(2, 0);
        x2 = Complex(a, P-1)/Complex(2, 0);
        alpha = (Complex(f2, 0)-x2*Complex(f1, 0))/(x1-x2);
        beta = (x1*Complex(f1, 0)-Complex(f2, 0))/(x1-x2);
        fn = alpha*qpow(x1, n-1)+beta*qpow(x2, n-1);
        sum1 = (qpow(x1*x1, n)-Complex(1, 0))/(x1*x1-Complex(1, 0))*alpha*alpha;
        sum2 = (qpow(x2*x2, n)-Complex(1, 0))/(x2*x2-Complex(1, 0))*beta*beta;
        sum3 = (qpow(x1*x2, n)-Complex(1, 0))/(x1*x2-Complex(1, 0))*beta*alpha*Complex(2, 0);
        ans = sum1+sum2+sum3;
        printf("%lld\n", ans.rat);
    }

    return 0;
}