Project Euler 492 Exploding sequence

Posted on | | 阅读数 次
Words count in article: 1.3k | Reading time ≈ 6

开启传送门

题目描述

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 \geq 1$.

Examples:

$a_3 $=$ 2359$

$a_6 $=$ 269221280981320216750489044576319$

$a_6 mod 1000000007 $=$ 203064689$

$a_{100} mod 1000000007 $=$ 456482974$

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

Examples:

$B(10^9, 10^3, 10^3)$ = $23674718882$

$B(10^9, 10^3, 10^{15})$ = $20731563854$

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

题意

给你一个二次递推序列,问你这个数列的第$10^{15}$项膜上一堆质数之后求和是多少.

分析

看了大佬的题解之后才做出来的,还是太菜了$\dots$

下面是分析.

首先我们令$b_n = 6a_n+5$

于是有

$b_1$ = $11$

$b_{n+1} $=$ b_n^2-2$

然后我们令　$ x + \frac{1}{x}$ = $11$

可以显然发现　$b_n$ = $x^{2^{n-1}} + \frac{1}{x^{2^{n-1}}}$

于是不妨设　$x$ = $\frac{11+\sqrt{117}}{2}$

$\therefore b_n$ = $(\frac{11+\sqrt{117}}{2})^{2^{n-1}} + (\frac{11-\sqrt{117}}{2})^{2^{n-1}}$

然后我们可以直接算就行了

我们搞出递推式，然后用矩阵乘法的做法来做

我们考虑一个更一般的问题，

令$f_n = (\frac{11+\sqrt{117}}{2})^n + (\frac{11-\sqrt{117}}{2})^n$

然后我们考虑$f_n$如何求解

首先我们假设$f_n = x\ast f_{n-1}+y\ast f_{n-2}$

$$令 (\frac{11+\sqrt{117}}{2})^{n-2} = a$$

$$令 (\frac{11-\sqrt{117}}{2})^{n-2} = b$$

则有
\begin{cases}
\begin{eqnarray}
f_n &=& a\ast (\frac{11+\sqrt{117}}{2})^2+b\ast (\frac{11-\sqrt{117}}{2})^2\\
&=& a\ast (\frac{238+22\sqrt{117}}{4})+b\ast (\frac{238-22\sqrt{117}}{4})\\
f_{n-1}&=&a\ast (\frac{11+\sqrt{117}}{2})+b\ast (\frac{11-\sqrt{117}}{2})\\
f_{n-2}&=&a+b \\
\end{eqnarray}
\end{cases}

\begin{eqnarray}
&\therefore& a\ast (\frac{238+22\sqrt{117}}{4})+b\ast (\frac{238-22\sqrt{117}}{4})\\
&=& x\ast a\ast (\frac{11+\sqrt{117}}{2})+x\ast b\ast (\frac{11-\sqrt{117}}{2}) + y\ast (a+b)\\
\end{eqnarray}

将非根号项提出来，可以得到

$\therefore a\ast \frac{119}{2}+b\ast \frac{119}{2}$=$xa\ast \frac{11}{2}+xb\ast \frac{11}{2}+y(a+b) $

同理将根号项提出可以得到

$a\ast \frac{11}{2}+b\ast \frac{-11}{2}$=$xa\ast \frac{1}{2}+xb\ast \frac{-1}{2} $

然后可以解得

\begin{cases}
x = 11\\
y = -1
\end{cases}

也就是说$f_n$ = $11\ast f_{n-1} - f_{n-2}$

好了，我们得到了$f_n$的递推式，然后考虑原问题，也就是说

$$
\left(
\begin{matrix}
b_{n+1}\\
b_n
\end{matrix}
\right)
=\left(
\begin{matrix}
11 & -1\\
1 & 0
\end{matrix}
\right)^{2^{n-1}-1}
\left(
\begin{matrix}
b_2\\
b_1
\end{matrix}
\right)
=\left(
\begin{matrix}
11 & -1\\
1 & 0
\end{matrix}
\right)^{2^{n-1}-1}
\left(
\begin{matrix}
119\\
11
\end{matrix}
\right)
$$

$\because n$ 很大($10^{15}$)

$\therefore 2^{n-1}-1$很大,以至于我们不能方便的计算

然后我们考虑降低指数

具体分析看大佬的分析,可以得到循环节可以是$(p+1)(p-1)$

这个值我们可以接受,然后直接裸的矩阵ksm就行了$\dots$

代码

给大佬递上我奇丑无比的代码　(/ω＼)

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
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
#include<bits/stdc++.h>
using namespace std;
const long long n = 1e15;
const long long l = 1e9;
const long long r = l+1e7;
long long add(long long a,long long b,long long mod){ return (a+b)%mod; }
long long mul(long long a,long long b,long long mod){
	a = add(a,mod,mod);
	b = add(b,mod,mod);
	long long ret = 0;
	while(b){
		if(b&1) ret = add(ret,a,mod);
		a = add(a,a,mod);
		b>>=1;
	}
	return ret;
}
long long qpow(long long a,long long b,long long mod){
	long long ret = 1;
	while(b){
		if(b&1) ret = mul(ret,a,mod);
		a = mul(a,a,mod);
		b>>=1;
	}
	return ret;
}
bool is_prime(int a){
	if(a&1){
		for(int i = 3;i*i<=a;i++) if(a%i==0) return false;
		return true;
	}
	return false;
}
struct p{
	int mat[2][2];
	void show(){
		for(int i = 0;i<2;i++){
			for(int j = 0;j<2;j++) cout<<mat[i][j]<<' ';
			cout<<endl;
		}
	}
};
p mul(p a,p b,long long mod){
	p ret;
	memset(ret.mat,0,sizeof ret.mat);
	for(int i = 0;i<2;i++){
		for(int j  =0;j<2;j++){
			for(int k = 0;k<2;k++){
				ret.mat[i][j] = add(ret.mat[i][j],mul(a.mat[i][k],b.mat[k][j],mod),mod);
			}
		}
	}
	return ret;
}
p qpow(p a,long long b,long long mod){
	p ret;
	for(int i = 0;i<2;i++) for(int j = 0;j<2;j++) {
		if(i==j) ret.mat[i][j]=1;
		else ret.mat[i][j]=0;
	}
	while(b){
		if(b&1) ret = mul(ret,a,mod);
		a = mul(a,a,mod);
		b>>=1;
	}
	return ret;
}
long long solve(long long prime){
	const long long MOD = prime*prime-1;

	long long zhi = qpow(2,n-1,MOD);
	zhi = add(zhi,MOD-1,MOD);

	p base;
	base.mat[0][0]=11,base.mat[0][1]=-1;
	base.mat[1][0]=1 ,base.mat[1][1]=0;
	
	base = qpow(base,zhi,prime);
	long long ret1 = mul(base.mat[1][0],119,prime);
	long long ret2 = mul(base.mat[1][1],11,prime);
	return add(ret1,ret2,prime);
}
long long inv(long long a,long long mod){
	if(a==1) return 1;
	return mul(mod-mod/a,inv(mod%a,mod),mod);
}
int main(){
	long long ans = 0;
	for(int i = l;i<=r;i++){
		if(is_prime(i)){
			long long bn = solve(i);
			long long now1 = add(bn,i-5,i);
			long long an = mul(now1,inv(6,i),i);
			ans+=an;
		}
		if(i%100000==0) cout<<i<<endl;
	}
	cout<<ans<<endl;
	return 0;
}