Square Free Integers

An integer $n$ is said to be square-free, if it is not divisible by any perfect square other than $1$. That is, each prime that appears in its prime factorisation is raised to the power of one.

So, for example, $6$ is square free but $27$ is not since it is divisble by $3^2$.

Let's consider, the probabilties of choosing integers randomly. In general, the probability of a number being divisible by $n$ is $\frac{1}{n}$. This is because a multiple of $n$ comes every $n$ numbers. Of course, in this we are assuming that each number is equally likely to be chosen, which is reasonable.

$$\therefore\hspace{5pt}P\left(\hspace{2pt}p^2\text{ divides }x\hspace{2pt}\right)=\frac{1}{p^2}\hspace{10pt}\text{for }x\in\mathbb{N}$$

$$\implies P\left(\hspace{2pt}\text{all primes }p^2\text{ do not divide }x\hspace{2pt}\right)=\left(1-\frac{1}{p_1^2}\right)\cdot\left(1-\frac{1}{p_2^2}\right)\cdot\left(1-\frac{1}{p_3^2}\right)\dots$$

Now, consider the famous standard result A.K.A. the Basel problem:

$$\sum_{i=1}^{\infty}\frac{1}{i^2}=1+\frac{1}{4}+\frac{1}{9}+\dots=\frac{\pi^2}{6}$$

Let $S=\sum_{i=1}^{\infty}\frac{1}{i^2}=\frac{1}{1^2}+\frac{1}{2^2}+\frac{1}{3^2}+\dots$ Since we are trying to involve $\left(1-\frac{1}{{p_k}^2}\right)$ into our calculations, we will do the following:

$$\frac{1}{2^2}S=\frac{1}{2^2}+\frac{1}{4^2}+\frac{1}{6^2}+\dots$$

$$\implies S-\frac{1}{2^2}S=\left(1-\frac{1}{2^2}\right)S=\frac{1}{1^2}+\frac{1}{3^2}+\frac{1}{5^2}+\dots$$

Notice how we have now removed all fractions with multiples of $2$ in their denominator from the RHS of the equation above. We will continue in this way.

$$\frac{1}{3^2}\left(1-\frac{1}{2^2}\right)S=\frac{1}{3^2}+\frac{1}{9^2}+\frac{1}{15^2}+\dots$$

$$\implies\left(1-\frac{1}{2^2}\right)S-\frac{1}{3^2}\left(1-\frac{1}{2^2}\right)S=\left(1-\frac{1}{3^2}\right)\left(1-\frac{1}{2^2}\right)S=\frac{1}{1^2}+\frac{1}{5^2}+\frac{1}{7^2}+\frac{1}{11^2}+\dots$$

So all fractions which have denominators with factors of 2 and 3 are now removed from the RHS. We can carry on for all primes, removing all the numbers from the RHS, until we are only left with 1.

$$\therefore\hspace{5pt}\dots\left(1-\frac{1}{7^2}\right)\left(1-\frac{1}{5^2}\right)\left(1-\frac{1}{3^2}\right)\left(1-\frac{1}{2^2}\right)S=1$$

$$\implies \prod_{k=1}^{\infty}\left(1-\frac{1}{{p_k}^2}\right)=\frac{1}{S}=\boxed{\frac{6}{\pi^2}}\hspace{10pt}\square$$

Thus, we have demonstrated that $P($all primes $p^2$ do not divide $x)=\frac{6}{\pi^2}$ i.e. the probability of a positive integer being square free is $\frac{6}{\pi^2}$ !

$$\frac{6}{\pi^2}\approx0.60792710185\dots$$

We can confirm the validity of this elegant result iteratively through a Python algorithm. The following code defines a function isSquareFree to determine whether a given number is square-free or not.Note that the for loop above the function definition is used for the iteration later on.


          b=np.array([])
          for k in range(0,100):
              def isSquareFree(n):
                  if n % 2 == 0:
                      n = n / 2
                  if n % 2 == 0:
                      return False
                  for i in range(3, int(np.sqrt(n) + 1)):
                      if n % i == 0:
                          n = n / i
                      if n % i == 0:
                          return False
                  return True

We then define an array $a$, which contains 10000 ranomdly generated integers between 0 and 1000. From this, we determine the proportion of numbers which are square-free and store this in $b$, previously defined.


              a=np.random.randint(1000,size=10000)

              count=0
              for i in range(0, len(a)):
                  if isSquareFree(a[i]) == True:
                      count=count+1
              x = count/len(a)
              b = np.append(b,x)

As this whole process is repeated 100 times (see the second line of code in the first code block above), the array $b$ contains 100 elements of proportions from 100 different sets of 10000 numbers.

We can output this data in various formats; for instance, we can print the array $b$ and work out its mean.

$$b = [0.6024, 0.6177, 0.6091, 0.6123, 0.607, 0.6146, 0.6031, 0.6124, 0.614, 0.6136, 0.6034, 0.6096, 0.6047, 0.6131, 0.6048, 0.614, 0.617, 0.607, 0.6115, 0.6096, 0.6054, 0.6033, 0.6104, 0.6067, 0.6133, 0.6142, 0.607, 0.6069, 0.6046, 0.6032, 0.6086, 0.604, 0.6079, 0.6172, 0.6042, 0.6127, 0.6096, 0.5994, 0.6076, 0.6008, 0.613, 0.6108, 0.6082, 0.6055, 0.6087, 0.6124, 0.6024, 0.6045, 0.6021, 0.6065, 0.6059, 0.6111, 0.6024, 0.607, 0.6113, 0.6125, 0.5997, 0.6122, 0.6122, 0.6063, 0.6083, 0.6124, 0.6122, 0.6094, 0.6145, 0.6075, 0.607, 0.613, 0.6065, 0.6018, 0.604, 0.604, 0.6045, 0.601, 0.6033, 0.6096, 0.6037, 0.6111, 0.6112, 0.61, 0.6093, 0.6074, 0.6073, 0.5981, 0.6186, 0.6059, 0.6097, 0.6087, 0.6065, 0.6044, 0.5977, 0.6141, 0.6012, 0.6095, 0.6007, 0.608, 0.6066, 0.5999, 0.6132, 0.609]$$

$$\therefore\hspace{5pt}\text{mean}(b)=\underline{0.607934}$$

Observe how extremely close the mean value computed above is to the true value of $\frac{6}{\pi^2}$. We can also graph all the iterations:

── KN7811 ──