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

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.

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

## 11l

Translation of: Python: Alternative
`F partitions(n)   V p = [BigInt(1)] [+] [BigInt(0)] * n   L(i) 1 .. n      V k = 0      L         k++         V j = (k * (3 * k - 1)) I/ 2         I j > i            L.break         I k [&] 1            p[i] += p[i - j]         E            p[i] -= p[i - j]         j = (k * (3 * k + 1)) I/ 2         I j > i            L.break         I k [&] 1            p[i] += p[i - j]         E            p[i] -= p[i - j]   R p[n] print(‘Partitions: ’(0.<15).map(x -> partitions(x))) V start = time:perf_counter()print(partitions(6666))print(time:perf_counter() - start)`
Output:
```Partitions: [1, 1, 2, 3, 5, 7, 11, 15, 22, 30, 42, 56, 77, 101, 135]
193655306161707661080005073394486091998480950338405932486880600467114423441282418165863
0.598528
```

## C

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++

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

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```

## Elixir

Loosely based on the Erlang version.

` use Bitwise, skip_operators: true defmodule Partition do   def init(), do:      :ets.new :pN, [:set, :named_table, :private]    def gpentagonals(), do: Stream.unfold {1, 0}, &next/1   defp next({m, n}) do      a = case rem m, 2 do             0 -> div m, 2             1 -> m          end      {n, {m + 1, n + a}}   end    def p(0), do: 1   def p(n) do      case :ets.lookup :pN, n do         [{^n, val}] -> val         [] ->            {val, _} = gpentagonals()                       |> Stream.drop(1)                       |> Stream.take_while(fn m -> m <= n end)                       |> Stream.map(fn g -> p(n - g) end)                       |> Enum.reduce({0, 0},                              fn n, {a, sgn} -> {                                          a + (if sgn < 2, do: n, else: -n),                                          band(sgn + 1, 3)                                          }                              end)            :ets.insert :pN, {n, val}            val      end   endend Partition.initIO.puts Partition.p 6666 `
Output:
```\$ time ./partition6666.ex
193655306161707661080005073394486091998480950338405932486880600467114423441282418165863

real	0m1.106s
user	0m1.191s
sys	0m0.116s
```

## Erlang

` -mode(compile). main(_) ->    ets:new(pN, [set, named_table, protected]),    io:format("~w~n", [p(6666)]). p(0) -> 1;p(N) ->    case ets:lookup(pN, N) of        [{N, Pn}] -> Pn;        [] ->            Terms = [p(N - G) || G <- gpentagonals(N)],            Pn = sum_partitions(Terms),            ets:insert(pN, {N, Pn}),            Pn    end. sum_partitions(Terms) -> sum_partitions(Terms, 0, 0).sum_partitions([], _, Sum) -> Sum;sum_partitions([N|Ns], Sgn, Sum) ->    Summand = case Sgn < 2 of        true  -> N;        false -> -N    end,    sum_partitions(Ns, (Sgn+1) band 3, Sum + Summand). gpentagonals(Max) -> gpentagonals(1, Max, [0]).gpentagonals(M, Max, Ps = [N|_]) ->    GP = N + case M rem 2 of                0 -> M div 2;                1 -> M             end,    if        GP > Max -> tl(lists:reverse(Ps));        true -> gpentagonals(M + 1, Max, [GP|Ps])    end. `
Output:
```\$ time ./partition6666.erl
193655306161707661080005073394486091998480950338405932486880600467114423441282418165863

real	0m0.480s
user	0m0.490s
sys	0m0.080s
```

## F#

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., 2021let 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

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
```

## FreeBASIC

### Unsiged 64bit version

Translation of: Python
`Function PartitionsP(n As UInteger) As ULongInt    ' if n > 416, the result becomes to large for a unsigned 64bit integer    Dim As ULongInt p(n)    Dim As UInteger k, j     p(0) = 1    For i As UInteger = 1 To n        k = 0        While TRUE            k += 1            j = (k * (3*k - 1)) \ 2            If (j > i) Then Exit While            If (k And 1) Then                p(i) += p(i - j)            Else                p(i) -= p(i - j)            End If            'j = (k * (3*k + 1)) \ 2            j += k            If (j > i) Then Exit While            If (k And 1) Then                p(i) += p(i - j)            Else                p(i) -= p(i - j)            End If        Wend    Next i    Return p(n)End Function Print !"\nPartitionsP: ";For x As UInteger = 0 To 12    Print PartitionsP(x);"  ";Next x Print !"\n\ndone"Sleep`
Output:
`PartitionsP: 1  1  2  3  5  7  11  15  22  30  42  56  77`

### Big numbers version

