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.

Project Euler Problem 779

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_term stores

$$\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