Project Euler 368 A Kempner-like series

开启传送门

题目描述

The harmonic series $1+\frac{1}{2}+\frac{1}{3}+\frac{1}{4}+\dots$ is well known to be divergent.

If we however omit from this series every term where the denominator has a $9$ in it, the series remarkably enough converges to approximately $22.9206766193$.

This modified harmonic series is called the Kempner series.

Let us now consider another modified harmonic series by omitting from the harmonic series every term where the denominator has $3$ or more equal consecutive digits. One can verify that out of the first $1200$ terms of the harmonic series, only $20$ terms will be omitted.
These $20$ omitted terms are:

$$\frac{1}{111},\frac{1}{222},\frac{1}{333},\frac{1}{444},\frac{1}{555},\frac{1}{666},\frac{1}{777},\frac{1}{888},\frac{1}{999},\frac{1}{1000},\frac{1}{1110}, \frac{1}{1111},\frac{1}{1112},\frac{1}{1113},\frac{1}{1114},\frac{1}{1115},\frac{1}{1116},\frac{1}{1117},\frac{1}{1118},\frac{1}{1119}$$

This series converges as well.

Find the value the series converges to.
Give your answer rounded to $10$ digits behind the decimal point.

题意

就是,调和级数不是发散的嘛,然后让你删掉那些分母含有连续三个相同数字的,然后可以证明剩下的级数是收敛的,然后问你收敛于多少.

分析

做法比较巧妙,下面一点点分析.

首先我们令$S_1(n,d)$表示一个包含所有这样的$n$位数的集合,集合中所有的$n$位数最后一位是$n$,并且倒数第二位不是$n$.

$S_2(n,d)$表示一个包含所有这样的$n$位数的集合,集合中所有的$n$位数最后两位是$n$,并且倒数第三位不是$n$.

然后我们令

$f_1(n,d,j)$ = $\sum_{x\in S_1(n,d)} \frac{1}{x^j}$

$f_2(n,d,j)$ = $\sum_{x\in S_2(n,d)} \frac{1}{x^j}$

然后答案显然就是

$\sum_{i=1}^{99} \frac{1}{i}+\sum_{n=3}^{\infty } \sum_{d=0}^{9} [f_1(n,d,1)+f_2(n,d,1)]$

关键是我们怎么算$f_1(n,d,j)$和$f_2(n,d,j)$

二者分析方法一样,这里给出$f_2(n,d,j)$怎么推出来的.

首先根据定义我们有

$f_2(n,d,j)$ = $\sum_{x\in S_2(n,d)} \frac{1}{x^j}$

我们不妨枚举每一个$x$,然后累加就是答案,所以子问题就是如何快速计算$\frac{1}{x^j} where x \in S_2(n,d)$

因为最后一位是$d$,所以我们不妨设$x$ = $y\ast 100+d\ast 10+d$

然后可以得到

$$
(\frac{1}{x\ast 100+d\ast 10+d})^j=\frac{1}{10^j} \cdot \frac{1}{(x\ast 10+d)^j} \cdot (\frac{1}{1+\frac{d}{x\ast 100+d\ast 10}})^j
$$

然后我们将第三个分式化简.

考虑到

$(1+a)^{-b} = \sum_{i=0}^{\infty } C_{b+i-1}^i \cdot (-a)^b$

所以上式可以化简为

\begin{eqnarray}
(\frac{1}{x\ast 100+d\ast 10+d})^j
&=& \frac{1}{10^j} \cdot \frac{1}{(x\ast 10+d)^j} \cdot \sum_{i=0}^{\infty } C_{j+i-1}^{i}(\frac{-d}{x\ast 100+d\ast 10})^i\\
&=& \frac{1}{10^j} \cdot \frac{1}{(x\ast 10+d)^j} \cdot \sum_{i=0}^{\infty } C_{j+i-1}^{i}(\frac{-d}{10})^i\cdot(\frac{1}{x\ast 10+d})^i\\
&=& \frac{1}{10^j} \sum_{i=0}^{\infty } C_{j+i-1}^{i} \cdot (\frac{-d}{10})^i \cdot (\frac{1}{x\ast 10+d})^{i+j}\\
\end{eqnarray}

然后我们考虑累加所有的$x$也就是原来的$f_2(n,d,j)$

可以得到

$$f_2(n,d,j) = \frac{1}{10^j} \sum_{i=0}^{\infty } (\frac{-d}{10})^i \cdot C_{j+i-1}^{i} \cdot f_1(n-1,d,i+j)$$

同理可以推得 $f_1(n,d,j)$

然后就是代码如下啦啦啦$\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
#include<bits/stdc++.h>
using namespace std;
const int maxm = 40;
double C[maxm*3][maxm*3];
double f1[2][10][maxm+1];
double f2[2][10][maxm+1];
double ans = 0;
void init(){//预处理出组合数和前99项的答案
for(int i = 0;i<maxm*3;i++) for(int j = 0;j<=i;j++){
if(j==0) C[i][j] = 1;
else C[i][j] = C[i-1][j-1]+C[i-1][j];
}
for(int i = 1;i<100;i++) ans+=1.0/i;
}
int main(){
init();
for(int d = 0;d<=9;d++){//预处理出n = 3的情况
for(int j = 1;j<=maxm;j++){
for(int pre = 10;pre<100;pre++){
int fir = pre/10;
int sec = pre%10;
if(fir==sec&&sec==d) continue;
if(sec==d) f2[1][d][j]+=1.0/pow(pre*10+d,j);
else f1[1][d][j]+=1.0/pow(pre*10+d,j);
}
}
ans+=f1[1][d][1]+f2[1][d][1];
}
for(int nouse = 4;nouse<=10000;nouse++){//滚动的算n>=4的时候
int now = nouse&1;
int las = now^1;
for(int d = 0;d<=9;d++){
for(int j = 1;j<=maxm;j++){
f1[now][d][j] = f2[now][d][j] = 0;//别忘了初始化
double pre = 1.0/pow(10.0,j);
for(int i = 0;i<=maxm;i++){
if(i+j<=maxm) f2[now][d][j]+=pre*C[j+i-1][i]*f1[las][d][i+j];
for(int x = 0;x<10;x++){
if(x==d) continue;
if(i+j<=maxm) f1[now][d][j]+=pre*C[j+i-1][i]*(f1[las][x][i+j]+f2[las][x][i+j]);
}
pre*=-d/10.0;
}
}
ans+=f1[now][d][1]+f2[now][d][1];
}
printf("%.10f\n",ans);
}
return 0;
}