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

Divisor functions - Continued..

This post continues the discussion we had in my previous post.

I discussed only $k=0$ case before dismissing the others as uninteresting. They aren't after all. After some further work, I think I've figured out how the first DGF given above generalized to the other cases as well.

Following the same notations given before we have,

$\displaystyle\sum_{n=1}^\infty\frac{\sigma_k(m\cdot n)}{n^s}=\zeta(s)\zeta(s-k) \prod_{p_j|m}\left(1+(p_j^{k}+p_j^{2k}+p_j^{3k}+\cdots+p_j^{m_jk})(1-p_j^{-s})\right)$

Quite simple after all!! But I don't think the second DGF simplifies as nicely. That should be one for another post. But there is something else.

The first formula given in the previous post works even when $n$ is replaced by $n^2$. We are gonna make that replacement because of the closed form DGF we have for $\sigma_k(n^2)$.

Again following the same notations, this time we have,

$\displaystyle\sum_{n=1}^\infty\frac{\sigma_k(m\cdot n^2)}{n^s}=\frac{\zeta(s)\zeta(s-2k)\zeta(s-k)}{\zeta(2s-2k)} \prod_{p_j|m}\left(1+(p_j^{k}+p_j^{2k}+\cdots+p_j^{m_jk})\left(\frac{1-p_j^{-s}}{1+p_j^{-(s-k)}}\right)\right)$

As before, these DGFs can be used in combination with Dirichlet's Hyperbola method to find their summatory functions. Hope you find these interesting and worth your while. See ya soon.

Till then
Yours Aye,
Me

Divisor functions

Oh... I didn't realize I'vnt posted anything this year. Anyway, here's my first post of the year - my thoughts on the Divisor function.

The Divisor functions are well known and have been studied extensively in Mathematics. Their relation to the Partition function is one of best identities I've seen. I was playing around with them as part of my effort to solve a problem and came up with some identities that I can't find anywhere on the Web. So this post will be about them.

The Divisor function is known to be multiplicative, but not completely multiplicative. But I needed an some formulas where I wanted to find the value of the Divisor functions of product of two numbers that are not co-prime. After some extensive calculations, I ended with the following:

Let $m=p_1^{m_1}p_2^{m_2}p_3^{m_3}\cdots$ and $n=p_1^{n_1}p_2^{n_2}p_3^{n_3}\cdots$. Then,

$\displaystyle\frac{\sigma_k(m\cdot n)}{\sigma_k(n)}=\prod_{p_j|m}\left(1+\frac{p_j^k+p_j^{2k}+p_j^{3k}+\cdots+p_j^{m_jk}}{1+p_j^{-k}+p_j^{-2k}+p_j^{-3k}+\cdots+p_j^{-n_jk}}\right)$

We can easily check a couple of easy cases. For example, when $n=1$, all of the $n_j$'s become $0$ and the above formulae reduces to the definition of Divisor function. Similarly, when $m$ and $n$ are co-prime, for any prime that divides $m$, the exponent of that prime in $n$ will be $0$, again reducing to a form that shows that divisor functions are multiplicative.

I found the case $k=0$ especially very interesting. We'll see why in a short while. First of all, when $k=0$, the above formula reduces to,

$\displaystyle\frac{\sigma_0(m\cdot n)}{\sigma_0(n)}=\left(1+\frac{m_1}{n_1+1}\right)\left(1+\frac{m_2}{n_2+1}\right)\left(1+\frac{m_3}{n_3+1}\right)\cdots$

As $\sigma_0(n)=(1+n_1)(1+n_2)(1+n_3)\cdots$, we can write

$\sigma_0(m\cdot n)=\sigma_0(n)+m_1\sigma_{0,p_1}(n)+m_2\sigma_{0,p_2}(n)+\cdots+m_1m_2\sigma_{0,p_1p_2}(n)+\cdots+m_1m_2\cdots m_k\sigma_{0,p_1p_2\cdots p_k}(n)+\cdots$

where $\sigma_{0,B}(n)$ is (defined only for square-free $B$) a function that counts the number of divisors of $n$ that are not divisible by any of the prime factors of $B$. Note that we have eliminated  everything related to prime factorization of $n$.

Actually this is why I said the case $k=0$ is interesting. Had $k$ been greater than $0$, we would have those annoying prime factors of $n$ in the above expression. Now with what we have above we can make an interesting Dirichlet Generating function.

Given $m=p_1^{m_1}p_2^{m_2}p_3^{m_3}\cdots$, we have

$\displaystyle\sum_{n=1}^\infty\frac{\sigma_0(m\cdot n)}{n^s}=\zeta^2(s) \prod_{p_j|m}\left(1+m_j(1-p_j^{-s})\right)$

Even better, for a given $m$, if we define a function $f(n)$ as

$f(n)=\displaystyle\sum_{d|m}\sigma_0(d\cdot n)$

