Partition function P
You are encouraged to solve this task according to the task description, using any language you may know.
The Partition Function P, often notated P(n) is the number of solutions where n∈ℤ can be expressed as the sum of a set of positive integers. For example, P(4)=5 because 4=Σ(4)=Σ(3,1)=Σ(2,2)=Σ(2,1,1)=Σ(1,1,1,1).
P(n) can be expressed as the recurrence relation: P(n) = P(n-1) +P(n-2) -P(n-5) -P(n-7) +P(n-12) +P(n-15) -P(n-22) -P(n-26) +P(n-35) +P(n-40)...
The successive numbers in the above equation have the differences: 1, 3, 2, 5, 3, 7, 4, 9, 5, 11, 6, 13, 7, 15, 8...
This task may be of popular interest because Mathologer made the video, The hardest "What comes next?" (Euler's pentagonal formula), where he asks the programmers among his viewers to calculate P(666). The video has been viewed more than 100,000 times in the first couple of weeks since its release.
In Wolfram Language, this function has been implemented as PartitionsP.
- Task
Write a function which returns the value of PartitionsP(n). Solutions can be iterative or recursive. Bonus task: show how long it takes to compute PartitionsP(6666).
- References
- The hardest "What comes next?" (Euler's pentagonal formula) The explanatory video by Mathologer that makes this task a popular interest.
- Partition Function P Mathworld entry for the Partition function.
- Partition function (number theory) Wikipedia entry for the Partition function.
- Related tasks
C
<lang c>#include <stdint.h>
- include <stdlib.h>
- include <stdio.h>
- include <time.h>
- include <gmp.h>
mpz_t* partition(uint64_t n) { mpz_t *pn = (mpz_t *)malloc((n + 2) * sizeof(mpz_t)); mpz_init_set_ui(pn[0], 1); mpz_init_set_ui(pn[1], 1); for (uint64_t i = 2; i < n + 2; i ++) { mpz_init(pn[i]); for (uint64_t k = 1, penta; ; k++) { penta = k * (3 * k - 1) >> 1; if (penta >= i) break; if (k & 1) mpz_add(pn[i], pn[i], pn[i - penta]); else mpz_sub(pn[i], pn[i], pn[i - penta]); penta += k; if (penta >= i) break; if (k & 1) mpz_add(pn[i], pn[i], pn[i - penta]); else mpz_sub(pn[i], pn[i], pn[i - penta]); } } mpz_t *tmp = &pn[n + 1]; for (uint64_t i = 0; i < n + 1; i ++) mpz_clear(pn[i]); free(pn); return tmp; }
int main(int argc, char const *argv[]) { clock_t start = clock(); mpz_t *p = partition(6666); gmp_printf("%Zd\n", p); printf("Elapsed time: %.04f seconds\n", (double)(clock() - start) / (double)CLOCKS_PER_SEC); return 0; }</lang>
- Output:
193655306161707661080005073394486091998480950338405932486880600467114423441282418165863 Elapsed time: 0.0136 seconds
C++
<lang cpp>#include <chrono>
- include <iostream>
- include <vector>
- include <gmpxx.h>
using big_int = mpz_class;
big_int partitions(int n) {
std::vector<big_int> p(n + 1); p[0] = 1; for (int i = 1; i <= n; ++i) { for (int k = 1, sign = 1;; ++k, sign = -sign) { int m = (k * (3*k - 1))/2; if (m > i) break; if (sign == 1) p[i] += p[i - m]; else p[i] -= p[i - m]; m = (k * (3*k + 1))/2; if (m > i) break; if (sign == 1) p[i] += p[i - m]; else p[i] -= p[i - m]; } } return p[n];
}
int main() {
auto start = std::chrono::steady_clock::now(); auto result = partitions(6666); auto end = std::chrono::steady_clock::now(); std::chrono::duration<double, std::milli> ms(end - start); std::cout << result << '\n'; std::cout << "elapsed time: " << ms.count() << " milliseconds\n";
}</lang>
- Output:
193655306161707661080005073394486091998480950338405932486880600467114423441282418165863 elapsed time: 9.06659 milliseconds
Factor
<lang factor>USING: kernel lists lists.lazy math sequences sequences.extras ;
! Compute the nth pentagonal number.
- penta ( n -- m ) [ sq 3 * ] [ - 2/ ] bi ;
! An infinite lazy list of indices to add and subtract in the ! sequence of partitions to find the next P.
- seq ( -- list )
1 lfrom [ penta 1 - ] <lazy-map> 1 lfrom [ neg penta 1 - ] <lazy-map> lmerge ;
! Reduce a sequence by adding two, subtracting two, adding two, ! etc...
- ++-- ( seq -- n ) 0 [ 2/ odd? [ - ] [ + ] if ] reduce-index ;
! Find the next member of the partitions sequence.
- step ( seq pseq -- seq 'pseq )
dup length [ < ] curry pick swap take-while over <reversed> nths ++-- suffix! ;
- partitions ( m -- n )
[ seq swap [ <= ] curry lwhile list>array ] [ V{ 1 } clone swap [ step ] times last nip ] bi ;</lang>
- Output:
IN: scratchpad [ 6666 partitions ] time . Running time: 0.084955341 seconds 193655306161707661080005073394486091998480950338405932486880600467114423441282418165863
Go
I also tried using Euler's generating function but it was about 20 times slower than this version. <lang go>package main
import (
"fmt" "math/big" "time"
)
var p []*big.Int var pd []int
func partDiffDiff(n int) int {
if n&1 == 1 { return (n + 1) / 2 } return n + 1
}
func partDiff(n int) int {
if n < 2 { return 1 } pd[n] = pd[n-1] + partDiffDiff(n-1) return pd[n]
}
func partitionsP(n int) {
if n < 2 { return } psum := new(big.Int) for i := 1; i <= n; i++ { pdi := partDiff(i) if pdi > n { break } sign := int64(-1) if (i-1)%4 < 2 { sign = 1 } t := new(big.Int).Mul(p[n-pdi], big.NewInt(sign)) psum.Add(psum, t) } p[n] = psum
}
func main() {
start := time.Now() const N = 6666 p = make([]*big.Int, N+1) pd = make([]int, N+1) p[0], pd[0] = big.NewInt(1), 1 p[1], pd[1] = big.NewInt(1), 1 for n := 2; n <= N; n++ { partitionsP(n) } fmt.Printf("p[%d)] = %d\n", N, p[N]) fmt.Printf("Took %s\n", time.Since(start))
}</lang>
- Output:
p[6666)] = 193655306161707661080005073394486091998480950338405932486880600467114423441282418165863 Took 54.82683ms
Julia
Recursive
<lang Julia>using Memoize
function partDiffDiff(n::Int)::Int
isodd(n) ? (n+1)÷2 : n+1
end
@memoize function partDiff(n::Int)::Int
n<2 ? 1 : partDiff(n-1)+partDiffDiff(n-1)
end
@memoize function partitionsP(n::Int)
T=BigInt if n<2 one(T) else psum = zero(T) for i ∈ 1:n pd = partDiff(i) if pd>n break end if ((i-1)%4)<2 psum += partitionsP(n-pd) else psum -= partitionsP(n-pd) end end psum end
end
n=6666 @time println("p($n) = ", partitionsP(n))</lang>
- Output:
p(6666) = 193655306161707661080005073394486091998480950338405932486880600467114423441282418165863 0.260310 seconds (3.58 M allocations: 77.974 MiB, 8.54% gc time)
Phix
Not exactly super-fast, but this sort of stuff is not really what Phix does best. <lang Phix>function partDiffDiff(integer n)
return (n+1)/(1+and_bits(n,1))
end function
sequence pd = {1} function partDiff(integer n)
while n>length(pd) do pd &= pd[$] + partDiffDiff(length(pd)) end while return pd[max(1,n)]
end function
include mpfr.e
sequence pn = {mpz_init(1)} function partitionsP(integer n)
mpz res = mpz_init(1) while n>length(pn) do integer nn = length(pn)+1 mpz psum = mpz_init(0) for i=1 to nn do integer pd = partDiff(i) if pd>nn then exit end if integer sgn = iff(remainder(i-1,4)<2 ? 1 : -1) mpz pnmpd = pn[max(1,nn-pd)] if sgn=-1 then mpz_sub(psum,psum,pnmpd) else mpz_add(psum,psum,pnmpd) end if end for pn = append(pn,psum) end while return pn[max(1,n)]
end function
atom t0 = time() integer n=6666 printf(1,"p(%d) = %s (%s)\n",{n,mpz_get_str(partitionsP(n)),elapsed(time()-t0)})</lang>
- Output:
p(6666) = 193655306161707661080005073394486091998480950338405932486880600467114423441282418165863 (0.8s)
Python
This follows the algorithm from the Mathloger video closely
<lang python>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)])</lang>
- Output:
(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]
- Stretch goal
From command line after running the above
In [2]: %time part(6666) Wall time: 243 ms Out[2]: 193655306161707661080005073394486091998480950338405932486880600467114423441282418165863
Python: Mathloger video prime generator
Add the following code after that above <lang python>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)))</lang>
- Output:
Primes: [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]
Raku
Not the fastest, but it gets the job done. <lang perl6>my @P = 1, { p(++$) } … *; my @i = lazy [\+] flat 1, ( (1 .. *) Z (1 .. *).map: * × 2 + 1 ); sub p ($n) { sum @P[$n X- @i[^(@i.first: * > $n, :k)]] Z× (flat (1, 1, -1, -1) xx *) }
put @P[^26]; put @P[6666];</lang>
- Output:
1 1 2 3 5 7 11 15 22 30 42 56 77 101 135 176 231 297 385 490 627 792 1002 1255 1575 1958 193655306161707661080005073394486091998480950338405932486880600467114423441282418165863
REXX
This REXX version is recursive. <lang rexx>/*REXX program calculates and displays a specific value (or a range of) partitionsP(N).*/ numeric digits 1000 /*able to handle some ginormous numbers*/ parse arg lo hi . /*obtain optional arguments from the CL*/ if lo== | lo=="," then lo= 0 /*Not specified? Then use the default.*/ if hi== | hi=="," then hi= lo /* " " " " " " */ @.= 0; @.0= 1; @.1= 1; @.2= 2; @.3= 3; @.4= 5 /*assign default value and low values. */ w= length( commas(hi) ) /*W: is used for aligning the index. */
do j=lo to hi /*compute a range of partitionsP. */ say right( commas(j), w) ' ' commas( partP(j) ) end /*j*/
exit 0 /*stick a fork in it, we're all done. */ /*──────────────────────────────────────────────────────────────────────────────────────*/ commas: parse arg ?; do jc=length(?)-3 to 1 by -3; ?=insert(',', ?, jc); end; return ? /*──────────────────────────────────────────────────────────────────────────────────────*/ partP: procedure expose @.; parse arg n /*obtain number (index) for computation*/
if @.n\==0 then return @.n /*Is it already computed? Return it. */ #= 0 /*initialize part P number.*/ do k=1 for n; z= n - (k*3 - 1) * k % 2 /*compute the partition P num*/ if z<0 then leave /*Is Z negative? Then leave.*/ if @.z==0 then x= partP(z) /*use recursion if not known.*/ else x= @.z /*use the pre─computed number*/ z= z - k /*subtract index (K) from Z. */ if z<0 then y= 0 /*Is Z negative? Then set Y=0*/ else if @.z==0 then y= partP(z) /*use recursion if not known.*/ else y= @.z /*use the pre─computed number*/ if k//2 then #= # + x + y /*Odd? Then sum X and Y.*/ else #= # - x - y /*Even? " subtract " " ".*/ end /*k*/ @.n= #; return # /*define and return partitionsP of N. */</lang>
- output when using the input of: 6666
6,666 193,655,306,161,707,661,080,005,073,394,486,091,998,480,950,338,405,932,486,880,600,467,114,423,441,282,418,165,863
- output when using the input of: 66666
66,666 931,283,431,869,095,717,238,416,063,534,148,471,363,928,685,651,267,074,563,360,050,977,549,251,436,058,076,515,262,033,789,845,457,356,081,278,451,246,751,375,974,038,318,319,810,923,102,769,109,446,980,055,567,090,089,060,954,748,065,131,666,952,830,623,286,286,824,837,240,058,805,177,319,865,230,913,417,587,625,830,803,675,380,262,765,598,632,742,903,192,085,254,824,621
Rust
<lang rust>// [dependencies] // rug = "1.11"
use rug::Integer;
fn partitions(n: usize) -> Integer {
let mut p = Vec::with_capacity(n + 1); p.push(Integer::from(1)); for i in 1..=n { let mut num = Integer::from(0); let mut k = 1; let mut sign = 1; loop { let mut j = (k * (3 * k - 1))/2; if j > i { break; } if sign == 1 { num += &p[i - j]; } else { num -= &p[i - j]; } j += k; if j > i { break; } if sign == 1 { num += &p[i - j]; } else { num -= &p[i - j]; } k += 1; sign = -sign; } p.push(num); } p[n].clone()
}
fn main() {
use std::time::Instant; let n = 6666; let now = Instant::now(); let result = partitions(n); let time = now.elapsed(); println!("P({}) = {}", n, result); println!("elapsed time: {} microseconds", time.as_micros());
}</lang>
- Output:
P(6666) = 193655306161707661080005073394486091998480950338405932486880600467114423441282418165863 elapsed time: 9291 microseconds
Wren
Although it may not look like it, this is actually a decent time for Wren which is interpreted and the above module is written entirely in Wren itself. <lang ecmascript>import "/big" for BigInt
var p = [] var pd = []
var partDiffDiff = Fn.new { |n| (n&1 == 1) ? (n + 1)/2 : n + 1 }
var partDiff = Fn.new { |n|
if (n < 2) return 1 pd[n] = pd[n-1] + partDiffDiff.call(n-1) return pd[n]
}
var partitionsP = Fn.new { |n|
if (n < 2) return var psum = BigInt.zero for (i in 1..n) { var pdi = partDiff.call(i) if (pdi > n) break var sign = (i-1)%4 < 2 ? 1 : -1 psum = psum + p[n-pdi] * sign } p[n] = psum
}
var start = System.clock var N = 6666 p = List.filled(N+1, null) pd = List.filled(N+1, 0) p[0] = BigInt.one p[1] = BigInt.one pd[0] = 1 pd[1] = 1 for (n in 2..N) partitionsP.call(n) System.print("p[%(N)] = %(p[N])") System.print("Took %(System.clock - start) seconds")</lang>
- Output:
p[6666] = 193655306161707661080005073394486091998480950338405932486880600467114423441282418165863 Took 1.428762 seconds