I'm working on modernizing Rosetta Code's infrastructure. Starting with communications. Please accept this time-limited open invite to RC's Slack.. --Michael Mol (talk) 20:59, 30 May 2020 (UTC)

Partition function P

From Rosetta Code
Task
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.


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


Related tasks




C[edit]

Library: GMP
#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;
}
Output:
193655306161707661080005073394486091998480950338405932486880600467114423441282418165863
Elapsed time: 0.0136 seconds


C++[edit]

Library: GMP
#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;; ++k) {
int j = (k * (3*k - 1))/2;
if (j > i)
break;
if (k & 1)
p[i] += p[i - j];
else
p[i] -= p[i - j];
j = (k * (3*k + 1))/2;
if (j > i)
break;
if (k & 1)
p[i] += p[i - j];
else
p[i] -= p[i - j];
}
}
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";
}
Output:
193655306161707661080005073394486091998480950338405932486880600467114423441282418165863
elapsed time: 8.99497 milliseconds

Delphi[edit]

Translation of: Go
 
program Partition_function_P;
 
{$APPTYPE CONSOLE}
 
uses
System.SysUtils,
Velthuis.BigIntegers,
System.Diagnostics;
 
var
p: TArray<BigInteger>;
pd: TArray<Integer>;
 
function PartiDiffDiff(n: Integer): Integer;
begin
if n and 1 = 1 then
exit((n + 1) div 2);
Result := n + 1;
end;
 
function partDiff(n: Integer): Integer;
begin
if n < 2 then
exit(1);
 
pd[n] := pd[n - 1] + PartiDiffDiff(n - 1);
Result := pd[n];
end;
 
procedure partitionP(n: Integer);
begin
if n < 2 then
exit;
 
var psum: BigInteger := 0;
for var i := 1 to n do
begin
var pdi := partDiff(i);
if pdi > n then
Break;
 
var sign: Int64 := -1;
 
if (i - 1) mod 4 < 2 then
sign := 1;
 
var t: BigInteger := BigInteger(p[n - pdi]) * BigInteger(sign);
psum := psum + t;
end;
p[n] := psum;
end;
 
begin
var stopwatch := TStopwatch.Create;
const n = 6666;
SetLength(p, n + 1);
SetLength(pd, n + 1);
stopwatch.Start;
p[0] := 1;
pd[0] := 1;
p[1] := 1;
pd[1] := 1;
for var i := 2 to n do
partitionP(i);
stopwatch.Stop;
writeln(format('p[%d] = %s', [n, p[n].ToString]));
writeln('Took ', stopwatch.ElapsedMilliseconds, 'ms');
Readln;
end.
Output:
p[6666] = 193655306161707661080005073394486091998480950338405932486880600467114423441282418165863
Took 131ms

F#[edit]

An implementation of the formula in the task description. P(123456) is included for comparison with the largest value in the related task.

 
// PartionP: Nigel Galloway. April 12th., 2021
let pP g=let rec fN i g e l=seq{yield(l,e+i);yield(-l,e+i+g);yield! fN(i+1)(g+2)(e+i+g)(-l)}
let N,G=Array.create(g+1) 1I,seq{yield (1I,1);yield! fN 1 3 1 1I}|>Seq.takeWhile(fun(_,n)->n<=g)|>List.ofSeq
seq{2..g}|>Seq.iter(fun p->N.[p]<-G|>List.takeWhile(fun(_,n)->n<=p)|>Seq.fold(fun Σ (n,g)->Σ+n*N.[p-g]) 0I); N.[g]
printfn "666->%A\n\n6666->%A\n\n123456->%A" (pP 666) pP(6666) (pP 123456)
 
Output:
666->11956824258286445517629485

6666->193655306161707661080005073394486091998480950338405932486880600467114423441282418165863
Real: 00:00:00.096

123456->30817659578536496678545317146533980855296613274507139217608776782063054452191537379312358383342446230621170608408020911309259407611257151683372221925128388387168451943800027128045369650890220060901494540459081545445020808726917371699102825508039173543836338081612528477859613355349851184591540231790254269948278726548570660145691076819912972162262902150886818986555127204165221706149989

Factor[edit]

Works with: Factor version 0.99 2020-08-14
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 ;
Output:
IN: scratchpad [ 6666 partitions ] time .

Running time: 0.084955341 seconds

193655306161707661080005073394486091998480950338405932486880600467114423441282418165863

Go[edit]

Translation of: Julia

I also tried using Euler's generating function but it was about 20 times slower than this version.

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))
}
Output:
p[6666)] = 193655306161707661080005073394486091998480950338405932486880600467114423441282418165863
Took 54.82683ms