Library: GMP
`' version 26-06-2021' compile with: fbc -s console #Include Once "gmp.bi" Sub PartitionsP(max As ULong, p() As MpZ_ptr)    ' based on Numericana code example    Dim As ULong a, b, i, k    Dim As Long j     Dim As Mpz_ptr s = Allocate(Len(__mpz_struct)) : Mpz_init(s)     Mpz_set_ui(p(0), 1)     For i = 1 To max        j = 1 : k = 1 : b = 2 : a = 5        While j > 0            ' j = i - (3*k*k+k) \ 2            j = i - b : b = b + a : a = a + 3            If j >= 0 Then                If k And 1 Then Mpz_add(s, s, p(j)) Else Mpz_sub(s, s, p(j))            End If            j = j + k            If j >= 0 Then                If k And 1 Then Mpz_add(s, s, p(j)) Else Mpz_sub(s, s, p(j))            End If            k = k +1        Wend        Mpz_swap(p(i), s)    Next     Mpz_clear(s) End Sub ' ------=< MAIN >=------ #Define max 6666 Dim As UInteger nDim As ZString Ptr ansDim As Double t = Timer ReDim big_p(max) As Mpz_ptrFor n = 0 To max    big_p(n) = Allocate(Len(__mpz_struct)) : Mpz_init(big_p(n))Next PartitionsP(max, big_p())ans = Mpz_get_str (0, 10, big_p(max))Print "PartitionsP("; Str(max); ") = "; "  "; *ans For n = 0 To max    Mpz_clear(big_p(n))Next Print Using "time = ###.## ms"; (Timer - t) * 1000 ' empty keyboard bufferWhile InKey <> "" : WendPrint : Print "hit any key to end program"SleepEnd`
Output:
```PartitionsP(6666) =   193655306161707661080005073394486091998480950338405932486880600467114423441282418165863
time =  32.97 ms```

## Go

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.Intvar 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

`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
```

## JavaScript

` function p(n){    var a = new Array(n+1)    a[0] = 1n     for (let i = 1; i <= n; i++){        a[i] = 0n        for (let k = 1, s = 1; s <= i;){            a[i] += (k & 1 ? a[i-s]:-a[i-s])            k > 0 ? (s += k, k = -k):(k = -k+1, s = k*(3*k-1)/2)        }    }     return a[n]} var t = Date.now()console.log("p(6666) = " + p(6666))console.log("Computation time in ms: ", Date.now() - t) `
Output:
```p(6666) = 193655306161707661080005073394486091998480950338405932486880600467114423441282418165863
Computation time in ms: 26
```

`{-# language DeriveFunctor #-} -------------------------------------------------------------- memoization utilities data Memo a = Node a (Memo a) (Memo a)  deriving Functor  memo :: Integral a => Memo p -> a -> pmemo (Node a l r) n  | n == 0 = a  | odd n = memo l (n `div` 2)  | otherwise = memo r (n `div` 2 - 1) nats :: Memo Intnats = Node 0 ((+1).(*2) <\$> nats) ((*2).(+1) <\$> nats) -------------------------------------------------------------- calculating partitions partitions :: Memo Integerpartitions = partitionP <\$> nats partitionP :: Int -> IntegerpartitionP n  | n < 2 = 1  | otherwise = sum \$ zipWith (*) signs terms  where    terms = [ memo partitions (n - i)            | i <- takeWhile (<= n) ofsets ]    signs = cycle [1,1,-1,-1] ofsets = scanl1 (+) \$ mix [1,3..] [1,2..]  where    mix a b = concat \$ zipWith (\x y -> [x,y]) a b  main = print \$ partitionP 6666`
```*Main> partitionP <\$> [1..10]
[1,2,3,5,7,11,15,22,30,42]

*Main> partitionP 100
190569292

*Main> partitionP 1000
24061467864032622473692149727991

*Main> partitionP 6666
193655306161707661080005073394486091998480950338405932486880600467114423441282418165863```

## J

Solution stolen verbatim from the J Wiki. Note the use of memoization (M.) for efficiency:

