CTime_Pup_314

不能再颓废了

0%

数论出题组比赛 数列

发表于 更新于 分类于 阅读次数： Valine：

题目描述

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

题解

由于这道题意要求，下文数列下标都从 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)\)

写出来就是

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
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}\) 即可

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
#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;
}