again we end up with a nice DGF. This we have,

$\displaystyle\sum_{n=1}^\infty\frac{f(n)}{n^s}=\sigma_0(m)\zeta^2(s)\prod_{p_j|m}\left(1+\frac{m_j}{2}(1-p_j^{-s})\right)$

All of the above DGFs in combination with Dirichlet's Hyperbola method can be used to find their respective summatory function.

We can see how this generalizes in the next post. See ya soon.

Till then
Yours Aye,
Me

SquareFree numbers with exactly $k$ factors

A positive integer is squarefree if and only if it is not divisible by any perfect square other than 1. That makes them interesting. This post will be about finding the DGF squarefree numbersthat has exactly $k$ factors. We'll use symmetric polynomials to our advantage.

First, the Dirichlet Generating function of square numbers is quite straight forward. By the definition of square free numbers, every prime can be present atmost once in its prime factorization. Let $\xi_2(s)$ be the DGF of square free numbers. Then it is easy to see that,

$\xi_2(s)=\displaystyle\prod_p(1+p^{-s})$

Now, let's review about symmetric polynomials. We need two of the most fundamental symmetric polynomials. The Elementary symmetric polynomials (denoted by $e_k(x_1,x_2,\cdots)$) and the Power sum symmetric polynomials (denoted by $p_k(x_1,x_2,\cdots)$.

Consider these polynomials with finitely many variables. Then,

$e_1=x_1+x_2+x_3+\cdots$
$e_2=x_1x_2+x_1x_3+\cdots+x_2x_3+x_2x_4+\cdots$ and so on.

Similarly,

$p_k=x_1^k+x_2^k+\cdots$

Now notice something interesting. Let $x_k$ be the inverse of the $k$th prime raised to the power $s$. That is,

$x_1=2^{-s}$, $x_2=3^{-s}$, $x_3=5^{-s}$, and so on.

It is now very easy to see that $e_k =e_k(s)$ becomes the Dirichlet generating function of the square free numbers with exactly $k$ factors. It is also clear that the

$p_k =p_k(s)=\zeta_P(ks)$

where $\zeta_P(s)$ is the Prime Zeta function.

Therefore by the fundamental realtion connecting the Elementary symmetric polymials and the Power sum symmetric polynomials, we know

$e_k=\displaystyle\frac{1}{k}\sum_{j=1}^k(-1)^{j-1}e_{k-j}p_j$

This equation can now be directly translated as,

$e_k(s)=\displaystyle\frac{1}{k}\sum_{j=1}^k(-1)^{j-1}e_{k-j}(s)\zeta_P(js)$

which gives the sought after DGF in a recursive fashion. We'll use this DGF in combination with the Dirichlet Hyperbola method in the next post to find the number of squarefree numbers with exactly $k$ factors.


Until then
Yours Aye
Me

Irrelevance of Bias

There were quite some places in mathematical problems in which a condition imposed to make the problem supposedly difficult actually turns out to be irrelevant. The first instance where I saw something like this in the 'Compass and Straight edge constructions' in which the concept of collapsible compass actually turns out to be irrelevant. This fact is called the Compass Equivalence Theroem.

In this post, we see one such 'irrelevance' in probability theory. I've often wondered whether actual coin tosses are really fair and, if not, how to make them fair. I did some searches and quite surprised to see that this turns out be an easy problem with a simple solution (atleast in theory).

Consider you've a coin with a probability $p$ of turning up heads. Say you want to create a game out of this coin such that each player has a $50\%$ probability of winning. This can be done by the following procedure:

Toss the coin twice. If the results of both the tosses are same, ignore the two tosses and toss again. Repeat until you've different results in both the tosses. Now if the result is $HT$, the first player wins. Else, the second player does.

This works because $\mathbb{P}(HT)=p(1-p)$ whereas $\mathbb{P}(TH)=(1-p)p$.

