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} &a_3 = 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);
}