`   pn =: -/@(+/)@:(\$:"0)@rec ` (x:@(0&=)) @. (0>:]) M.   rec=: - (-: (*"1) _1 1 +/ 3 * ]) @ (>:@[email protected]>[email protected]%:@((2%3)&*))`
Output:
```   pn 6666
193655306161707661080005073394486091998480950338405932486880600467114423441282418165863
```

## jq

Translation of: Python:Alternative

`def partitions(\$n):  def div2: (. - (.%2)) / 2;  reduce range(1; \$n + 1) as \$i ( {p: ([1] + [range(0;\$n)|0])};    . + {k: 0, stop: false}    | until(.stop;        .k += 1        | (((.k * (3*.k - 1)) | div2) ) as \$j        | if \$j > \$i then .stop=true  	  else if (.k % 2) == 1  	       then .p[\$i] = .p[\$i] + .p[\$i - \$j]               else .p[\$i] = .p[\$i] - .p[\$i - \$j]               end          | (((.k * (3*.k + 1)) | div2)) as \$j          | if \$j > \$i then .stop=true            elif (.k % 2) == 1            then .p[\$i] = .p[\$i] + .p[\$i - \$j]            else .p[\$i] = .p[\$i] - .p[\$i - \$j]    	  end  	end ))    | .p[\$n] ; [partitions(range(1;15))]`
Output:
`[1,2,3,5,7,11,15,22,30,42,56,77,101,135]`

jq's built-in integer precision is insufficient for computing ``partitions(6666)``, but more as a test of the BigInt.jq library for jq than anything else, here are the results of using it in conjunction with a trivially-modified version of the partitions implementation above. That is, after modifying the lines that refer to "p" (or ".p"), we see that partitions(6666) yields:

```   "193655306161707661080005073394486091998480950338405932486880600467114423441282418165863"
```

The user+sys time is 7m3s as the BigInt.jq library is written in jq.

## Recursive

Translation of: Julia
with memoization
`def partDiffDiff(\$n):  if (\$n % 2) == 1 then (\$n+1) / 2 else \$n+1 end; # in:  {n, partDiffMemo} # out: object with possibly updated memoizationdef partDiff:  .n as \$n  | if .partDiffMemo[\$n] then .    elif \$n<2 then .partDiffMemo[\$n]=1    else ((.n=(\$n-1)) | partDiff)    | .partDiffMemo[\$n] = .partDiffMemo[\$n-1] + partDiffDiff(\$n-1)    end; # in:  {n, memo, partDiffMemo}#      where `.memo[i]` memoizes partitions(i)#      and   `.partDiffMemo[i]` memoizes partDiff(i) # out: object with possibly updated memoizationdef partitionsM:  .n as \$n  | if .memo[\$n] then .    elif \$n<2 then .memo[\$n] = 1    else label \$out    | foreach range(1; \$n+2) as \$i (.emit = false | .psum = 0;        if \$i > \$n then .emit = true        else ((.n = \$i) | partDiff)	| .partDiffMemo[\$i] as \$pd        | if \$pd > \$n then .emit=true, break \$out          else {psum, emit} as \$local  # for restoring relevant state	  |  ((.n = (\$n-\$pd)) | partitionsM)	  | .memo[\$n-\$pd] as \$increment	  | . + \$local                 # restore          | if ((\$i-1)%4)<2            then .psum += \$increment            else .psum -= \$increment            end	  end        end;        select(.emit) )    | .memo[\$n] = .psum    end ; def partitionsP:  . as \$n  | {n: \$n, memo:[], partDiffMemo:[]}  | partitionsM  | .memo[\$n]; # Stretch goal:6666 | partitionsP `
Using gojq, the above program takes 41.35 seconds (u+s) on a 3GHz Mac Mini to produce:
`193655306161707661080005073394486091998480950338405932486880600467114423441282418165863`

## Julia

### Recursive

`using Memoize function partDiffDiff(n::Int)::Int    isodd(n) ? (n+1)÷2 : n+1end @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    endend 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

`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;  send: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`

## Mathematica / Wolfram Language

`PartitionsP /@ Range[15]PartitionsP[666]PartitionsP[6666]`
Output:
```{1,2,3,5,7,11,15,22,30,42,56,77,101,135,176}
11956824258286445517629485
193655306161707661080005073394486091998480950338405932486880600467114423441282418165863```

## Nim

Translation of: C++
Library: bignum
`import sequtils, strformat, timesimport 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

`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

Library: Phix/mpfr

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

```with javascript_semantics
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
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

` /* Picat 3.0#5 *//* Author: Hakan Kjellerstrand */tablepartition1(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

` /* 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 = 193655306161707661080005073394486091998480950338405932486880600467114423441282418165863. `

## Python

### Python: Mathloger

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

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

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

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

These three REXX versions are recursive.

### version 1

`/*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. */!.= @.;   !.1= 1; !.3= 1; !.5= 1; !.7= 1; !.9= 1 /*assign default value and even digits.*/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
```

### version 2

This REXX version is about   25%   faster than REXX version 1.

The biggest part of the improvement was using the expression     k+k+k     instead of     k*3.

`/*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  /*default values for some low numbers. */!.= @.;  !.1= 1; !.3= 1; !.5= 1; !.7= 1; !.9= 1  /*   "       "    "  all the 1─digit #s*/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+k+k - 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*/               parse var   k   ''  -1  _                   /*obtain K's last decimal dig*/               if !._     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   is identical to the 1st REXX version.

### version 3

This REXX version is about   90%   faster than REXX version 1.

The biggest part of the improvement was using memoization of the expressions     (k+k+k - 1) * k % 2     for all values of (positive)   k   up to   hi.

`/*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  /*default values for some low numbers. */!.= @.;  !.1= 1; !.3= 1; !.5= 1; !.7= 1; !.9= 1  /*   "       "    "  all the 1─digit #s*/w= length( commas(hi) )                          /*W:  is used for aligning the index.  */       do i=1  for hi;  a.i= (i+i+i - 1) * i % 2 /*calculate HI expressions (for partP).*/       end   /*i*/        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 @. !. a.;  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 - a.k                  /*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*/               parse var   k   ''  -1  _                   /*obtain K's last decimal dig*/               if !._     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   is identical to the 1st REXX version.

## 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;        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

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

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.clockvar N = 6666p = List.filled(N+1, null)pd = List.filled(N+1, 0)p[0] = BigInt.onep[1] = BigInt.onepd[0] = 1pd[1] = 1for (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
```