Java[edit]

import java.math.BigInteger;
 
public class PartitionFunction {
public static void main(String[] args) {
long start = System.currentTimeMillis();
BigInteger result = partitions(6666);
long end = System.currentTimeMillis();
System.out.println("P(6666) = " + result);
System.out.printf("elapsed time: %d milliseconds\n", end - start);
}
 
private static BigInteger partitions(int n) {
BigInteger[] p = new BigInteger[n + 1];
p[0] = BigInteger.ONE;
for (int i = 1; i <= n; ++i) {
p[i] = BigInteger.ZERO;
for (int k = 1; ; ++k) {
int j = (k * (3 * k - 1))/2;
if (j > i)
break;
if ((k & 1) != 0)
p[i] = p[i].add(p[i - j]);
else
p[i] = p[i].subtract(p[i - j]);
j += k;
if (j > i)
break;
if ((k & 1) != 0)
p[i] = p[i].add(p[i - j]);
else
p[i] = p[i].subtract(p[i - j]);
}
}
return p[n];
}
}
Output:
P(6666) = 193655306161707661080005073394486091998480950338405932486880600467114423441282418165863
elapsed time: 59 milliseconds

Julia[edit]

Recursive[edit]

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))
Output:
p(6666) = 193655306161707661080005073394486091998480950338405932486880600467114423441282418165863
  0.260310 seconds (3.58 M allocations: 77.974 MiB, 8.54% gc time)


Maple[edit]

p:=proc(n)
option remember;
local k,s:=0,m;
for k from 1 while (m:=iquo(k*(3*k-1),2))<=n do
s-=(-1)^k*p(n-m);
od;
for k from 1 while (m:=iquo(k*(3*k+1),2))<=n do
s-=(-1)^k*p(n-m);
od;
s
end:
p(0):=1:
 
time(p(6666));
# 0.796
 
time(combinat[numbpart](6666));
# 0.406
 
p~([$1..20]);
# [1, 2, 3, 5, 7, 11, 15, 22, 30, 42, 56, 77, 101, 135, 176, 231, 297, 385, 490, 627]
 
combinat[numbpart]~([$1..20]);
# [1, 2, 3, 5, 7, 11, 15, 22, 30, 42, 56, 77, 101, 135, 176, 231, 297, 385, 490, 627]
 
p(1000)
# 24061467864032622473692149727991
 
combinat[numbpart](1000);
# 24061467864032622473692149727991

Nim[edit]

Translation of: C++
Library: bignum
import sequtils, strformat, times
import bignum
 
func partitions(n: int): Int =
var p = newSeqWith(n + 1, newInt())
p[0] = newInt(1)
for i in 1..n:
var k = 1
while true:
var j = k * (3 * k - 1) div 2
if j > i: break
if (k and 1) != 0:
inc p[i], p[i - j]
else:
dec p[i], p[i - j]
j = k * (3 * k + 1) div 2
if j > i: break
if (k and 1) != 0:
inc p[i], p[i - j]
else:
dec p[i], p[i - j]
inc k
result = p[n]
 
let t0 = cpuTime()
echo partitions(6666)
echo &"Elapsed time: {(cpuTime() - t0) * 1000:.2f} ms"
Output:
193655306161707661080005073394486091998480950338405932486880600467114423441282418165863
Elapsed time: 9.73 ms

Perl[edit]

use strict;
use warnings;
no warnings qw(recursion);
use Math::AnyNum qw(:overload);
use Memoize;
 
memoize('partitionsP');
memoize('partDiff');
 
sub partDiffDiff { my($n) = @_; $n%2 != 0 ? ($n+1)/2 : $n+1 }
 
sub partDiff { my($n) = @_; $n<2 ? 1 : partDiff($n-1) + partDiffDiff($n-1) }
 
sub partitionsP {
my($n) = @_;
return 1 if $n < 2;
 
my $psum = 0;
for my $i (1..$n) {
my $pd = partDiff($i);
last if $pd > $n;
if ( (($i-1)%4) < 2 ) { $psum += partitionsP($n-$pd) }
else { $psum -= partitionsP($n-$pd) }
}
return $psum
}
 
print partitionsP($_) . ' ' for 0..25; print "\n";
print partitionsP(6666) . "\n";
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

Phix[edit]

Not exactly super-fast, but this sort of stuff is not really what Phix does best.

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)})
Output:
p(6666) = 193655306161707661080005073394486091998480950338405932486880600467114423441282418165863 (0.8s)


