题解 P5517 【[MtOI2019]幻想乡数学竞赛】

· · 题解

真·数学&卡常题

官方题解中似乎是用生成函数，我这里就讲另一个神奇的科技——特征方程。

首先，我们注意到这个递推式中有个很讨厌的3^n，如果把这个给去掉就好了。注意到，每一项都带一个3^n，我们考虑把这个东西求出来，把递推式边个形，有

3^n=a_n-3a_{n-1}-a_{n-2}+3a_{n-3}

即

3^{n-1}=a_{n-1}-3a_{n-2}-a_{n-3}+3a_{n-4}

那么有

\begin{aligned}3^n&=3\cdot3^{n-1}\\&=3a_{n-1}-9a_{n-2}-3a_{n-3}+9a_{n-4}\end{aligned}

所以

\begin{aligned}a_n&=3a_{n-1}+a_{n-2}-3a_{n-3}+3^n\\&=3a_{n-1}+a_{n-2}-3a_{n-3}+3a_{n-1}-9a_{n-2}-3a_{n-3}+9a_{n-4}\\&=6a_{n-1}-8a_{n-2}-6a_{n-3}+9a_{n-4}\end{aligned}

这样就变成了一个齐次的递推式了。顺便一提，这是四阶的，所以要先知道前四项，用最初的递推式手算一下得到a_3=-6。

现在来求这个东西的通项公式。先引入一个东西——特征方程。顺便再介绍一下，这个东西很有用，出了在解递推数列中有用，在矩阵、微分方程、积分方程等中都是一个黑科技。下面就介绍一下齐次递推数列中的特征方程。

首先，我们把递推式写成这样的形式：

a_{n+4}=6a_{n+3}-8a_{n+2}-6a_{n+1}+9a_n

然后，把a_{n+k}替换成x^k，得到

x^4=6x^3-8x^2-6x+9

这有什么好处呢？不妨设x_1,x_2,x_3,x_4为它的四个根，那么就有

\begin{cases}x_1^4=6x_1^3-8x_1^2-6x_1+9\\x_2^4=6x_2^3-8x_2^2-6x_2+9\\x_3^4=6x_3^3-8x_3^2-6x_3+9\\x_4^4=6x_4^3-8x_4^2-6x_4+9\end{cases}

两边同时乘以x_k^n，得

\begin{cases}x_1^{n+4}=6x_1^{n+3}-8x_1^{n+2}-6x_1^{n+1}+9x_1^n&(1)\\x_2^{n+4}=6x_2^{n+3}-8x_2^{n+2}-6x_2^{n+1}+9x_2^n&(2)\\x_3^{n+4}=6x_3^{n+3}-8x_3^{n+2}-6x_3^{n+1}+9x_3^n&(3)\\x_4^{n+4}=6x_4^{n+3}-8x_4^{n+2}-6x_4^{n+1}+9x_4^n&(4)\end{cases}

那么，把(1)式乘以c_1（一个系数），把(2)式乘以c_2，把(3)式乘以c_3，把(4)式乘以c_4，然后相加得到

\begin{aligned}&c_1x_1^{n+4}+c_2x_2^{n+4}+c_3x_3^{n+4}+c_4x_4^{n+4}\\=&6(c_1x_1^{n+3}+c_2x_2^{n+3}+c_3x_3^{n+3}+c_4x_4^{n+3})\\-&8(c_1x_1^{n+2}+c_2x_2^{n+2}+c_3x_3^{n+2}+c_4x_4^{n+2})\\-&6(c_1x_1^{n+1}+c_2x_2^{n+1}+c_3x_3^{n+1}+c_4x_4^{n+1})\\+&8(c_1x_1^n+c_2x_2^n+c_3x_3^n+c_4x_4^n)\end{aligned}

当a_n=c_1x_1^n+c_2x_2^n+c_3x_3^n+c_4x_4^n时就有

a_{n+4}=6a_{n+3}-8a_{n+2}-6a_{n+1}+9a_n

然后，我们知道4个a_n的值，就可以把c_k解出来了。当然，这个过程是可以推导到任意次数的，用的时候可以直接记结论，就是求根，根据前几项解a_n=c_1x_1^n+c_2x_2^n+c_3x_3^n+c_4x_4^n，然后写出通项公式。

再多说几句，当一个特征多项式只能看出来几个因式的时候，或出现了不爽的根的时候，比如x^2+1，则可以把递推式化简成a_n+a_{n-2}=c_1(a_{n-1}+a_{n-3})+\cdots，有时可以便于求解。

好，看起来要结束了，我们先把方程解出来，如果都是有理根就很舒服了。我们移项得

x^4-6x^3+8x^2+6x-9=0

因式分解（可以用试根法）得

(x-3)^2(x-1)(x+1)=0

那么根就是x_1=x_2=3,x_3=1,x_4=-1。等一下，好像有点问题，这里出现了重根，那么解c_k时c_1和c_2前面的系数时完全相同的（因为x_1=x_2）。那怎么办呢？

