Project Euler 712 - Exponent Difference
Official link: https://projecteuler.net/problem=712
Note: Code takes around ~560 seconds to run, but I am very happy with my solution
Thought Process
Having seen many problems like this, I had a feeling the problem would have some trick with primes below and above the square root of N.
But first, I needed to have a better understanding of D(n,m) and S(N).
Lets exam S(10) more deeply
The prime less than 10 are: p = 2, 3, 5, 7
In the range [1, 10] we have:
floor(10/2^0) = 10 numbers divisible by 2^0 = 1
floor(10/2^1) = 5 numbers divisible by 2^1 = 2
floor(10/2^2) = 2 numbers divisible by 2^2 = 4
floor(10/2^3) = 1 numbers divisible by 2^3 = 8
Therefore it is clear that we have:
floor(10/2^0) - floor(10/2^1) = 5 numbers such that v2(n) = 0
They are: 1, 3, 5, 7, 9
floor(10/2^1) - floor(10/2^2) = 3 numbers such that v2(n) = 1
They are: 2, 6, 10
floor(10/2^2) - floor(10/2^3) = 1 numbers such that v2(n) = 2
They are: 4
floor(10/2^3) - floor(10/2^4) = 1 numbers such that v2(n) = 3
They are: 8
Using this information we can now understand S(N) much easier.
I continue with the p = 2 case:
Let's say n = 1, 3, 5, 7, 9. That is v2(n) = 0
Now if m = 1, 3, 5, 7, 9. Then v2(m) = 0 and D(n, m) = 0
Now if m = 2, 6, 10. Then v2(m) = 1 and D(n, m) = 1
Now if m = 4. Then v2(m) = 2 and D(n, m) = 2
Now if m = 8. Then v2(m) = 3 and D(n, m) = 3
Let's say n = 2, 6, 10. That is v2(n) = 1
Now if m = 2, 6, 10. Then v2(m) = 1 and D(n, m) = 0
Now if m = 4. Then v2(m) = 2 and D(n, m) = 1
Now if m = 8. Then v2(m) = 3 and D(n, m) = 2
Let's say n = 4. That is v2(n) = 2
Now if m = 4. Then v2(m) = 2 and D(n, m) = 1
Now if m = 8. Then v2(m) = 3 and D(n, m) = 2
Let's say n = 8. That is v2(n) = 3
Now if m = 8. Then v2(m) = 3 and D(n, m) = 0
All of these possibilities can be easily summed. For example if n = 1, 3, 5, 7, 9, m = 2, 6, 10, we have 5 choices for n, 3 choices for m such that D(n, m) = 1, so total contribution is 5*3*1 = 15.
We also don't forget than n, m are interchangeable so the total is multiplied by 2. Now we can do the exact same for the other primes and we now have a neat formula to calculate S(N)
Now, this method works great for up till 10^8, but we are limited because we need to generate primes up till N.
This is not feasible for N = 10^12, hence my intuition was correct and we need a seperate method for big primes.
By experimenting like above for p = 5, 7, I noticed something.
Let's say we have a prime p between N/2 and N, this means that floor(N/p) = 1, this means that floor(N/p^0) - floor(N/p^1) = N - 1 and floor(N/p^1) - floor(N/p^2) = 1.
Hence using the same method as above we find that the contribution of this prime is (N - 1)*1*1 (Since D(n, m) = 1) therefore we do not need to know what the primes are we just need to know how many primes there are between N/2 and N, which is where the Prime Counting Function comes in.
Knowing all this we can easily find how to calculate the contribution of all primes greater than the square root of N, and we use the original method to calculate contribution of primes less than the square root of N.
Now that the maths has been figured out the problem becomes how to calculate 𝜋(n) efficiently.
In my search I found the Meissel-Lehmer Algorithm which uses Legendre's Formula which was then adapted by Meissel and then one last time by Lehmer to get the algorithm I implemented.
The following article by Adithya Ganesh helped me a lot in understanding and implementing the algorithm efficiently.
This primepi function has been added to my mathslib python package for those interested to read the code.
Interactive Code
Enter an integer (yourinput)
Code will output S(yourinput)