Cheenta

Academy for Gifted Students

How Cheenta works to ensure student success?

Explore the Back-Story** Author**:

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

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

In number theory, a* Partition of a positive integer *$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*:

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)$?

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:

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.

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)$

$$ \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 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.

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.

** Author**:

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

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

In number theory, a* Partition of a positive integer *$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*:

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)$?

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:

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.

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)$

$$ \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 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.

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.

Cheenta is a knowledge partner of Aditya Birla Education Academy

Advanced Mathematical Science. Taught by olympians, researchers and true masters of the subject.

JOIN TRIALAcademic Programs

Free Resources

Why Cheenta?

Online Live Classroom Programs

Online Self Paced Programs [*New]

Past Papers

More

kudos. that is some very efficient code right there. mine is way faster than the recustion at the beginning but also way slower than you final version.

UPDATE: I've now implemented the mathologer algorithm myself and it runs even more efficent than yours. I think that's because mine takes about half as many lines and accordingly has way fewer loops in it.

Nice!!...

I was actually not looking into efficiency but only in the algorithm. Maybe you can share your code with others.. if it's alright

Sorry, I totally overlooked your reply. Here's my code:

def Part(n):

# initializing an indexlist IL of the form [1, 2, 5, 7, 12, ...], using L = [1, 3, 2, 5, 4, ...]

L = [1]; IL = [1]; k = 2

while IL[-1] n:

IL = IL[:-1]

# intitializing a list P of partition numbers, then calculating new ones using IL

P = [1, 1]; k = 2; sublist = IL[:2]

while len(P) < n+1:

if k < len(IL) and IL[k] <= len(P):

sublist.append(IL[k])

k += 1

summ = 0; count = 0

for i in sublist:

if count%4 < 2:

summ += P[-i]

else:

summ -= P[-i]

count += 1

P.append(summ)

return(P[-1])

For n = 50 000 it's about 10 times as fast as yours on my machine (8s vs 80s). My code hits that mark at about n = 200 000 at which point the number of partitions of n has gone up from about 3.6*10^243 to 1.1*10^492.

Also, I apparently can't reply directly to your comment. 🙁

Aaaand I just noticed that my code and indentation didn't make it through the comment processing alive. The first few rows should be

def PartE(n):

# initializing an indexlist IL of the form [1, 2, 5, 7, 12, ...], using L = [1, 3, 2, 5, 4, ...]

L = [1]; IL = [1]; k = 2

while IL[-1] n:

IL = IL[:-1]

... and so on. I have no idea how that happened - I just copy-pasted it.

Here are the indentation-levels of each line of the complete code (starting with def() and excluding comments):

0, 1, 1, 2, 2, 2, 2, 2, 1, 2, 1, 1, 2, 3, 3, 2, 2, 3, 4, 3, 4, 3, 2, 1,.

I hope that's all.

Goddamnit! There's a few lines it just doesn't show when I post it. One more try:

def PartE(n):

A = timer()

# initializing an indexlist IL of the form [1, 2, 5, 7, 12, ...], using L = [1, 3, 2, 5, 4, ...]

L = [1]; IL = [1]; k = 2

while IL[-1] n:

IL = IL[:-1]

... and so on. There should be a total of 11 lines of code above (not counting comments) if it works correctly.

I resign.