这里再稍微拓展一下，讲一下x_1是它的k重根时的做法（当然这里只有二重根，但为了形象地描述，还是用这个方程）。设x^{n+4}=6x^{n+3}-8x^{n+2}-6x^{n+1}+9x^n移项后为f(x)=0，则(x-x_1)^k是它的因式，不妨设f(x)=(x-x_1)^kg(x)，显然f(x_1)=0，而由洛必达法，对于1\leqslant m<k，对m从小到大归纳，有

\begin{aligned}0&=(x_1-x_1)^{k-m}g(x_1) \\&=\lim_{x\to x_1}(x-x_1)^{k-m}g(x)\\&=\lim_{x\to x_1}\frac{(x-x_1)^kg(x)}{(x-x_1)^m}\\&=\lim_{x\to x_1}\frac{f(x)}{(x-x_1)^m}\\&=\lim_{x\to x_1}\frac{f'(x)}{((x-x_1)^m)'}\\&=\lim_{x\to x_1}\frac{f'(x)}{m(x-x_1)^{m-1}}\\&=\begin{cases}f'(x_1)&&m=1\\\lim\limits_{x\to x_1}\frac{f^{(2)}(x)}{m(m-1)(x-x_1)^{m-2}}=\cdots=\lim\limits_{x\to x_1}\frac{f^{(m)}(x)}{m!}=\frac{f^{(m)}(x_1)}{m!}&(\texttt{由归纳},f^{(m')}(x_1)=0(m'<m))&m>1\end{cases}\end{aligned} \therefore f^{(m)}(x_1)=0

而移项前得x^{n+4}=6x^{n+3}-8x^{n+2}-6x^{n+1}+9x^n，对它两边同时求m\ (1\leqslant m<k)次导应该也是成立的，所以有

(n+4)^\frac m{}x^{n+4-m}=6(n+3)^\frac m{} x^{n+3-m}-8(n+2)^\frac m{}x^{n+2-m}-6(n+1)^\frac m{}x^{n+1-m}+9n^\frac m{}x^{n-m}

稍微写好看一点，用n+m替代n：

(n+4+m)^\frac m{}x^{n+4}=6(n+3+m)^\frac m{} x^{n+3}-8(n+2+m)^\frac m{}x^{n+2}-6(n+1+m)^\frac m{}x^{n+1}+9(n+m)^\frac m{}x^n

发现每个x^n前得系数都是一个固定的m次的多项式，所以用0\leqslant m<k这里的m应该能表示所有的k-1次多项式，即

h(n+4)x^{n+4}=6h(n+3)x^{n+3}-8h(n+2)x^{n+2}-6h(n+1)x^{n+1}+9h(n)x^n

当h为小于k次的多项式时成立。

所以，和前面同理，这里待定的系数就是这个k-1次的多项式的系数。这又是一个可以直接记的结论：当有k重根x_1时，待定(c_{1,0}+c_{1,1}n+\cdots+c_{1,k-1}n^{k-1})x_1^n，其它项和原来的相同。

所以，在这道题中，不妨设a_i=(c_1n+c_2)3^n+c_3+c_4(-1)^n，根据a_0=-3,a_1=-6,a_2=-12,a_3=-6，解得c_1=\frac98,c_2=-\frac{117}{32},c_3=\frac98,c_4=-\frac{15}{32}

然而，这题需要卡常数。根据a^p\equiv a\pmod p，因为3\neq 0\pmod p，所以3^{p-1}\equiv1\pmod p，所以计算3^n时，可以计算3^{n\mod (p-1)}，这样，可以预处理3^{0,1,\cdots,65535}，和3^{0,65536,2\cdot65536,\cdots,65536\cdot65536}，然后 \Theta(1)算快速幂。当然还可以结合其他卡常技巧。

代码：

#include <bits/stdc++.h>
using namespace std;
typedef unsigned long long ull;
typedef long long ll;
typedef unsigned int uint;
ull sd;int op;
#define sz 65537
#define kcz 1000000007
ll t0[sz],t1[sz];
int main()
{
    int T;
    register int i,ans=0;
    ull n;
    scanf("%d%llu%d",&T,&sd,&op);
    for(i=1,t0[0]=1;i<sz;i++) t0[i]=t0[i-1]*3%kcz;
    for(i=1,t1[0]=1;i<sz;i++) t1[i]=t1[i-1]*t0[65536]%kcz; // 预处理
    for(;T--;)
    {
        sd ^= sd << 43,sd ^= sd >> 29,sd ^= sd << 34;
        if(op==2) n=sd;
        else if (op == 0) n = sd%USHRT_MAX + 1;
        else n = sd%UINT_MAX + 1; // 把Mker的内容放到主函数里
        ans^=(((ll)(n%kcz)*125000002+93749997)%kcz*t0[(n%(kcz-1))&65535]%kcz*t1[((n%(kcz-1))>>16)&65535]+((n&1)?343750004:906250007))%kcz; // 提前算9/8等，而9/8-15/32*(-1)^n是51/32或21/32（根据n的奇偶），也可以提前算
    }
    printf("%d\n",ans);
    return 0;
}