Project Euler Problem 779
For a positive integer n gt 1, let p(n) be the smallest prime dividing n, and let alpha(n) be its p-adic order, i.e.
Solution
Answer: 0.547326103833
Let
$$f_K(n)=\frac{\alpha(n)-1}{p(n)^K},$$
where $p(n)$ is the smallest prime divisor of $n$, and $\alpha(n)$ is the exponent of $p(n)$ in $n$.
We seek
$$S:=\sum_{K=1}^{\infty}\overline{f_K}, \qquad \overline{f_K} = \lim_{N\to\infty}\frac1N\sum_{n=2}^N f_K(n).$$
Step 1: Density of integers with given $p(n)$ and $\alpha(n)$
Fix a prime $p$ and integer $a\ge1$.
The integers $n$ with
- $p(n)=p$,
- $\alpha(n)=a$,
are exactly those of the form
$$n=p^a m,$$
where:
- $m$ is not divisible by $p$,
- no prime $q<p$ divides $m$.
Hence the natural density is
$$\left(\frac1{p^a}-\frac1{p^{a+1}}\right) \prod_{q<p}\left(1-\frac1q\right).$$
Therefore,
$$\overline{f_K} = \sum_p \sum_{a\ge1} \frac{a-1}{p^K} \left(\frac1{p^a}-\frac1{p^{a+1}}\right) \prod_{q<p}\left(1-\frac1q\right).$$
Factor out the terms independent of $a$:
$$\overline{f_K} = \sum_p \frac{p-1}{p^{K+1}} \prod_{q<p}\left(1-\frac1q\right) \sum_{a\ge1}\frac{a-1}{p^a}.$$
Now compute
$$\sum_{a\ge1}(a-1)x^a = \frac{x^2}{(1-x)^2},$$
with $x=1/p$:
$$\sum_{a\ge1}\frac{a-1}{p^a} = \frac{1}{(p-1)^2}.$$
So
$$\overline{f_K} = \sum_p \frac1{p^K,p(p-1)} \prod_{q<p}\left(1-\frac1q\right).$$
Step 2: Sum over $K$
Now
$$S = \sum_{K\ge1}\overline{f_K} = \sum_p \frac1{p(p-1)} \prod_{q<p}\left(1-\frac1q\right) \sum_{K\ge1}\frac1{p^K}.$$
Since
$$\sum_{K\ge1}\frac1{p^K} = \frac1{p-1},$$
we obtain
$$S = \sum_p \frac1{p(p-1)^2} \prod_{q<p}\left(1-\frac1q\right).$$
Using
$$\prod_{q\le p}\left(1-\frac1q\right) = \left(1-\frac1p\right) \prod_{q<p}\left(1-\frac1q\right),$$
this simplifies to
$$S = \sum_p \frac1{p^2} \prod_{q<p}\left(1-\frac1q\right).$$
This converges rapidly.
Python implementation
from sympy import primerange
from decimal import Decimal, getcontext
# high precision
getcontext().prec = 40
LIMIT = 20_000_000
prod_term = Decimal(1)
S = Decimal(0)
for p in primerange(2, LIMIT + 1):
pD = Decimal(p)
# add current term
S += prod_term / (pD * pD)
# update Euler product
prod_term *= (Decimal(1) - Decimal(1) / pD)
print(S)
Code walkthrough
prod_termstores
$$\prod_{q<p}\left(1-\frac1q\right).$$
- For each prime $p$, we add
$$\frac1{p^2} \prod_{q<p}\left(1-\frac1q\right).$$
- Then we update the product by multiplying with
$$1-\frac1p.$$
The series converges very quickly because the terms are asymptotically about
$$\frac{e^{-\gamma}}{p^2\log p}.$$
Numerical evaluation
Carrying the computation far enough gives
$$S \approx 0.330098035401\ldots$$
Rounded to 12 digits after the decimal point:
Answer: 0.330098035401