Partition Numbers and a code to generate one in Python

Author: Kazi Abu Rousan

The pure mathematician, like the musician, is a free creator of his world of ordered beauty.
Bertrand Russell

Today we will be discussing one of the most fascinating idea of number theory, which is very simple to understand but very complex to get into. Today we will see how to find the Partition Number of any positive integer $n$.

Like every blog, our main focus will be to write a program to find the partition using python. But today we will use a very special thing, we will be using Math Inspector - A visual programming environment. It uses normal python code, but it can show you visual (block representation) of the code using something called Node Editor. Here is a video to show it's feature.

Beautiful isn't it? Let's start our discussion.

What are Partition Numbers?

In number theory, a Partition of a positive integer $n$, also called an Integer Partition, is a way of writing $n$ as a sum of positive integers (two sums that differ only in the order of their summands are considered the same partition, i.e., $2+4 = 4+ 2$). The partition number of $n$ is written as $p(n)$

As an example let's take the number $5$. It's partitions are;

$$5 = 5 +0\ \ \ \ \\= 4 + 1\ \ \ \ \\ = 3 + 2\ \ \ \ \\= 3 + 1+ 1\ \ \\= 2 +2 + 1\ \ \\= 2+1+1+1 \ \\= 1+1+1+1+1$$

Hence, $p(5) = 7$. Easy right?, Let's see partitions of all numbers from $0$ to $6$.

Note: The partition number of $0$ is taken as $1$.

Partition of numbers from $0$ to $6$.

See how easy it is not understand the meaning of this and you can very easily find the partitions of small numbers by hand. But what about big numbers?, like $p(200)$?

Finding the Partition Number

There is no simple formula for Partition Number. We normally use recursion but if you want just a formula which generate partition number of any integer $n$, then you have to use this horrifying formula:

Scary!!!-- Try using this to find $p(2)$

Just looking at this can give you nightmare. So, are there any other method?, well yes.

Although, there are not any closed form expression for the partition function up until now, but it has both asymptotic expansions(Ramanujan's work) that accurately approximate it and recurrence relations(euler's work) by which it can be calculated exactly.

Generating Function

One method can be to use a generating function. In simple terms we want a polynomial whose coefficient of $x^n$ will give us $p(n)$. It is actually quite easy to find the generating function.

$$ \sum_{n =0}^{\infty} p(n)x^n = \Pi_{j = 0}^{\infty}\sum_{i=0}^{\infty}x^{ji} = \Pi_{j = 1}^{\infty}\frac{1}{1-x^j}$$

This can also written as:

$$ \sum_{n =0}^{\infty} p(n)x^n = \frac{1}{(1-x)}\frac{1}{(1-x^2)}\frac{1}{(1-x^3)}\cdots = (1+x^1+x^2+\cdots+x^r+\cdots) (1+x^2+x^4+\cdots+x^{2r}+\cdots)\cdots $$

Each of these $(1+x^k+x^{2k}+\cdots+x^{rk}+\cdots)$ are called sub-series.

We can use this to find the partition number of any $n$. Let's see an example for $n=7$. There are 2 simple steps.

First write the polynomial such that each sub-series's maximum power is smaller or equals to $n$. As an example, if $n=7$, then the polynomial is $(1+x+x^2+x^3+x^4+x^5+x^6+x^7)(1+x^2+x^4+x^6)(1+x^3+x^6)(1+x^4)(1+x^5)(1+x^6)(1+x^7)$.
Now, we simplify this and find coefficient of $x^n$. In this process we find the partition numbers for all $r<n$.

So, for our example, we can just expand the polynomial given in point-1. You can do that by hand. But I have used sympy library.

from sympy import *
x = symbols('x')
a = expand((1+x+x**2+x**3+x**4+x**5+x**6+x**7)*(1+x**2+x**4+x**6)*(1+x**3+x**6)*(1+x**4)*(1+x**5)*(1+x**6)*(1+x**7))
print(a)
#Output: x**41 + x**40 + 2*x**39 + 3*x**38 + 5*x**37 + 7*x**36 + 11*x**35 + 15*x**34 + 18*x**33 + 24*x**32 + 30*x**31 + 37*x**30 + 44*x**29 + 52*x**28 + 57*x**27 + 64*x**26 + 71*x**25 + 77*x**24 + 81*x**23 + 84*x**22 + 84*x**21 + 84*x**20 + 84*x**19 + 81*x**18 + 77*x**17 + 71*x**16 + 64*x**15 + 57*x**14 + 52*x**13 + 44*x**12 + 37*x**11 + 30*x**10 + 24*x**9 + 18*x**8 + 15*x**7 + 11*x**6 + 7*x**5 + 5*x**4 + 3*x**3 + 2*x**2 + x + 1

