Prime Partition

Question

A prime partition of a number is a set of primes that sum to the number. For instance, {2 3 7} is a prime partition of $12$ because $2 + 3 + 7 = 12$. In fact, there are seven prime partitions of $12$: {2 2 2 2 2 2}, {2 2 2 3 3}, {3 3 3 3}, {2 2 3 5}, {2 5 5}, {2 3 7}, and {5 7}. The number of prime partitions of a number is given by A000607.

I want to calculate the number of prime partitions of a number. From reading the text of A000607 it seems to me that the formula $\prod\limits_{p \in \pi(n)} \frac1{1 - n^p}$ should work. But it doesn't. Consider the case of $n=12$. The primes less than $12$ are $2, 3, 5, 7$, and $11$. Thus the formula computes $\frac1{1-12^2} \times \frac1{1-12^3} \times \frac1{1-12^5} \times \frac1{1-12^7} \times \frac1{1-12^{11}}$ = $\frac1{-143} \times \frac1{-1727} \times \frac1{-248831} \times \frac1{-35831807} \times \frac1{-743008370687}$ = $\frac{-1}{1636045119596820253743372240719}$ which is obviously incorrect.

How can I compute the number of prime partitions of a number?

Answer 1

You need to learn a bit about generating functions. The text associated with A000607 means the following. For each prime $p$ expand the function $1/(1-x^p)$ as a power series: $$ \frac1{1-x^p}=1+x^p+x^{2p}+x^{3p}\cdots=\sum_{k=0}^\infty x^{kp}. $$ Call that series $f_p(x)$. Then you multiply these power series together, and identify the coefficient of $x^n$. That coefficient is then the desired value of the prime partition function. Let's do this for $n=12$. To that end we can throw away all the terms of degree $>12$. I denote those with three dots. So start with $$ f_2(x)=1+x^2+x^4+x^6+x^8+x^{10}+x^{12}+\cdots $$ Multiplying this with $f_3(x)=1+x^3+x^6+x^9+x^{12}+\cdots$ gives $$ \begin{aligned} f_2(x)f_3(x)=&f_2(x)+x^3f_2(x)+x^6f_2(x)+x^9f_2(x)+x^{12}f_2(x)+\cdots\\ =&1+x^2+x^3+x^4+x^5+2x^6+x^7+2x^8+2x^9+2x^{10}+2x^{11}+3x^{12}+\cdots \end{aligned} $$ At this point you should check that the coefficient of $x^k$ counts the number of ways of writing $k$ as a sum of twos and threes.

Next we add $p=5$ to the mix, and multiply the above with $f_5(x)=1+x^5+x^{10}+\cdots$ and get $$ \begin{aligned} &f_2(x)f_3(x)f_5(x)\\ =&1+x^2+x^3+x^4+2x^5+2x^6+2x^7+3x^8+3x^9+4x^{10}+4x^{11}+5x^{12}+\cdots \end{aligned} $$

Next we multiply this with $f_7(x)=1+x^7+\cdots$ and get $$ \begin{aligned} &f_2(x)f_3(x)f_5(x)f_7(x)\\ =&1+x^2+x^3+x^4+x^5+2x^6+3x^7+3x^8+4x^9+5x^{10}+5x^{11}+7x^{12}+\cdots\\ \end{aligned} $$

As a laste step we multiply this with $f_{11}(x)=1+x^{11}+\cdots$ to end with $$ \begin{aligned} &f_2(x)f_3(x)f_5(x)f_7(x)f_{11}(x)\\ =&1+x^2+x^3+x^4+x^5+2x^6+3x^7+3x^8+4x^9+5x^{10}+6x^{11}+7x^{12}+\cdots\\ \end{aligned} $$

Here the term $7x^{12}$ appears. Primes $p>12$ won't affect the term of degree $12$, so at long last that tells us that there are seven prime partitions of $n=12$.

Answer 2

I've decided to turn one of my comments into a full answer, in the interest of keeping this thread self-contained. The formulae needed are from here (formulae 5-10).

First, we define the sum of prime factors function, $\mathrm{sopf}(n)$ (sequence A008472 in the OEIS):

$$\mathrm{sopf}(n)=\sum_{p\in\mathbb P,\;p\mid n}p$$

In Mathematica: PrimeDivisorSum[n_Integer] := DivisorSum[n, Identity, PrimeQ]. (On some other language and assuming you have a pre-cached array of primes, just check for the primes $\leq n$ that are divisors, and add them all up.)

The prime partition counting function, which I'll denote in this answer by $\kappa(n)$, then satisfies the following recursion:

$$\kappa(n)=\frac1{n}\left(\mathrm{sopf}(n)+\sum_{j=1}^{n-1} \mathrm{sopf}(j)\kappa(n-j)\right)$$

