The core function is described here, however I am unhappy with this as in the description of the OEIS it specifically links to Project Euler Problem 795 which kind of gives away that I should stick with this formula (Although there were other solutions).

The implementation of this formula is extremely simple if you have the prime factorization of all the d's which divide n, the only tricky part is to make that part efficient which I did by using a smallest prime factor sieve and quickly building a table with the prime factorization of every number from 1 to n.

After that we just apply the formula (Lots of simplification can be made if you look carefully). Now that the problem is done, with a sour taste left in my mouth I moved on to the much more challenging and fun part for me, proving why this identity is true!

Under the Wikipedia GCD Properties section there is an interesting one which we will be using

Now let's continue

Most of the hard work is done now, we just look (as the original identity suggests) at the 2 cases where n is even or odd. I start with the easier one first, n is odd case

And lastly the meatiest part, the case where n is even