Division of Dirichlet Series with Python

We saw how to do Dirichlet Series multiplication in the previous post. This post is about dividing two Dirichlet series. Using pretty much the same notation as before, here we are concerned about the following.

$H(s)=\displaystyle\frac{F(s)}{G(s)}$

We can use the following Python code to get the Summatory function $H(n)$.

Note: We know that Dirichlet inverse exists if and only if $g(1)\ne0$. In fact, for most of the multiplicative functions $g(1)=1$. So I have made this assumption in the code. If it is not, changes have to made accordingly.

from timeit import default_timer

start = default_timer()

n = 1000000000
s = int(n ** 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 DirDiv(Lowsn, Highsn, Lowsd, Highsd, n):
    s = int(n ** 0.5)
    Lows, Highs = [0] * (1 + s), [0] * (1 + s)
    Lows[1] = 1 // (Lowsd[1] - Lowsd[0])
    for k in range(2, 1 + s):
        for j in Divs[k][:-1]:
            Lows[k] -= (Lows[j] - Lows[j - 1]) * (Lowsd[k // j] - Lowsd[k // j - 1])
        Lows[k] += Lowsn[k] - Lowsn[k - 1]
        Lows[k] *= Lows[1]
        Lows[k] += Lows[k - 1]
    for k in range(s, 0, -1):
        tempn = n // k
        res = Highsn[k]
        u = int(tempn ** 0.5)
        v = tempn // (1 + u)
        for i in range(1, 1 + v):
            res -= (Lows[i] - Lows[i - 1]) * (Highsd[k * i] if k * i <= s else Lowsd[tempn // i])
        for i in range(2, 1 + u):
            res -= (Lowsd[i] - Lowsd[i - 1]) * (Highs[k * i] if k * i <= s else Lows[tempn // i])
        i = 1 + u
        res += Lows[tempn // i] * Lowsd[u]
        Highs[k] = res
    return Lows, Highs


Lowsn, Highsn = [0] + [1] * s, [0] + [1] * s
Lowsd, Highsd = [k for k in range(1 + s)], [0] + [n // k for k in range(1, 1 + s)]
Lows, Highs = DirDiv(Lowsn, Highsn, Lowsd, Highsd, n)
print(Highs[1])
print(default_timer() - start)

As in the multiplication case, here again the 'Lows' and 'Highs' mean the same. In the sample code given above, we have $F(s)=1$ and $G(s)=\zeta(s)$. This means $H(s)$ is the Dirichlet series of the Moebius function.

As a result, we get the values of the Merten's function as the output. Here again, the runtime of the algorithm is about $O(n^{\frac{3}{4}})$.

Until then
Yours Aye
Me

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