with the initial condition $\kappa(1)=0$. In effect, we are saying that $\kappa(n)$ is the Euler transform of the function $[n\in\mathbb P]$, where $[p]$ is an Iverson bracket.

In Mathematica:

PrimePartitionCount[1] = 0; 
PrimePartitionCount[n_Integer] := PrimePartitionCount[n] =
   (PrimeDivisorSum[n] +
    Sum[PrimeDivisorSum[k] PrimePartitionCount[n - k], {k, n - 1}])/n

---

If you want to see the actual partitions themselves, however, you will need to solve the Frobenius coin problem.

Answer 3

You are misunderstanding the text of OEIS's sequence A000607, I think. The function $$ g(x) = {\prod}_{p}\frac{1}{1-x^{p}} = \frac{1}{1-x^2}\times\frac{1}{1-x^3}\times\frac{1}{1-x^5}\times{...} $$ is the generating function for the sequence. The number of prime partitions of $n$ is given by the coefficient of $x^n$ in the Taylor series expansion of $g(x)$ around $x=0$: that is, $$ a_n = \frac{1}{n!}\frac{\partial^n}{\partial{x}^n}g(x)\Big\vert_{x=0}. $$ This doesn't give a practical way to actually calculate the coefficients, though. For that, a dynamic programming approach can easily be formulated (e.g., by recursively evaluating $b_{n,p}$, the number of prime partitions of $n$ using primes no greater than $p$, and setting $a_n=b_{n,n}$).

Answer 4

This is a misunderstanding. The formula you quote is the generating function for the numbers you're looking for. It would more usually be written as

$$\prod_p\frac1{1-x^p}\;,$$

where the product runs over all primes $p$, and then the number of partitions of $n$ into prime parts is the coefficient of $x^n$ in the power series for this function. However, this is not meant as a practical way to compute these numbers, more like an efficient way to organize the information about these numbers.

The MathWorld article on prime partitions describes how to compute the numbers of prime partitions using an Euler transform.

Answer 5

https://www.dropbox.com/s/6diyuy3dijqkjvd/image1.jpeg?dl=0

See above link for image of a table for calculating recursively a(n)=A000607(n), the number of partitions of n into primes. It is easy to use and involves no complicated formulae. The first column is n, the second a(n), and subsequent columns record the contributions to a(n) associated with different primes, denoted in the top row by p'. Thus a(n) is the sum of the numbers recorded in the boxes (n,p') in the nth row, p denoting primes between 2 and n. As worked example here is how to calculate a(31). Note that the box at (31,31) is already filled in with number 1 (denoting that 31 is a prime partition of itself). To find the value for (31,p), count up to the box p squares above and add together the value in that box with those in all boxes in the same row to the right. Thus in box (31,2) the value is 67+14+4+1+1=87. Likewise in box (31,3) the value is 12+3+1+1=17. Similarly (31,5)=4, (31,7)=2, (31,31)=1 and the rest are zeroes. The sum of these numbers is a(n)= 87+17+4+2+1=111. Bingo! Enter these numbers in the table and they contribute to calculations for higher values of a(n).

If the count up from any box runs up through the top row then that value is 0. The 1 in the top row is "floating", that is to say can move to any box on the first row to attribute its value in any count ending on that square.

I don't know if this method is generally known, contact me if you are interested in how it works David.

Answer 6

Let $b(n,p)$ be the number of partitions of $n$ into pieces each of which is positive, is prime, and is less than or equal to $p$. Let $pp(k)$ be largest prime less than or equal to $k$.
$$ pp(k) = \begin{cases} k , & k \text{ prime} \\ \max (\mathrm{Primes} \cap [1,k]) , & k \text{ composite} \end{cases} $$ Since the only numbers that can be written as a partition into primes $\leq 2$ are even numbers, \begin{align*} b(n,2) = \begin{cases} 0 , & n \text{ odd} \\ 1 , & n \text{ even} \end{cases} \text{.} \end{align*} Now suppose $p > 2$ is prime. Then \begin{align*} b(n,p) &= b(n - 0p,pp(p)) + b(n-1p, pp(p)) + \cdots + b \left( n- \left\lfloor \frac{n}{p} \right\rfloor p, pp(p) \right) \\ &= \sum_{k=0}^{\left\lfloor n/p \right\rfloor} b(n - k p, pp(p)) \text{.} \end{align*} Otherwise $p > 2$ is composite and $$ b(n,p) = b(n,pp(p)) \text{.} $$

An implementation in Mathematica:

b[n_?NumericQ, 2] := If[EvenQ[n], 1, 0];
b[n_?NumericQ, p_ /; PrimeQ[p]]       := 
    b[n, p] = Sum[ b[n - k p, NextPrime[p, -1] ], {k, 0, Floor[n/p]}]
