Multiplication of Dirichlet Series with Python

It is well known that polynomial multiplication can be much more efficiently using Fast Fourier Transform. I was searching for similar to use in case Dirichlet series. For example, if I have two functions $f(n)$ and $g(n)$, with Dirichlet series $F(s)$ and $G(s)$ respectively, what would be efficient way to get the summartory function of the Dirichlet convolution of the two functions?

Again, we could simply use the Dirichlet's Hyperbola method to get a very generic algorithm to solve this.

Let $f(n)$ and $g(n)$ be two arithmetic functions. Define $h(n)=(f*g)(n)$ as the Dirichlet convolution of $f$ and $g$. Also, let $F(n)$, $G(n)$ and $H(n)$ be the summatory functions for the respective functions.

Now using the formulas we saw in Dirichlet's Hyperbola method, I wrote the following Python code that does the job. The runtime of the following algorithm is about $O(n^{\frac{3}{4}})$

from timeit import default_timer

start = default_timer()

Lim = 1000000000
s = int(Lim ** 0.5)

Divs = [0] * (1 + s)

for k in range(1, 1 + s):
    temp = [1] if k == 1 else [1, k]
    stemp = int(k ** 0.5)
    for j in range(2, 1 + stemp):
        if k % j == 0:
            if j * j == k:
                temp += [j]
            else:
                temp += [j, k // j]
    Divs[k] = sorted(temp)


def DirMult(Losf, Highsf, Losg, Highsg, n):
    s = int(n ** 0.5)
    Highs = [0]
    Los = [0] * (1 + s)
    for k in range(1, 1 + s):
        temp = 0
        for j in Divs[k]:
            temp += (Losf[j] - Losf[j - 1]) * (Losg[k // j] - Losg[k // j - 1])
        Los[k] = Los[k - 1] + temp
        temp, tempn = 0, Lim // k
        u = int(tempn ** 0.5)
        v = tempn // (1 + u)
        mini, maxi = min(u, v), max(u, v)
        for j in range(1, 1 + mini):
            temp += (Losf[j] - Losf[j - 1]) * (Highsg[k * j] if k * j <= s else Losg[tempn // j]) + (Losg[j] - Losg[
                j - 1]) * (Highsf[k * j] if k * j <= s else Losf[tempn // j])
        for j in range(1 + mini, 1 + maxi):
            temp += (Losg[j] - Losg[j - 1]) * (Highsf[k * j] if k * j <= s else Losf[tempn // j])
        j = 1 + u
        temp -= Losg[u] * Losf[tempn // j]
        Highs += [temp]
    return Los, Highs


Losf, Highsf = [k for k in range(1 + s)], [0] + [Lim // k for k in range(1, 1 + s)]
Losg, Highsg = [k for k in range(1 + s)], [0] + [Lim // k for k in range(1, 1 + s)]

Lows, Highs = DirMult(Losf, Highsf, Losg, Highsg, Lim)

print(Highs[1])
print(default_timer() - start)

We'll see how to use this code. To use the function $\text{DirMult}$, we need four parameters. These parameters defined as below.

$\text{Losf}[k]=F(k)$, $\text{Highsf}[k]=F\left(\left\lfloor\frac{n}{k}\right\rfloor\right)$, $\text{Losg}[k]=G(k)$, and $\text{Highsg}[k]=G\left(\left\lfloor\frac{n}{k}\right\rfloor\right)$

We get two outputs and they too are interpreted the same way.

$\text{Los}[k]=H(k)$, $\text{Highs}[k]=H\left(\left\lfloor\frac{n}{k}\right\rfloor\right)$

In the example code above, we have $F(s)=G(s)=\zeta(s)$. Therefore, the output we get is the Divisor Summatory function.

Note that as the algorithm is made for an almost no-brainer approach. It doesn't take advantage of the properties of the functions that we are multiplying and hence the increased runtime.

I'll see you in the next post with Division of Dirichlet Series.


Until then
Yours Aye
Me