Project Euler 492

题目

Exploding sequence

Define the sequence $a_1, a_2, a_3, \dots$ as:

• $a_1 = 1$

• $a_{n+1} = 6a_n^2 + 10a_n + 3$ for $n \ge 1$.

Examples:

$\begin{aligned}
&a3 = 2359\
&a_6 = 269221280981320216750489044576319\
&a_6 \bmod 1 000 000 007 = 203064689\
&a{100} \bmod 1 000 000 007 = 456482974
\end{aligned}$

Define $B(x,y,n)$ as $\sum (a_n \bmod p)$ for every prime $p$ such that $x \le p \le x+y$.

Examples:

$\begin{aligned}
&B(10^9, 10^3, 10^3) = 23674718882\
&B(10^9, 10^3, 10^{15}) = 20731563854
\end{aligned}$

Find $B(10^9, 10^7, 10^{15})$.

解决方案

构造一个新序列$b_n=6a_n+5$，并且回代到$a_n$的递推式中，那么得到$b_n$的递推式：

$$\left\{\begin{aligned} &b_1=11\\ &b_{n+1}=b_n^2-2 \end{aligned}\right.$$

可以发现$b_n$的式子比起$a_n$消去了一次项，并且$b_n^2$的项系数恰好为$1$。

根据$b_n$的递推式，再联想到如下关于$x$的等式：

$$x^2+\dfrac{1}{x^2}=\left(x+\dfrac{1}{x}\right)^2-2$$

如果将$b_1$写成形如$b_1=k+\dfrac{1}{k}$的形式，那么$b_n$的通项公式可以很轻易地写出来：

$$b_n=k^{2^{n-1}}+\dfrac{1}{k^{2^{n-1}}}$$

求解$k+\dfrac{1}{k}=11$，得到$k=\dfrac{11\pm 3\sqrt{13}}{2}$。

接下来参考两种情况：

• 如果$13$是$p$的二次剩余，那么计算出$13$的二次剩余后，直接代入计算$k$的值即可。
• 如果$13$不是$p$的二次剩余，那么接下来这个地方参考了Thread的一些内容：考虑在域$\mathbb{Z}_p[\sqrt{13}]$进行运算，这个域的阶为$p^2-1$。这个域的乘法运算性质为：$(x_1,y_1)\times(x_2,y_2)=(x_1x_2+13y_1y_2,x_1y_2+x_2y_1)$，相当于是$x_1+y_1\sqrt{13}$和$x_2+y_2\sqrt{13}$这两个数相乘。那么根据这个运算性质，$\left(\dfrac{11}{2},\dfrac{3}{2}\right)^{2^{n-1}}$很容易就能够计算出来。$\left(\dfrac{11}{2},\dfrac{3}{2}\right)$在$\mathbb{Z}_p[\sqrt{13}]$上的逆元是$\left(\dfrac{11}{2},-\dfrac{3}{2}\right)$，具有共轭的性质。最终$\left(\dfrac{11}{2},\dfrac{3}{2}\right)^{2^{n-1}}+\left(\dfrac{11}{2},-\dfrac{3}{2}\right)^{2^{n-1}}$的$y$值一定为$0$，求出的$x$值即为最终的$b_n$。

计算出$b_n$后，回代得到$a_n=\dfrac{b_n-5}{6}$，直接计算即可。

代码

# include <bits/stdc++.h>
# include "tools.h"
using namespace std;
typedef long long ll;


const int L=1000000000;
const int M=10000000;
const ll N=1000000000000000;

struct Zpsqrt{
    static const ll D=13;
    ll x,y,p;
    Zpsqrt operator * (Zpsqrt a){
        return Zpsqrt{(x*a.x+y*a.y%p*D)%p,(x*a.y+y*a.x)%p,p};
    }
    Zpsqrt qpow(ll m){
        Zpsqrt a{1,0,p},n=*this;
        for(;m;m>>=1){
            if(m&1) a=a*n;
            n=n*n;
        }
        return a;
    }
};
ll cal(ll n,ll p){
    ll s;
    ll inv2=(p+1)>>1;
    if(qpow(13,(p-1)>>1,p)==1){
        ll e=qpow(2,n-1,p-1);
        ll u=3ll*sqrt_mod_prime(13,p)+11;
        ll k=u*inv2%p;
        ll x=qpow(k,e,p);
        ll y=mod_inverse(x,p);
        s=x+y;
    }
    else{
        Zpsqrt a{inv2*11%p,inv2*3%p,p};
        ll e=qpow_mpz(2,n-1,p*p-1);
        a=a.qpow(e);
        s=a.x<<1;
    }
    s%=p;
    return (s-5+p)*mod_inverse(6,p)%p;
}
int main(){
    ll ans=0;
    for(int p=next_prime(L-1);p<=L+M;p=next_prime(p)){
        ans+=cal(N,p);
    }
    printf("%lld\n",ans);
}