Partition function P: Difference between revisions
m (→{{header|REXX}}: added another example calculation.) |
m (→{{header|REXX}}: added HTML tags.) |
||
Line 196:
</pre>
{{out|output|text= when using the input of: <tt> 66666 </tt>}}
<pre>
66666 931283431869095717238416063534148471363928685651267074563360050977549251436058076515262033
789845457356081278451246751375974038318319810923102769109446980055567090089060954748065131666952
830623286286824837240058805177319865230913417587625830803675380262765598632742903192085254824621
</pre>
|
Revision as of 17:16, 28 October 2020
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.
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 sign = ((i-1)%4)<2 ? 1 : -1 psum += sign*partitionsP(n-pd) end psum end
end
n=6666 @time println("p($n) = ", partitionsP(n))</lang>
- Output:
p(6666) = 193655306161707661080005073394486091998480950338405932486880600467114423441282418165863 0.405178 seconds (5.34 M allocations: 114.048 MiB, 13.34% gc time)
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(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 " " " " " */</lang>
- output when using the input of: 6666
6666 193655306161707661080005073394486091998480950338405932486880600467114423441282418165863
- output when using the input of: 66666
66666 931283431869095717238416063534148471363928685651267074563360050977549251436058076515262033 789845457356081278451246751375974038318319810923102769109446980055567090089060954748065131666952 830623286286824837240058805177319865230913417587625830803675380262765598632742903192085254824621