As both has the same probability of happening and these are the only two outcomes. (Note that, we have toss the coin in 'twos' and we can't combine the intermediate parts of two tosses).

The beauty of this procedure is two-fold. First, the entire argument boils down to the commutative property of multiplication and the second is that we don't even have to know what the unfair probability is. Incredible!!

Now that we have dealt with making an unfair coin to come up with a fair result, let's look at the other way. Now we have a fair coin; how can we make game in which the winning probability of someone is $\alpha$.

This again has a simple solution. Cast the real number $\alpha$ in binary terms. Flip the fair coin until it gives the first Head. Now the player wins, if the $n$th digit in the binary expansion is $1$ and loses otherwise.

Say for example, $\alpha=\frac{3}{4}=(0.11)_2$. Now the player wins if the first Heads appears in first or the second toss and loses otherwise.

$\mathbb{P}(Win)=\mathbb{P}(\text{First Head in }1^{st}\text{ toss})+\mathbb{P}(\text{First Head in }2^{nd}\text{ toss})=\frac{1}{2}+\frac{1}{4}=\frac{3}{4}$.

All this means is that we can take any (un)fair coin with arbitrary probability $p$ and use it to create an arbitrary probability $\alpha$ of winning for some player. In other words, the bias doesn't matter.

The only difference would be that would be needing more and more tosses for the game to conclude depending the unfairness of the coin and winning probability. With that in mind, we will do something more.

Say we want to create a equal probability of winning between three players with an unfair coin. We can directly generalize the first argument and achieve this with the following procedure.

Toss the coin thrice. If all the results are the same, repeat. Else, either exactly one Heads has appeared or exactly one Tails. Assume WLOG, exactly one Heads. Then the the First player wins if the Heads happened in the first toss, the Second player wins if it did in the second toss else the third player wins.

This procedure can be applied to $4$ players, $5$ players, .. Oh wait, for $6$ players, we can do better. Instead of making six tosses every time and waiting for exactly one Heads (or one Tails), we can do it with four tosses.

Toss the coin $4$ times until exactly two Heads appear. Discard other batches of four tosses. Now the First player wins when the Heads occur in positions $(1,2)$, Second wins in positions $(1,3)$, Third wins in $(1,4)$, Fourth in $(2,3)$, Fifth in $(2,4)$ and the Sixth in $(3,4)$.

Similarly, among $10$ players, instead of tossing the coin in batches of ten, we can toss it in batches of five and wait for exactly two Heads (or two Tails) to appear.

By now one can see that we can use the unfair coin and create a game with equal probability of winning among $\binom{n}{k}$ players. Though I'vnt done much analysis further, a judicious choice depending on $n$, $k$ and $p$ will lead to shorter games than the other choices available. Quite neat, don't you think?

Yours Aye

Me

Bridge between Generating functions

I've always wondered about the question of going from OGF to the EGF or vice versa. Turns out this is not too complicated, though it needs some precise mathematical justification. That aside, I find the 'transition' between the generating functions very useful in some cases and very beautiful in others. Anyways, that will be discussion in this post. First, the simpler case of moving from EGF to OGF.

Note that these are just simple manipulations which anyone can do and I do it here just for the fun of doing it.

Let $f_e(x)$ be the EGF of a sequence. That is, 

$f_e(x)=\displaystyle a_0+a_1\frac{x}{1!}+a_2\frac{x^2}{2!}+a_3\frac{x^3}{3!}+\cdots$

Now using the definition of the Gamma integral, one can easily see that,

$f_o(x)=\displaystyle\int\limits_{0}^{\infty}e^{-t}f_e(tx)\, dt$

Technically, we should say that this equality holds provided either the sum on the LHS or the integral on the RHS converges.

For the simplest non-trivial sequence of $a_k=1, k \geq 0$, we know that $f_e(x)=e^x$ and $f_o(x)=1/(1-x)$ which readily satisfies the above equation.

Now for the opposite direction, we make use of the inverse Laplace transform of $t^n$. That is,

$\displaystyle\mathcal{L}^{-1}\left\{\frac{1}{s^{n+1}}\right\}=\frac{t^n}{n!}$

Again simple manipulations yield, $\displaystyle\mathcal{L}^{-1}\left\{\frac{1}{s}f_o\left(\frac{1}{s}\right)\right\}=f_e(t)$

Here again, the trivial case can be seen to work well here. However, we can one another nice application. Consider the OGF of the number of ways to partition a number $n$ with each part atmost $3$. That is,

$f_o(x)=\displaystyle\frac{1}{(1-x)(1-x^2)(1-x^3)}$

This generating function by itself does not help us in finding the number of partitions directly. Let's find the EGF of the same and see how useful it is.

$\displaystyle\frac{1}{s}f_o\left(\frac{1}{s}\right)=\frac{s^5}{(s-1)(s^2-1)(s^3-1)}$

Using Wolfram Alpha to compute the Inverse Laplace Transform of the above we get,

$f_e(t)=\displaystyle \frac{1}{12}t^2e^t+\frac{7}{12}te^t+\frac{47}{72}e^t+\frac{2}{9}Re[e^{\omega t}]+\frac{1}{8}e^{-t}$

where $\omega$ is the cubic root of unity. Knowing that $n![t^n]t^ke^t=n(n-1)\cdots(n-k+1)$, we have

$n![t^n]f_e(t)=\displaystyle \frac{n(n-1)}{12}+\frac{7n}{12}+\frac{47}{72}+\frac{2}{9}Re[\omega^n]+\frac{(-1)^n}{8}$

This already is a cool closed form expression for partition of $n$ with each part atmost $3$. If we also note that $Re[\omega^n]$ is periodic with period $3$, we now have an expression which can also be shown to count the number of triangles with perimeter $n$.

Hope you enjoyed this. I've a couple of ideas for my next post.

Until then

Yours Aye

Me