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