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

(Redirected from Partition Function P)
Partition function P is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.
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.

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

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

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

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

## Phix

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=6666printf(1,"p(%d) = %s (%s)\n",{n,mpz_get_str(partitionsP(n)),elapsed(time()-t0)})`
Output:
```p(6666) = 193655306161707661080005073394486091998480950338405932486880600467114423441282418165863 (0.8s)
```

## Python

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 [2]: %time part(6666)
Wall time: 243 ms
Out[2]: 193655306161707661080005073394486091998480950338405932486880600467114423441282418165863
```

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

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

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(hi)                                    /*W:  is used for aligning the index.  */             do j=lo  to hi                      /*compute a range of  partitionsP.     */             say right(j, w)  partitionsP(j)     /*display index and it's  partitionsP. */             end     /*j*/exit 0                                           /*stick a fork in it,  we're all done. *//*──────────────────────────────────────────────────────────────────────────────────────*/partitionsP: procedure expose @.; parse arg n    /*obtain number (index) for computation*/             if @.n\==0  then return @.n         /*Is it already computed?   Return it. */             #= 0                    do k=1  for n;  _= n - (k*3 - 1) * k % 2;    if _<0  then leave                    if @._==0  then x= partitionsP(_)                    /*use recursion*/                               else x= @._                               /*use found #. */                    _= _ - k                    if _<0     then y= 0                                 /*Negative?    */                               else if @._==0  then y= partitionsP(_)    /*use recursion*/                                               else y= @._               /*use found #. */                    if k//2    then #= # + x + y          /*Odd?  Then   sum    X and Y */                               else #= # - x - y          /*Even?   "  subtract "  "  " */                    end   /*k*/             @.n= #                              /*define the computed partitionsP of N.*/             return #                            /*return  "      "          "      " " */`
output   when using the input of:     6666
```6666 193655306161707661080005073394486091998480950338405932486880600467114423441282418165863
```
output   when using the input of:     66666
```66666 931283431869095717238416063534148471363928685651267074563360050977549251436058076515262033
789845457356081278451246751375974038318319810923102769109446980055567090089060954748065131666952
830623286286824837240058805177319865230913417587625830803675380262765598632742903192085254824621
```

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