I watched a great Mathologer video on youtube which examined the number of integer partitions of natural numbers. Here's the video:

.A partition of a natural number is number of ways that positive integers can be summed to equal the number, so for 4:

4 = 4 4 = 1 + 3 4 = 2 + 2 4 = 1 + 1 + 2 4 = 3 + 1 4 = 1 + 2 + 1 4 = 2 + 1 + 1 4 = 1 + 1 + 1 + 1

Counting shows there are 8 "partitions". Normally duplicates when disregarding order are removed so, for example, 3 + 1 is a duplicate of 1 + 3. These are the actual partitions that the video then concentrates on. E.g again, for four the actual partitions are:

4 = 4 4 = 2 + 2 4 = 1 + 3 4 = 1 + 1 + 2 4 = 1 + 1 + 1 + 1

giving a count of the partitions as 5.

## Partition finder.

The video shows that by considering the n-1 "gaps" between n blocks as being either open or closed you can generate all the actual partitions of n things (with duplicates).

This lead to me writing the following code:

def part_by(count, n): " Opening/closing n - 1 gaps between n blocks according to binary count" p, num = [], 1 for _ in range(n - 1): if count & 1: num += 1 else: p.append(num) num = 1 count //= 2 p.append(num) return p def partitions_with_dups(n): "All ways of opening/closing n - 1 gaps between n blocks" return [part_by(i, n) for i in range(2**(n - 1) if n else 1)] def partitions(n): "partitions removing duplications in order of terms" return sorted(set(tuple(sorted(part_by(i, n))) for i in range(2**(n - 1) if n else 1))) def partition(n): "Count of different partitions irrespective of order" return len(partitions(n)) if __name__ == '__main__': print("partition(n) for n from 0... = ", [partition(n) for n in range(11)])

Running the code produces the correct sequence of partiton(n)

partition(n) for n from 0... = [1, 1, 2, 3, 5, 7, 11, 15, 22, 30, 42]

## Eulers reduction/formula

The video in chapter 2 tries a modified fibonacci-like algorithm to generate the sequence. I coded the algorithm likes this, and submitted it to a rosetta Code task written by someone else who was inspired by the video.

from itertools import islice def posd(): "diff between position numbers. 1, 2, 3... interleaved with 3, 5, 7..." count, odd = 1, 3 while True: yield count yield odd count, odd = count + 1, odd + 2 def pos_gen(): "position numbers. 1 3 2 5 7 4 9 ..." val = 1 diff = posd() while True: yield val val += next(diff) def plus_minus(): "yield (list_offset, sign) or zero for Partition calc" n, sign = 0, [1, 1] p_gen = pos_gen() out_on = next(p_gen) while True: n += 1 if n == out_on: next_sign = sign.pop(0) if not sign: sign = [-next_sign] * 2 yield -n, next_sign out_on = next(p_gen) else: yield 0 def part(n): "Partition numbers" p = [1] p_m = plus_minus() mods = [] for _ in range(n): next_plus_minus = next(p_m) if next_plus_minus: mods.append(next_plus_minus) p.append(sum(p[offset] * sign for offset, sign in mods)) return p[-1] print("(Intermediaries):") print(" posd:", list(islice(posd(), 10))) print(" pos_gen:", list(islice(pos_gen(), 10))) print(" plus_minus:", list(islice(plus_minus(), 15))) print("\nPartitions:", [part(x) for x in range(15)])

The narrator goes on to show how the slightly modified series can be used to generate prime numbers, which I coded as the following addition to the last chunk of code:

def par_primes(): "Prime number generator from the partition machine" p = [1] p_m = plus_minus() mods = [] n = 0 while True: n += 1 next_plus_minus = next(p_m) if next_plus_minus: mods.append(next_plus_minus) p.append(sum(p[offset] * sign for offset, sign in mods)) if p[0] + 1 == p[-1]: yield p[0] p[0] += 1 yield p print("\nPrimes:", list(islice(par_primes(), 15)))

The output from the previous two code sections is this:

(Intermediaries): posd: [1, 3, 2, 5, 3, 7, 4, 9, 5, 11] pos_gen: [1, 2, 5, 7, 12, 15, 22, 26, 35, 40] plus_minus: [(-1, 1), (-2, 1), 0, 0, (-5, -1), 0, (-7, -1), 0, 0, 0, 0, (-12, 1), 0, 0, (-15, 1)] Partitions: [1, 1, 2, 3, 5, 7, 11, 15, 22, 30, 42, 56, 77, 101, 135] Primes: [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]

This Euler algorithm is very fast as the old one needed generating 2**(n - 1) actual partitions.

## Other (recursive) method.

I had originally thought to solve the Rosetta Code task by searching for an easier algorithm and implementing that.

I found this stackoverflow answer which gave a recursive formula that seemed more straight-forward.

I did find one issue with the explanation but coded the algorithm up as follows:

from functools import lru_cache import sys sys.setrecursionlimit(10_000_000) @lru_cache(maxsize=None) def _p(n, m): if n < 0: return 0 elif n == 0: return 1 elif n / 2 < m < n: return 1 else: return sum(_p(n - i, i) for i in range(m, n+1)) def p(n): return _p(n, 1) #%% for n in [666] + list(range(1000, 6667, 50)) + [6666]: print(n, p(n))

**
Ooh such short and succinct code!**

You might be able to tell though that it too has efficiency issues. It will calculate p(666) but the Rosetta Code task asks for the calculation of p(6666) - with four sixes not three. If I just ask for p(6666) then Python crashes, silently!

If you I gradually jump up to p(6666) then it will eventually give the answer of 193655306161707661080005073394486091998480950338405932486880600467114423441282418165863 I think the cache needs to grow to limit recursion depth. This routine took hours to get to p(6666) where Eulers algorithm took milliseconds.

Note: The following cache information after calculating p(6666)

>>> _p.cache_info() CacheInfo(hits=21_921_546_630, misses=22_221_112, maxsize=None, currsize=22_221_112) >>>

**END.**

## No comments:

## Post a Comment