Here’s kind of an unusual question: What is the density of integers that have an even number of prime factors with an exponent greater than 1?
To define the density, you take the proportion up to an integer N then take the limit as N goes to infinity.
It’s not obvious that the limit should exist. But if it does exist, you might guess that it’s 1/2; it’s a question involving whether something is even or odd, and so even and odd occurrences might balance each other out.
But the result, known as the Feller-Tonier constant, is bigger than 1/2. This seems more reasonable when you consider that zero is an even number and a lot of numbers may not have any prime factors with exponent bigger than 1.
Here’s a little Python code to compute the ratio for N = 106.
from sympy.ntheory import factorint def f(n): v = factorint(n).values() n = len([x for x in v if x > 1]) return n % 2 == 0 c = 0 N = 1000000 for n in range(N): if f(n): c += 1 print (c/N)
This gives a ratio of 0.661344. The limiting value is 0.661370494….
Computing the Feller-Tornier constant
The code above gives an estimate for the Feller-Tornier constant, but how accurate is it? Is there a more efficient way to estimate it that doesn’t require factoring a lot of numbers?
The constant is given analytically by
where pn is the nth prime.
Here’s some Python code to estimate the Feller-Tornier constant using the expression above.
N = 1000000 product = 1 for p in primerange(2, N): product *= 1 - 2*p**-2 print(0.5*(1+product))
Note that we used
primerange rather than
prime(n). If you need to generate a range of prime numbers, the former is much more efficient.
The code gives us an estimate of 0.66131707, which is good to seven decimal places. Furthermore, we know this is an upper bound on the exact value because we’ve left out terms in the infinite product that are less than 1.
One could calculate a lower bound as well and thus know an interval that must contain the true value. The lower bound is left as an exercise to the reader.