Computing a Number Theoretic Sum

Published Jan 4, 2025

While I was attempting this problem, I misread it after writing it down on paper and instead solved

S = \sum_{xy \leq n} xy \gcd(x, y)

Here the sum is over positive integers $(x, y)$ satisfying $xy \leq n$ , rather than $x, y \leq n$ . In my opinion it turned out more interesting than the original one, so I wanted to share.

Compute $S = \sum_{xy \leq n} xy \gcd(x, y)$ in sublinear time. Here I present $\mathcal{O}(n^{2/3})$ but it can be much faster.

—Source: Accidental misread

Solution

Setting up notation first, let $M$ be set of multiplicative functions, and * be the Dirichlet product, then $(M, *)$ is an abelian group. Let $\epsilon$ , $1$ , $I$ , $\tau$ , $\sigma$ , $\mu$ , $\phi$ be the $*$ -identity ( $1$ -indicator), constant $1$ , identity, divisor count, divisor sum, Mobius, and totient functions respectively. Note that they are all multiplicative.

Write $k = xy$ so that we can iterate over $k$ as follows:

S = \sum_{k \leq n} k \sum_{d \mid k} \gcd(d, k/d)

Key Lemma $\sum_{d^2 \mid k} \tau(k/d^2) \phi(d) = \sum_{d \mid k} \gcd(d, k/d)$
Note here that $d$ does not represent the same iterator, but I kept the variable name anyways.

Proof of Lemma. A cool fact is that any sum of a product of multiplicative functions with the form of the left side is multiplicative. We can define here

f(d^2) = \begin{cases} 0 & \text{$d$ is not a square} \\ \phi(d) & \text{otherwise} \end{cases}

Note that $f$ is multiplicative, and the left side is actually $f * \tau$ , making it multiplicative.

We can verify the right side is multiplicative too, with just the textbook approach. One can notice that if $\gcd(p, q) = 1$ , $d \mid pq \implies d = rs$ , where $r \mid p$ and $s \mid q$ . Also, $\gcd(rs, pq/(rs)) = \gcd(r, p/r) \gcd(s, q/s)$ . Therefore,

\sum_{d \mid pq} \gcd(d, pq/d) = \sum_{r \mid p} \gcd(r, p/r) \sum_{s \mid q} \gcd(s, q/s)

It remains to simply check both sides coincide at prime powers:

\sum_{i \leq \lfloor a/2\rfloor} \tau(p^{a - 2i}) \phi(p^i) = \tau(p^a) + \sum_{1 \leq i \leq \lfloor a/2\rfloor} (a - 2i + 1) p^{i - 1}(p - 1)

= \sum_{i \leq \lfloor a/2\rfloor - 1} 2p^i + (a - 2\lfloor a/2\rfloor + 1) p^{\lfloor a/2\rfloor} = \sum_{i \leq a} p^{\min(i, a - i)}

So we’re done. The sum telescopes in the middle. $\square$

Therefore

S = \sum_{k \leq n} k \sum_{d^2 \mid k} \tau(k/d^2) \phi(d)

The goal here was to turn the sum into one consisting of multiplicative functions, so we can try applying Dirichlet’s Hyperbola Method. We can quickly make the observation $d^2 \leq k \leq n$ , and it looks promising to try iterating over $d \leq \sqrt{n}$ instead:

S = \sum_{d \leq \sqrt{n}} \phi(d) \sum_{d^2 \mid k} k \tau(k/d^2)

Now $k$ has a really simple characterization: it is just $d^2 \alpha$ for $\alpha \leq n/d^2$ . Let’s iterate over $\alpha$ ,

S = \sum_{d \leq \sqrt{n}} \phi(d) d^2 \sum_{\alpha \leq n/d^2} \alpha \tau(\alpha)

So we would be happy if we can compute prefix sums of $I \tau$ quickly, luckily it is multiplicative and simple enough, so let’s try to write it as a Dirichlet product:

I\tau = I * I = I * \mu * \sigma = \phi * \sigma

At this point we are almost at the end, as many methods to calculate prefix sums of $\phi$ and $\sigma$ are well known:

$\sigma = I * 1$ , which decomposes it into two functions whose prefix sums can be computed in $\mathcal{O}(1)$ . With the Hyperbola Method we can compute prefix sums of $\sigma$ in $\sqrt{n}$ .

Prefix sums of $\phi$ can be computed in $\mathcal{O}(n^{2/3})$ by linear sieving the first $n^{2/3}$ values. Then using a recursion with another Hyperbola Method on the rest, with $\phi = 1 * \mu$ , while making use of the fact that there are $\Theta(\sqrt{n})$ unique values of $\lfloor n/i \rfloor$ . See here! This balances out to $\mathcal{O} \left(\sum^{x^{1/3}} \sqrt{x/i} \right) = \mathcal{O}(x^{2/3})$ . In fact, through the recursion we get all the “important values”, the prefix sums up to $\phi(\lfloor n/i \rfloor)$ , for free.

So now using the Hyperbola Method to consolidate $I \tau$ , choose splitting point $pq = s$ :

\sum^s I\tau = \sum^{p} \phi(i) \sum^{s/i} \sigma(j) + \sum^{q} \sigma(i) \sum^{s/i} \phi(j) -\sum^{q} \sigma(i) \sum^{p} \phi(j)

As the computation of $\phi$ ’s prefix sums dominate, we can simply pick $p = q = \sqrt{s}$ while also not worrying about quickly obtaining all the “important values” of $\sigma$ , which only contributes $\mathcal{O}(\sqrt{n} \log n)$ .

Now to resolve

S = \sum_{d \leq \sqrt{n}} \phi(d) d^2 \sum_{\alpha \leq n/d^2} \alpha \tau(\alpha)

We can linear sieve $\phi$ for the other part in the $S$ sum at no asymptotic cost. Thus the total cost of this is

\sum_{d \leq \sqrt{n}} \mathcal{O}((n/d^2)^{2/3}) = \mathcal{O}(n^{2/3})

Although this is asymptotically optimal, in practice for values like $n \approx 10^{10}$ the constant factors come into play. Some experimentation shows that it’s better to decrease the precomputation of the linear sieves and rely more on the Hyperbola Method-induced recursion. $\blacksquare$

Further Improvements

I’m sure that many improvements can be made to this. I haven’t explored too much about it, but here’s some immediate ideas.

A quick skim of search results tell me prefix sums of $\phi$ and $\sigma$ can be calculated even quicker using more advanced methods that I can’t understand.
I have seen this remarkable blog, showing how to compute the “important values” of the Dirichlet product of two functions whose “important values” we already know in $\mathcal{O}(\sqrt{n} \text{ polylog } n)$ . Combined with the above the whole algorithm can also be $\mathcal{O}(\sqrt{n} \text{ polylog } n)$ .
A non-improvement: Although $\phi I^2$ , the other part of the multiplicand, is also multiplicative, the Hyperbola Method is limited by a $\sqrt{n}$ time. So trying to use it on the overall sum $S$ with a fast $\phi I^2$ prefix sum gains nothing as long as $I\tau = \phi * \sigma$ is used.
If anyone wants to correct me or share any improvements please do. With that being said I should implement blog commenting.