b[n_?NumericQ, p_ /; Not[PrimeQ[p]] ] := 
    b[n, p] = b[n, NextPrime[p, -1]];
b[n_?NumericQ] := b[n, n]

The idiom b[n_,m_] := b[n,m] = ... memoizes computed values of $b$. It causes Mathematica to store the particular values as short-circuit evaluations benefiting future recursions that require them. That is, $b(11,11) = 6$ is only ever computed by recursion once until the kernel is restarted (or you Clear the assignment b[11,11] = 6). The definition b[n_] := b[n,n] makes a unary version of $b$ that computes the number of partitions of $n$ using primes no larger then $n$.

Table[{n, b[n]}, {n, 2, 20}] // TableForm
(*  2   1
    3   1
    4   1
    5   2
    6   2
    7   3
    8   3
    9   4
    10  5
    11  6
    12  7
    13  9
    14  10
    15  12
    16  14
    17  17
    18  19
    19  23
    20  26  *)

b[100]
(*  40899  *)

b[200]
(*  9845164  *)

For large arguments, we need to prevent the mild attempts to prevent infinite loops.

Block[{
    $RecursionLimit = 5000, 
    $IterationLimit = 5000
  }, 
  DiscretePlot[{Log[n, b[n]]}, {n, 2, 2000}]
]

TableForm[Table[{100 k, b[100 k]}, {k, 1, 30}]]

$$ \begin{array}{rr} n & b(n) \\ \hline 100 & 40\,899 \\ 200 & 9\,845\,164 \\ 300 & 627\,307\,270 \\ 400 & 20\,075\,018\,700 \\ 500 & 414\,270\,104\,287 \\ 600 & 6\,267\,622\,640\,718 \\ 700 & 75\,023\,235\,861\,131 \\ 800 & 746\,579\,898\,675\,387 \\ 900 & 6\,392\,242\,051\,193\,026 \\ 1000 & 48\,278\,613\,741\,845\,757 \\ 1100 & 327\,744\,293\,190\,020\,336 \\ 1200 & 2\,029\,147\,863\,444\,364\,710 \\ 1300 & 11\,590\,550\,282\,322\,997\,472 \\ 1400 & 61\,654\,092\,345\,462\,767\,992 \\ 1500 & 307\,767\,211\,680\,878\,103\,266 \\ 1600 & 1\,450\,997\,588\,181\,557\,017\,431 \\ 1700 & 6\,495\,951\,378\,101\,560\,661\,960 \\ 1800 & 27\,743\,014\,519\,009\,825\,296\,856 \\ 1900 & 113\,481\,312\,640\,973\,435\,202\,574 \\ 2000 & 446\,121\,153\,521\,479\,463\,708\,832 \\ 2100 & 1\,690\,633\,708\,394\,918\,231\,947\,808 \\ 2200 & 6\,192\,537\,346\,397\,477\,180\,809\,944 \\ 2300 & 21\,975\,198\,457\,476\,870\,147\,875\,871 \\ 2400 & 75\,709\,986\,489\,697\,972\,886\,549\,803 \\ 2500 & 253\,714\,665\,648\,596\,332\,053\,234\,626 \\ 2600 & 828\,407\,291\,699\,814\,300\,659\,000\,601 \\ 2700 & 2\,639\,436\,199\,390\,992\,422\,102\,282\,380 \\ 2800 & 8\,217\,667\,977\,417\,902\,397\,061\,965\,136 \\ 2900 & 25\,032\,453\,880\,409\,897\,119\,453\,395\,767 \\ 3000 & 74\,692\,232\,370\,346\,735\,766\,630\,583\,120 \\ \end{array} $$

A semantically equivalent implementation of $b$ in Python (2.7) which is set up to produce the table of values shown above.

from numbers import Number
from sympy import isprime, prevprime

cache = dict()

def b(n, p = None):
    if p == None:
        p = n

    if not (( isinstance(n, Number) ) and (n >= 0)):
        raise TypeError('First argument to function b() must be an integer >= 0.')
    if not (( isinstance(p, Number) ) and (p >= 2)):
        raise TypeError('Second argument to function b() must be an integer >= 2.')

    if (n, p) in cache:
        return cache[(n,p)]

    if p == 2:
        # If n is even, return 1, else return 0.
        retVal = 1 - (n % 2)

    elif isprime(p):
        retVal = sum(( b(n - k * p, prevprime(p) ) for k in range(0, 1 + n//p)) )
    else:
        retVal = b(n, prevprime(p))

    cache[(n, p)] = retVal
    return retVal

print( "b(100) = " + str(b(100)) )
print( "b(200) = " + str(b(200)) )
for k in range(1, 31):
    print( "b(" + str(100*k) + ") = " + str(b(100*k)) )