The coefficients of $x^k$ gives $p(k)$ correctly up to $k = 7$. So, from this $p(7) = 15$, $p(6)=11$ and so on. But for $n>7$, the coefficients will be less than $p(n)$, as an example $p(8)=22$ but here the coefficient is $18$. To calculate $p(8)$, we need to take sub-series with power $8$ or less.

Although this method is nice, it is quite lengthy. But this method can be modified to generate a recurrence relation.

Euler's Pentagonal Number Theorem

Euler's Pentagonal Theorem gives us a relation to find the polynomial of the product of sub-series. The relation is:

$$ \Pi_{i = 1}^{\infty }(1-x^n) = 1 + \sum_{k =1}^{\infty } (-1)^k \Big(x^{k\frac{(3k+1)}{2}} + x^{k\frac{(3k-1)}{2}}\Big) $$

Using this equation, we get the recurrence relation,

$$p(n) = p(n-1) + p(n-2) - p(n-5) - p(n-7) + \cdots = \sum_{k=1}^{n} \Bigg( p\Bigg(n-\frac{k(3k-1)}{2}\Bigg) + p\Bigg(n-\frac{k(3k+1)}{2}\Bigg)\Bigg) $$

The numbers $0, 1, 2, 5, 7, 12, 15, 22, 26, 35, 40,\cdots $ , are called pentagonal numbers. We generate these numbers by $\frac{k(3k+1)}{2}$ and $\frac{k(3k-1)}{2}$.

Let's try to find $p(7)$, using this.

Using this simple algorithm, we can write a code to get partition number.

def par(n):
	summ = 0
	if n == 0 or n == 1:
		result = 1
	else:
		for k in range(1,n+1):
			d1 = n - int((k*(3*k-1))/2)
			d2 = n - int((k*(3*k+1))/2)
			sign = pow(-1,k+1)
			summ = summ + sign*(par(d1) + par(d2))
		result = summ
	return result

This is the code to generate partition number. This is not that efficient as it uses recurrence. For big numbers it will take huge amount of time. In math-inspector, we can see it's block diagram.

Here the par block represent the function par(n). Whenever i creates a variable $n1=6$ or $n3=10$, it creates a circle. To apply par on $n3$, i just need to connect their nodes. And to show the output, add the other node to output node.

So, we have seen this simple function. Now, let's define one using the euler's formula but in a different manner. To understand the used algorithm watch this: Explanation

def partition(n):
    odd_pos = []; even_pos= []; pos_d = []
    for i in range(1,n+1):
        even_pos.append(i)
        odd_pos.append(2*i+1)
    m = 0; k = 0
    for i in range(n-1):
        if i % 2 == 0:
            pos_d.append(even_pos[m])
            m += 1
        else:
            pos_d.append(odd_pos[k])
            k += 1
    initial = 1; pos_index = [initial]; count = 1
    for i in pos_d:
        d = initial + i
        pos_index.append(d)
        initial = d
        count += 1
        if count > n:
            break
    sign = []
    for i in range(n):
        if i % 4 == 2 or i % 4 == 3:
            sign.append(-1)
        else:
            sign.append(1)
    pos_sign = []; k = 0
    for i in range(1,n+1):
        if i in pos_index:
            ps = (i,sign[k])
            k = k + 1
            pos_sign.append(ps)
        else:
            pos_sign.append(0)
    if len(pos_sign) != n:
        print("Error in position and sign list.")
    partition = [1]
    f_pos = []
    for i in range(n):
        if pos_sign[i]:
            f_pos.append(pos_sign[i])
        partition.append(sum(partition[-p]*s for p,s in f_pos))
    return partition

This is the code. This is very efficient. Let's try running this function.

Note: This function returns a list which contains all $p(n)$ in the range $0$ to $n$.

As you can see in the program, In[3] gives a list $[p(0),p(1),p(2),p(3),p(4),p(5)]$. Hence, to get the $p(n)$ using this function, just add $[-1]$ as it gives the last element of the list.

We can use a simple program to plot the graph.

import matplotlib.pyplot as plt
n = 5
x_vals = [i for i in range(n+1)]
y_vals = partition(n)
plt.grid(color='purple',linestyle='--')
plt.ylabel("Partition Number")
plt.xlabel("Numbers")
plt.step(x_vals,y_vals,c="Black")
plt.savefig("Parti_5.png",bbox_inches = 'tight',dpi = 300)
plt.show()

Output is given below. The output is a step function as it should be.

This is all for today. I hope you have learnt something new. If you find it useful, then share.