Picat[edit]

 
/* Picat 3.0#5 */
/* Author: Hakan Kjellerstrand */
table
partition1(0) = 1.
partition1(N) = P =>
S = 0,
K = 1,
M = (K*(3*K-1)) // 2,
while (M <= N)
S := S - ((-1)**K)*partition1(N-M),
K := K + 1,
M := (K*(3*K-1)) // 2
end,
K := 1,
M := (K*(3*K+1)) // 2,
while (M <= N)
S := S - ((-1)**K)*partition1(N-M),
K := K + 1,
M := (K*(3*K+1)) // 2
end,
P = S.
 
Picat> time(println('p(6666)'=partition1(6666)))
p(6666) = 193655306161707661080005073394486091998480950338405932486880600467114423441282418165863
 
CPU time 0.206 seconds.
 

Prolog[edit]

 
/* SWI-Prolog 8.3.21 */
/* Author: Jan Burse */
:- table p/2.
p(0, 1) :- !.
p(N, X) :-
aggregate_all(sum(Z), (between(1,inf,K), M is K*(3*K-1)//2,
(M>N, !, fail; L is N-M, p(L,Y), Z is (-1)^K*Y)), A),
aggregate_all(sum(Z), (between(1,inf,K), M is K*(3*K+1)//2,
(M>N, !, fail; L is N-M, p(L,Y), Z is (-1)^K*Y)), B),
X is -A-B.
 
?- time(p(6666,X)).
% 13,962,294 inferences, 2.610 CPU in 2.743 seconds (95% CPU, 5350059 Lips)
X = 1936553061617076610800050733944860919984809503384
05932486880600467114423441282418165863.
 

Python[edit]

Python: Mathloger[edit]

This follows the algorithm from the Mathloger video closely

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)])
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 [1]: part(6666)
Out[1]: 193655306161707661080005073394486091998480950338405932486880600467114423441282418165863

In [2]: %timeit part(6666)
103 ms ± 477 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

Python: Mathloger video prime generator[edit]

Add the following code after that above

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)))
Output:
Primes: [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]

Python: Alternative[edit]

Translation of: C++
from typing import List
 
 
def partitions(n: int) -> int:
"""Count partitions."""
p: List[int] = [1] + [0] * n
for i in range(1, n + 1):
k: int = 0
while True:
k += 1
j: int = (k * (3*k - 1)) // 2
if (j > i):
break
if (k & 1):
p[i] += p[i - j]
else:
p[i] -= p[i - j]
j = (k * (3*k + 1)) // 2
if (j > i):
break
if (k & 1):
p[i] += p[i - j]
else:
p[i] -= p[i - j]
 
return p[n]
 
 
if __name__ == '__main__':
print("\nPartitions:", [partitions(x) for x in range(15)])
Output:
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 [3]: partitions(6666)
Out[3]: 193655306161707661080005073394486091998480950338405932486880600467114423441282418165863

In [4]: %timeit partitions(6666)
215 ms ± 1.84 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Raku[edit]

Works with: Rakudo version 2020.09

Not the fastest, but it gets the job done.

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];
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[edit]

This REXX version is recursive.

/*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. */
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[edit]

// [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;
loop {
let mut j = (k * (3 * k - 1)) / 2;
if j > i {
break;
}
if (k & 1) == 1 {
num += &p[i - j];
} else {
num -= &p[i - j];
}
j += k;
if j > i {
break;
}
if (k & 1) == 1 {
num += &p[i - j];
} else {
num -= &p[i - j];
}
k += 1;
}
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());
}
Output:
P(6666) = 193655306161707661080005073394486091998480950338405932486880600467114423441282418165863
elapsed time: 8912 microseconds

Sidef[edit]

Built-in:

say partitions(6666)   # very fast

User-defined:

func partitionsP(n) {
func (n) is cached {
 
n <= 1 && return n
 
var a = sum(1..floor((sqrt(24*n + 1) + 1)/6), {|k|
(-1)**(k-1) * __FUNC__(n - ((k*(3*k - 1)) >> 1))
})
 
var b = sum(1..ceil((sqrt(24*n + 1) - 7)/6), {|k|
(-1)**(k-1) * __FUNC__(n - ((k*(3*k + 1)) >> 1))
})
 
a + b
}(n+1)
}
 
var t = Time.micro
 
say partitionsP.map(0..25).join(' ')
say partitionsP(6666)
 
say ("Took %.4f seconds" % Time.micro-t)
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
Took 24.5225 seconds

Wren[edit]

Translation of: Julia
Library: Wren-big

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.

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")
Output:
p[6666] = 193655306161707661080005073394486091998480950338405932486880600467114423441282418165863
Took 1.428762 seconds