A numerical curiosity about the Kempner series I am writing some notes for a course in basic Arithmetics and Calculus: today I was sketching a brief paragraph about Benford's law and the following well-known fact, related to Cauchy condensation test:

Let $E_7$ be the set of positive natural numbers such that no $7$
   occurs in their decimal representation. Then $\mathcal{S}_7 =
 \sum_{n\in E_7}\frac{1}{n}$ is convergent to a real number $<89$.

Now my simple question: 

Q: What is an efficient method for finding an accurate approximation of the exact value of $\mathcal{S}_7$?

A proof of convergence is simply derived from
$$ \sum_{\substack{n\in E_7\\ 10^m \leq n < 10^{m+1}}}\!\!\!\!\!\!\!\frac{1}{n}\leq \frac{9^{m+1}}{9\cdot 10^m} $$
but such inequality only produces very crude bounds. The Schnirelmann density is very likely to be involved in a more "professional" approach.
 A: Here's a summary of Baillie's paper I found in the comments, since not everyone has JSTOR access:
Fix $m$; in your case, $m=7$.  Let $S_i$ be the set of numbers with no digit $m$ with $i$ digits.  Then,
$$S_{i+1}=\{10x+k|x \in S_i, k \in {0,...,9}, k \ne m\} $$
To compute the sum you desire, we define
$$s(i,j):=\sum_{x \in S_i}x^{-j}, $$
then, the Kempner sum is
$$\sum_{i=1}^{\infty}s(i,1). $$
But,
$$s(i+1,j)=\sum_{x \in S_{i-1}} \sum_{\substack{k=0\\k \ne m}}^9 (10x+k)^{-j} $$
and the inner term $(10x+k)^{-j}$ can be expanded to
$$(10x+k)^{-j}=(10x)^{-j} \sum_{n=0}^\infty \binom{n+j-1}{n}k^n(10x)^{-n} $$
so we can rearrange terms to get 
$$s(i+1,j)=\sum_{n=0}^\infty a(j,n) s(i,j+n) $$
with
$$a(j,n)=10^{-(j+n)}\sum_{\substack{k=0\\k \ne m}}^9 k^n \binom{n+j-1}{n} $$
For $j>1$ we can bound $s(i,j)$ by the corresponding sum over all integers $\ge 10^i$ which converges exponentially fast so if you want to calculate some $M$ digits:


*

*for $i$ sufficiently large we can approximate $s(i+1,j) \approx
   a(j,0) s(i,j)=(9/10)s(i,j)$.

*For $i$ fixed, $s(i,j)$ will also
decrease exponentially

*Baillie in 1979 calculated $s(i,1)$ explicitly for $i \le
   4$ but the Python one-liner sum(1/x for x in range(1,1000000) if '7' not in str(x)) computes $\sum_{i=1}^6 s(i,1)$ and runs in roughly one second.

