Partition function P: Difference between revisions
m (C++ - removed unnecessary variable) |
m (→{{header|Wren}}: Minor tidy) |
||
(62 intermediate revisions by 30 users not shown) | |||
Line 1: | Line 1: | ||
<includeonly>{{#set:is task=true}}{{#set:is draft=true}}</includeonly>{{infobox_begin}}'' '''{{PAGENAME}}''' '' 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:{{PAGENAME}}|talk page]].<includeonly>[[Category:Draft Programming Tasks]] [[Category:{{{1|Draft Programming Tasks}}}]]</includeonly>{{infobox_end}} |
|||
{{task|Recursion}} |
{{task|Recursion}} |
||
[[Category:Memoization]] |
[[Category:Memoization]] |
||
The [https://mathworld.wolfram.com/PartitionFunctionP.html Partition Function P] |
The [https://mathworld.wolfram.com/PartitionFunctionP.html Partition Function P] is the function P(n), where n∈ℤ, defined as the number of distinct ways in which n can be expressed as the sum of non-increasing positive integers. |
||
Line 16: | Line 15: | ||
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 ... |
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 [https://www.youtube.com/channel/UC1_uAIS3r8Vu6JjXWvastJg Mathologer] made the video, [https://www.youtube.com/watch?v=iJ8pnCO0nTY The hardest "What comes next?" (Euler's pentagonal formula)], where he asks the programmers among his viewers to calculate P(666). The video |
This task may be of popular interest because [https://www.youtube.com/channel/UC1_uAIS3r8Vu6JjXWvastJg Mathologer] made the video, [https://www.youtube.com/watch?v=iJ8pnCO0nTY The hardest "What comes next?" (Euler's pentagonal formula)], where he asks the programmers among his viewers to calculate P(666). The video was viewed more than 100,000 times in the first couple of weeks after its release. |
||
In Wolfram Language, this function has been implemented as PartitionsP. |
In Wolfram Language, this function has been implemented as PartitionsP. |
||
Line 38: | Line 37: | ||
<br><br> |
<br><br> |
||
=={{header|11l}}== |
|||
{{trans|Python: Alternative}} |
|||
<syntaxhighlight lang="11l">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)</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
Partitions: [1, 1, 2, 3, 5, 7, 11, 15, 22, 30, 42, 56, 77, 101, 135] |
|||
193655306161707661080005073394486091998480950338405932486880600467114423441282418165863 |
|||
0.598528 |
|||
</pre> |
|||
=={{header|C}}== |
=={{header|C}}== |
||
{{libheader|GMP}} |
{{libheader|GMP}} |
||
< |
<syntaxhighlight lang="c">#include <stdint.h> |
||
#include <stdlib.h> |
#include <stdlib.h> |
||
#include <stdio.h> |
#include <stdio.h> |
||
Line 77: | Line 114: | ||
(double)(clock() - start) / (double)CLOCKS_PER_SEC); |
(double)(clock() - start) / (double)CLOCKS_PER_SEC); |
||
return 0; |
return 0; |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 85: | Line 122: | ||
</pre> |
</pre> |
||
=={{header|C#|CSharp}}== |
|||
This can also be done using BigIntegers, but will take around five times longer. Since only adding and subtracting are required, some simple routines can be created to handle the tabulations. Also note that detecting odd and even numbers on each loop iteration is avoided by coding four increments per loop. |
|||
<syntaxhighlight lang="csharp">using System; |
|||
class Program { |
|||
const long Lm = (long)1e18; |
|||
const string Fm = "D18"; |
|||
// provides 5 x 18 = 90 decimal digits |
|||
struct LI { public long lo, ml, mh, hi, tp; } |
|||
static void inc(ref LI d, LI s) { // d += s |
|||
if ((d.lo += s.lo) >= Lm) { d.ml++; d.lo -= Lm; } |
|||
if ((d.ml += s.ml) >= Lm) { d.mh++; d.ml -= Lm; } |
|||
if ((d.mh += s.mh) >= Lm) { d.hi++; d.mh -= Lm; } |
|||
if ((d.hi += s.hi) >= Lm) { d.tp++; d.hi -= Lm; } |
|||
d.tp += s.tp; |
|||
} |
|||
static void dec(ref LI d, LI s) { // d -= s |
|||
if ((d.lo -= s.lo) < 0) { d.ml--; d.lo += Lm; } |
|||
if ((d.ml -= s.ml) < 0) { d.mh--; d.ml += Lm; } |
|||
if ((d.mh -= s.mh) < 0) { d.hi--; d.mh += Lm; } |
|||
if ((d.hi -= s.hi) < 0) { d.tp--; d.hi += Lm; } |
|||
d.tp -= s.tp; |
|||
} |
|||
static LI set(long s) { LI d; |
|||
d.lo = s; d.ml = d.mh = d.hi = d.tp = 0; return d; } |
|||
static string fmt(LI x) { // returns formatted string value of x |
|||
if (x.tp > 0) return x.tp.ToString() + x.hi.ToString(Fm) + x.mh.ToString(Fm) + x.ml.ToString(Fm) + x.lo.ToString(Fm); |
|||
if (x.hi > 0) return x.hi.ToString() + x.mh.ToString(Fm) + x.ml.ToString(Fm) + x.lo.ToString(Fm); |
|||
if (x.mh > 0) return x.mh.ToString() + x.ml.ToString(Fm) + x.lo.ToString(Fm); |
|||
if (x.ml > 0) return x.ml.ToString() + x.lo.ToString(Fm); |
|||
return x.lo.ToString(); |
|||
} |
|||
static LI partcount(int n) { |
|||
var P = new LI[n + 1]; P[0] = set(1); |
|||
for (int i = 1; i <= n; i++) { |
|||
int k = 0, d = -2, j = i; |
|||
LI x = set(0); |
|||
while (true) { |
|||
if ((j -= (d += 3) -k) >= 0) inc(ref x, P[j]); else break; |
|||
if ((j -= ++k) >= 0) inc(ref x, P[j]); else break; |
|||
if ((j -= (d += 3) -k) >= 0) dec(ref x, P[j]); else break; |
|||
if ((j -= ++k) >= 0) dec(ref x, P[j]); else break; |
|||
} |
|||
P[i] = x; |
|||
} |
|||
return P[n]; |
|||
} |
|||
static void Main(string[] args) { |
|||
var sw = System.Diagnostics.Stopwatch.StartNew (); |
|||
var res = partcount(6666); sw.Stop(); |
|||
Console.Write("{0} {1} ms", fmt(res), sw.Elapsed.TotalMilliseconds); |
|||
} |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre>193655306161707661080005073394486091998480950338405932486880600467114423441282418165863 12.9365 ms</pre> |
|||
Timing from Tio.run. |
|||
=={{header|C++}}== |
=={{header|C++}}== |
||
===GMP version=== |
|||
{{libheader|GMP}} |
{{libheader|GMP}} |
||
< |
<syntaxhighlight lang="cpp">#include <chrono> |
||
#include <iostream> |
#include <iostream> |
||
#include <vector> |
#include <vector> |
||
Line 100: | Line 202: | ||
for (int i = 1; i <= n; ++i) { |
for (int i = 1; i <= n; ++i) { |
||
for (int k = 1;; ++k) { |
for (int k = 1;; ++k) { |
||
int |
int j = (k * (3*k - 1))/2; |
||
if ( |
if (j > i) |
||
break; |
break; |
||
if (k & 1) |
if (k & 1) |
||
p[i] += p[i - |
p[i] += p[i - j]; |
||
else |
else |
||
p[i] -= p[i - |
p[i] -= p[i - j]; |
||
j = (k * (3*k + 1))/2; |
|||
if ( |
if (j > i) |
||
break; |
break; |
||
if (k & 1) |
if (k & 1) |
||
p[i] += p[i - |
p[i] += p[i - j]; |
||
else |
else |
||
p[i] -= p[i - |
p[i] -= p[i - j]; |
||
} |
} |
||
} |
} |
||
Line 126: | Line 228: | ||
std::cout << result << '\n'; |
std::cout << result << '\n'; |
||
std::cout << "elapsed time: " << ms.count() << " milliseconds\n"; |
std::cout << "elapsed time: " << ms.count() << " milliseconds\n"; |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 132: | Line 234: | ||
193655306161707661080005073394486091998480950338405932486880600467114423441282418165863 |
193655306161707661080005073394486091998480950338405932486880600467114423441282418165863 |
||
elapsed time: 8.99497 milliseconds |
elapsed time: 8.99497 milliseconds |
||
</pre> |
|||
===Non GMP version=== |
|||
{{trans|C#}} |
|||
<syntaxhighlight lang="cpp">#include <chrono> |
|||
#include <iostream> |
|||
using namespace std; |
|||
using namespace chrono; |
|||
const long long Lm = (long)1e18; |
|||
const int Fm = 18; |
|||
struct LI { long long lo, ml, mh, hi, tp; }; |
|||
LI set(long long s) { LI d; |
|||
d.lo = s; d.ml = d.mh = d.hi = d.tp = 0; return d; } |
|||
void inc(LI& d, LI s) { // d += s |
|||
if ((d.lo += s.lo) >= Lm) { d.ml++; d.lo -= Lm; } |
|||
if ((d.ml += s.ml) >= Lm) { d.mh++; d.ml -= Lm; } |
|||
if ((d.mh += s.mh) >= Lm) { d.hi++; d.mh -= Lm; } |
|||
if ((d.hi += s.hi) >= Lm) { d.tp++; d.hi -= Lm; } |
|||
d.tp += s.tp; |
|||
} |
|||
void dec(LI& d, LI s) { // d -= s |
|||
if ((d.lo -= s.lo) < 0) { d.ml--; d.lo += Lm; } |
|||
if ((d.ml -= s.ml) < 0) { d.mh--; d.ml += Lm; } |
|||
if ((d.mh -= s.mh) < 0) { d.hi--; d.mh += Lm; } |
|||
if ((d.hi -= s.hi) < 0) { d.tp--; d.hi += Lm; } |
|||
d.tp -= s.tp; |
|||
} |
|||
inline string sf(long long n) { |
|||
int len = Fm; |
|||
string result(len--, '0'); |
|||
for (long long i = n; len >= 0 && i > 0; --len, i /= 10) |
|||
result[len] = '0' + i % 10; |
|||
return result; |
|||
} |
|||
string fmt(LI x) { // returns formatted string value of x |
|||
if (x.tp > 0) return to_string(x.tp) + sf(x.hi) + sf(x.mh) + sf(x.ml) + sf(x.lo); |
|||
if (x.hi > 0) return to_string(x.hi) + sf(x.mh) + sf(x.ml) + sf(x.lo); |
|||
if (x.mh > 0) return to_string(x.mh) + sf(x.ml) + sf(x.lo); |
|||
if (x.ml > 0) return to_string(x.ml) + sf(x.lo); |
|||
return to_string(x.lo); |
|||
} |
|||
LI partcount(int n) { LI p[n + 1]; p[0] = set(1); |
|||
for (int i = 1; i <= n; i++) { |
|||
int k = 0, d = -2, j = i; LI x = set(0); while (true) { |
|||
if ((j -= (d += 3) - k) >= 0) inc(x, p[j]); else break; |
|||
if ((j -= ++k) >= 0) inc(x, p[j]); else break; |
|||
if ((j -= (d += 3) - k) >= 0) dec(x, p[j]); else break; |
|||
if ((j -= ++k) >= 0) dec(x, p[j]); else break; |
|||
} p[i] = x; } return p[n]; } |
|||
int main() { |
|||
auto start = steady_clock::now(); |
|||
auto result = partcount(6666); |
|||
auto end = steady_clock::now(); |
|||
duration<double, milli> ms(end - start); |
|||
cout << fmt(result) << " " << ms.count() << " ms"; |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
Timing from Tio.run, but execution time can't be directly compared to the GMP version, as GMP isn't accessible at Tio.run. |
|||
<pre>193655306161707661080005073394486091998480950338405932486880600467114423441282418165863 7.32706 ms</pre> |
|||
=={{header|Delphi}}== |
|||
{{libheader| System.SysUtils}} |
|||
{{libheader| Velthuis.BigIntegers}} |
|||
{{libheader| System.Diagnostics}} |
|||
{{Trans|Go}} |
|||
<syntaxhighlight lang="delphi"> |
|||
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.</syntaxhighlight> |
|||
{{out}} |
|||
<pre>p[6666] = 193655306161707661080005073394486091998480950338405932486880600467114423441282418165863 |
|||
Took 131ms</pre> |
|||
=={{header|Elixir}}== |
|||
Loosely based on the Erlang version. |
|||
<syntaxhighlight lang="elixir"> |
|||
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 |
|||
end |
|||
end |
|||
Partition.init |
|||
IO.puts Partition.p 6666 |
|||
</syntaxhighlight> |
|||
{{Out}} |
|||
<pre> |
|||
$ time ./partition6666.ex |
|||
193655306161707661080005073394486091998480950338405932486880600467114423441282418165863 |
|||
real 0m1.106s |
|||
user 0m1.191s |
|||
sys 0m0.116s |
|||
</pre> |
|||
=={{header|Erlang}}== |
|||
<syntaxhighlight lang="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. |
|||
</syntaxhighlight> |
|||
{{Out}} |
|||
<pre> |
|||
$ time ./partition6666.erl |
|||
193655306161707661080005073394486091998480950338405932486880600467114423441282418165863 |
|||
real 0m0.480s |
|||
user 0m0.490s |
|||
sys 0m0.080s |
|||
</pre> |
|||
=={{header|F_Sharp|F#}}== |
|||
An implementation of the formula in the task description. P(123456) is included for comparison with the largest value in the related task. |
|||
<syntaxhighlight lang="fsharp"> |
|||
// 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) |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
666->11956824258286445517629485 |
|||
6666->193655306161707661080005073394486091998480950338405932486880600467114423441282418165863 |
|||
Real: 00:00:00.096 |
|||
123456->30817659578536496678545317146533980855296613274507139217608776782063054452191537379312358383342446230621170608408020911309259407611257151683372221925128388387168451943800027128045369650890220060901494540459081545445020808726917371699102825508039173543836338081612528477859613355349851184591540231790254269948278726548570660145691076819912972162262902150886818986555127204165221706149989 |
|||
</pre> |
</pre> |
||
=={{header|Factor}}== |
=={{header|Factor}}== |
||
{{works with|Factor|0.99 2020-08-14}} |
{{works with|Factor|0.99 2020-08-14}} |
||
< |
<syntaxhighlight lang="factor">USING: kernel lists lists.lazy math sequences sequences.extras ; |
||
! Compute the nth pentagonal number. |
! Compute the nth pentagonal number. |
||
Line 158: | Line 527: | ||
: partitions ( m -- n ) |
: partitions ( m -- n ) |
||
[ seq swap [ <= ] curry lwhile list>array ] |
[ seq swap [ <= ] curry lwhile list>array ] |
||
[ V{ 1 } clone swap [ step ] times last nip ] bi ;</ |
[ V{ 1 } clone swap [ step ] times last nip ] bi ;</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 165: | Line 534: | ||
Running time: 0.084955341 seconds |
Running time: 0.084955341 seconds |
||
193655306161707661080005073394486091998480950338405932486880600467114423441282418165863 |
|||
</pre> |
|||
=={{header|FreeBASIC}}== |
|||
===Unsiged 64bit version=== |
|||
{{trans|Python}} |
|||
<syntaxhighlight lang="freebasic">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</syntaxhighlight> |
|||
{{out}} |
|||
<pre>PartitionsP: 1 1 2 3 5 7 11 15 22 30 42 56 77</pre> |
|||
===Big numbers version=== |
|||
{{libheader|GMP}} |
|||
[https://rosettacode.org/wiki/9_billion_names_of_God_the_integer#FreeBASIC From the 9_billion_names_of_God_the_integer entry] |
|||
<syntaxhighlight lang="freebasic">' 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 n |
|||
Dim As ZString Ptr ans |
|||
Dim As Double t = Timer |
|||
ReDim big_p(max) As Mpz_ptr |
|||
For 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 buffer |
|||
While InKey <> "" : Wend |
|||
Print : Print "hit any key to end program" |
|||
Sleep |
|||
End</syntaxhighlight> |
|||
{{out}} |
|||
<pre>PartitionsP(6666) = 193655306161707661080005073394486091998480950338405932486880600467114423441282418165863 |
|||
time = 32.97 ms</pre> |
|||
=={{header|Frink}}== |
|||
Frink has a built-in function for counting partitions that uses Euler's pentagonal method. It works for arbitrarily-large integers and caches results. |
|||
<syntaxhighlight lang="frink">println[partitionCount[6666]]</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
193655306161707661080005073394486091998480950338405932486880600467114423441282418165863 |
193655306161707661080005073394486091998480950338405932486880600467114423441282418165863 |
||
</pre> |
</pre> |
||
Line 171: | Line 662: | ||
{{trans|Julia}} |
{{trans|Julia}} |
||
I also tried using Euler's generating function but it was about 20 times slower than this version. |
I also tried using Euler's generating function but it was about 20 times slower than this version. |
||
< |
<syntaxhighlight lang="go">package main |
||
import ( |
import ( |
||
Line 229: | Line 720: | ||
fmt.Printf("p[%d)] = %d\n", N, p[N]) |
fmt.Printf("p[%d)] = %d\n", N, p[N]) |
||
fmt.Printf("Took %s\n", time.Since(start)) |
fmt.Printf("Took %s\n", time.Since(start)) |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 236: | Line 727: | ||
Took 54.82683ms |
Took 54.82683ms |
||
</pre> |
</pre> |
||
=={{header|Java}}== |
|||
<syntaxhighlight lang="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]; |
|||
} |
|||
}</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
P(6666) = 193655306161707661080005073394486091998480950338405932486880600467114423441282418165863 |
|||
elapsed time: 59 milliseconds |
|||
</pre> |
|||
=={{header|JavaScript}}== |
|||
<syntaxhighlight lang="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) |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
p(6666) = 193655306161707661080005073394486091998480950338405932486880600467114423441282418165863 |
|||
Computation time in ms: 26 |
|||
</pre> |
|||
=={{header|Haskell}}== |
|||
<syntaxhighlight lang="haskell">{-# LANGUAGE DeriveFunctor #-} |
|||
------------------------------------------------------------ |
|||
-- memoization utilities |
|||
data Memo a = Node a (Memo a) (Memo a) |
|||
deriving (Functor) |
|||
memo :: Integral a => Memo p -> a -> p |
|||
memo (Node a l r) n |
|||
| n == 0 = a |
|||
| odd n = memo l (n `div` 2) |
|||
| otherwise = memo r (n `div` 2 - 1) |
|||
nats :: Memo Int |
|||
nats = |
|||
Node |
|||
0 |
|||
((+ 1) . (* 2) <$> nats) |
|||
((* 2) . (+ 1) <$> nats) |
|||
------------------------------------------------------------ |
|||
-- calculating partitions |
|||
partitions :: Memo Integer |
|||
partitions = partitionP <$> nats |
|||
partitionP :: Int -> Integer |
|||
partitionP 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 :: [Int] |
|||
ofsets = scanl1 (+) $ mix [1, 3 ..] [1, 2 ..] |
|||
where |
|||
mix a b = concat $ zipWith (\x y -> [x, y]) a b |
|||
main :: IO () |
|||
main = print $ partitionP 6666</syntaxhighlight> |
|||
<pre>*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</pre> |
|||
=={{header|J}}== |
|||
Solution stolen [https://code.jsoftware.com/wiki/Essays/Partitions#The_Number_of_Partitions verbatim from the J Wiki]. Note the use of memoization (M.) for efficiency: |
|||
<syntaxhighlight lang="j"> pn =: -/@(+/)@:($:"0)@rec ` (x:@(0&=)) @. (0>:]) M. |
|||
rec=: - (-: (*"1) _1 1 +/ 3 * ]) @ (>:@i.@>.@%:@((2%3)&*))</syntaxhighlight> |
|||
{{out}} |
|||
<pre> pn 6 |
|||
11 |
|||
pn 66 |
|||
2323520 |
|||
pn 666 |
|||
11956824258286445517629485 |
|||
pn 6666 |
|||
193655306161707661080005073394486091998480950338405932486880600467114423441282418165863</pre> |
|||
=={{header|jq}}== |
|||
Translation of: Python:Alternative |
|||
<syntaxhighlight lang="jq">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))]</syntaxhighlight> |
|||
{{out}} |
|||
<pre>[1,2,3,5,7,11,15,22,30,42,56,77,101,135]</pre> |
|||
Using gojq 0.12.11, `partitions(6666)` yields (in about 12 minutes (u+s) on a 3GHz machine): |
|||
193655306161707661080005073394486091998480950338405932486880600467114423441282418165863 |
|||
The integer precision of the C implementation of jq is insufficient for computing ``partitions(6666)``, but as a test |
|||
of the BigInt.jq library for jq, 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" |
|||
Curiously, the u+s time is 7m3s, which is significantly less than the above-mentioned gojq time, even though the BigInt library is written in jq. |
|||
=== Recursive === |
|||
{{trans|Julia}} with memoization |
|||
<syntaxhighlight lang="jq">def partDiffDiff($n): |
|||
if ($n % 2) == 1 then ($n+1) / 2 else $n+1 end; |
|||
# in: {n, partDiffMemo} |
|||
# out: object with possibly updated memoization |
|||
def 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 memoization |
|||
def 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 |
|||
</syntaxhighlight> |
|||
Using gojq, the above program takes 41.35 seconds (u+s) on a 3GHz Mac Mini to produce:<pre> |
|||
193655306161707661080005073394486091998480950338405932486880600467114423441282418165863</pre> |
|||
=={{header|Julia}}== |
=={{header|Julia}}== |
||
===Recursive=== |
===Recursive=== |
||
< |
<syntaxhighlight lang="julia">using Memoize |
||
function partDiffDiff(n::Int)::Int |
function partDiffDiff(n::Int)::Int |
||
Line 271: | Line 1,008: | ||
n=6666 |
n=6666 |
||
@time println("p($n) = ", partitionsP(n))</ |
@time println("p($n) = ", partitionsP(n))</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>p(6666) = 193655306161707661080005073394486091998480950338405932486880600467114423441282418165863 |
<pre>p(6666) = 193655306161707661080005073394486091998480950338405932486880600467114423441282418165863 |
||
0.260310 seconds (3.58 M allocations: 77.974 MiB, 8.54% gc time)</pre> |
0.260310 seconds (3.58 M allocations: 77.974 MiB, 8.54% gc time)</pre> |
||
=={{header|Lingo}}== |
|||
Lingo natively only supports 32 bit integers, so P(6666) would be way too big. |
|||
<syntaxhighlight lang="lingo">-- returns number of partitions of n |
|||
on partitions(n, res_table) |
|||
if n < 2 then return 1 |
|||
if voidP(res_table) then |
|||
res_table = [] |
|||
res_table[n] = 0 |
|||
else if res_table[n] then |
|||
return res_table[n] |
|||
end if |
|||
res = 0 |
|||
i = 0 |
|||
param = 1 |
|||
repeat while param <= n |
|||
if i mod 4 < 2 then |
|||
res = res + partitions(n - param, res_table) |
|||
else |
|||
res = res - partitions(n - param, res_table) |
|||
end if |
|||
if i mod 2 then |
|||
param = param + i + 2 |
|||
else |
|||
param = param + i / 2 + 1 |
|||
end if |
|||
i = i + 1 |
|||
end repeat |
|||
res_table[n] = res |
|||
return res |
|||
end</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
ms = _system.milliseconds |
|||
put "P(121):", partitions(121) |
|||
put "Time (ms):", _system.milliseconds - ms |
|||
-- P(121): 2056148051 |
|||
-- Time (ms): 3 |
|||
</pre> |
|||
=={{header|Maple}}== |
|||
<syntaxhighlight lang="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; |
|||
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</syntaxhighlight> |
|||
=={{header|Mathematica}} / {{header|Wolfram Language}}== |
|||
<syntaxhighlight lang="mathematica">PartitionsP /@ Range[15] |
|||
PartitionsP[666] |
|||
PartitionsP[6666]</syntaxhighlight> |
|||
{{out}} |
|||
<pre>{1,2,3,5,7,11,15,22,30,42,56,77,101,135,176} |
|||
11956824258286445517629485 |
|||
193655306161707661080005073394486091998480950338405932486880600467114423441282418165863</pre> |
|||
=={{header|Nim}}== |
|||
{{trans|C++}} |
|||
{{libheader|bignum}} |
|||
<syntaxhighlight lang="nim">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"</syntaxhighlight> |
|||
{{out}} |
|||
<pre>193655306161707661080005073394486091998480950338405932486880600467114423441282418165863 |
|||
Elapsed time: 9.73 ms</pre> |
|||
=={{header|Perl}}== |
|||
<syntaxhighlight lang="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";</syntaxhighlight> |
|||
{{out}} |
|||
<pre>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 |
|||
</pre> |
|||
=={{header|Phix}}== |
=={{header|Phix}}== |
||
{{libheader|Phix/mpfr}} |
|||
Not exactly super-fast, but this sort of stuff is not really what Phix does best. |
Not exactly super-fast, but this sort of stuff is not really what Phix does best. |
||
<!--<syntaxhighlight lang="phix">(phixonline)--> |
|||
<lang Phix>function partDiffDiff(integer n) |
|||
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span> |
|||
return (n+1)/(1+and_bits(n,1)) |
|||
<span style="color: #008080;">function</span> <span style="color: #000000;">partDiffDiff</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">)</span> |
|||
end function |
|||
<span style="color: #008080;">return</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">n</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)/(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">+</span><span style="color: #7060A8;">and_bits</span><span style="color: #0000FF;">(</span><span style="color: #000000;">n</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">))</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
|||
sequence pd = {1} |
|||
function partDiff(integer n) |
|||
<span style="color: #004080;">sequence</span> <span style="color: #000000;">pd</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">1</span><span style="color: #0000FF;">}</span> |
|||
while n>length(pd) do |
|||
<span style="color: #008080;">function</span> <span style="color: #000000;">partDiff</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">)</span> |
|||
pd &= pd[$] + partDiffDiff(length(pd)) |
|||
<span style="color: #008080;">while</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">></span><span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pd</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span> |
|||
end while |
|||
<span style="color: #000000;">pd</span> <span style="color: #0000FF;">&=</span> <span style="color: #000000;">pd</span><span style="color: #0000FF;">[$]</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">partDiffDiff</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pd</span><span style="color: #0000FF;">))</span> |
|||
return pd[max(1,n)] |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">while</span> |
|||
end function |
|||
<span style="color: #008080;">return</span> <span style="color: #000000;">pd</span><span style="color: #0000FF;">[</span><span style="color: #7060A8;">max</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">n</span><span style="color: #0000FF;">)]</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
|||
include mpfr.e |
|||
<span style="color: #008080;">include</span> <span style="color: #004080;">mpfr</span><span style="color: #0000FF;">.</span><span style="color: #000000;">e</span> |
|||
sequence pn = {mpz_init(1)} |
|||
function partitionsP(integer n) |
|||
<span style="color: #004080;">sequence</span> <span style="color: #000000;">pn</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)}</span> |
|||
mpz res = mpz_init(1) |
|||
<span style="color: #008080;">function</span> <span style="color: #000000;">partitionsP</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">)</span> |
|||
while n>length(pn) do |
|||
<span style="color: #004080;">mpz</span> <span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span> |
|||
integer nn = length(pn)+1 |
|||
<span style="color: #008080;">while</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">></span><span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pn</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span> |
|||
mpz psum = mpz_init(0) |
|||
<span style="color: #004080;">integer</span> <span style="color: #000000;">nn</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pn</span><span style="color: #0000FF;">)+</span><span style="color: #000000;">1</span> |
|||
for i=1 to nn do |
|||
<span style="color: #004080;">mpz</span> <span style="color: #000000;">psum</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">mpz_init</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">)</span> |
|||
integer pd = partDiff(i) |
|||
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">nn</span> <span style="color: #008080;">do</span> |
|||
if pd>nn then exit end if |
|||
<span style="color: #004080;">integer</span> <span style="color: #000000;">pd</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">partDiff</span><span style="color: #0000FF;">(</span><span style="color: #000000;">i</span><span style="color: #0000FF;">)</span> |
|||
integer sgn = iff(remainder(i-1,4)<2 ? 1 : -1) |
|||
<span style="color: #008080;">if</span> <span style="color: #000000;">pd</span><span style="color: #0000FF;">></span><span style="color: #000000;">nn</span> <span style="color: #008080;">then</span> <span style="color: #008080;">exit</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
mpz pnmpd = pn[max(1,nn-pd)] |
|||
<span style="color: #004080;">integer</span> <span style="color: #000000;">sgn</span> <span style="color: #0000FF;">=</span> <span style="color: #008080;">iff</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">remainder</span><span style="color: #0000FF;">(</span><span style="color: #000000;">i</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">4</span><span style="color: #0000FF;">)<</span><span style="color: #000000;">2</span> <span style="color: #0000FF;">?</span> <span style="color: #000000;">1</span> <span style="color: #0000FF;">:</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span> |
|||
if sgn=-1 then |
|||
<span style="color: #004080;">mpz</span> <span style="color: #000000;">pnmpd</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">pn</span><span style="color: #0000FF;">[</span><span style="color: #7060A8;">max</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">nn</span><span style="color: #0000FF;">-</span><span style="color: #000000;">pd</span><span style="color: #0000FF;">)]</span> |
|||
mpz_sub(psum,psum,pnmpd) |
|||
<span style="color: #008080;">if</span> <span style="color: #000000;">sgn</span><span style="color: #0000FF;">=-</span><span style="color: #000000;">1</span> <span style="color: #008080;">then</span> |
|||
else |
|||
<span style="color: #7060A8;">mpz_sub</span><span style="color: #0000FF;">(</span><span style="color: #000000;">psum</span><span style="color: #0000FF;">,</span><span style="color: #000000;">psum</span><span style="color: #0000FF;">,</span><span style="color: #000000;">pnmpd</span><span style="color: #0000FF;">)</span> |
|||
mpz_add(psum,psum,pnmpd) |
|||
<span style="color: #008080;">else</span> |
|||
<span style="color: #7060A8;">mpz_add</span><span style="color: #0000FF;">(</span><span style="color: #000000;">psum</span><span style="color: #0000FF;">,</span><span style="color: #000000;">psum</span><span style="color: #0000FF;">,</span><span style="color: #000000;">pnmpd</span><span style="color: #0000FF;">)</span> |
|||
end for |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span> |
|||
pn = append(pn,psum) |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span> |
|||
end while |
|||
<span style="color: #000000;">pn</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">append</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pn</span><span style="color: #0000FF;">,</span><span style="color: #000000;">psum</span><span style="color: #0000FF;">)</span> |
|||
return pn[max(1,n)] |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">while</span> |
|||
end function |
|||
<span style="color: #008080;">return</span> <span style="color: #000000;">pn</span><span style="color: #0000FF;">[</span><span style="color: #7060A8;">max</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">n</span><span style="color: #0000FF;">)]</span> |
|||
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span> |
|||
atom t0 = time() |
|||
integer n=6666 |
|||
<span style="color: #004080;">atom</span> <span style="color: #000000;">t0</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">time</span><span style="color: #0000FF;">()</span> |
|||
printf(1,"p(%d) = %s (%s)\n",{n,mpz_get_str(partitionsP(n)),elapsed(time()-t0)})</lang> |
|||
<span style="color: #004080;">integer</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">=</span><span style="color: #000000;">6666</span> |
|||
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"p(%d) = %s (%s)\n"</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">n</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">mpz_get_str</span><span style="color: #0000FF;">(</span><span style="color: #000000;">partitionsP</span><span style="color: #0000FF;">(</span><span style="color: #000000;">n</span><span style="color: #0000FF;">)),</span><span style="color: #7060A8;">elapsed</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">time</span><span style="color: #0000FF;">()-</span><span style="color: #000000;">t0</span><span style="color: #0000FF;">)})</span> |
|||
<!--</syntaxhighlight>--> |
|||
{{out}} |
{{out}} |
||
<pre> |
<pre> |
||
Line 322: | Line 1,213: | ||
</pre> |
</pre> |
||
=={{header|Picat}}== |
|||
<syntaxhighlight lang="picat"> |
|||
/* 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. |
|||
</syntaxhighlight> |
|||
=={{header|PicoLisp}}== |
|||
Based on the Erlang implementation. |
|||
<syntaxhighlight lang="picolisp"> |
|||
(de gpentagonals (Max) |
|||
(make |
|||
(let (N 0 M 1) |
|||
(loop |
|||
(inc 'N (if (=0 (& M 1)) (>> 1 M) M)) |
|||
(T (> N Max)) |
|||
(link N) |
|||
(inc 'M))))) |
|||
(de p (N) |
|||
(cache '(NIL) N |
|||
(if (=0 N) |
|||
1 |
|||
(let (Sum 0 Sgn 0) |
|||
(for G (gpentagonals N) |
|||
((if (< Sgn 2) 'inc 'dec) 'Sum (p (- N G))) |
|||
(setq Sgn (& 3 (inc Sgn)))) |
|||
Sum)))) |
|||
</syntaxhighlight> |
|||
{{Out}} |
|||
<pre> |
|||
: (bench (p 6666)) |
|||
0.959 sec |
|||
-> 193655306161707661080005073394486091998480950338405932486880600467114423441282418165863 |
|||
</pre> |
|||
=={{header|Prolog}}== |
|||
<syntaxhighlight lang="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 = 1936553061617076610800050733944860919984809503384 |
|||
05932486880600467114423441282418165863. |
|||
</syntaxhighlight> |
|||
=={{header|Python}}== |
=={{header|Python}}== |
||
===Python: Mathloger=== |
|||
This follows the algorithm from the Mathloger video closely |
This follows the algorithm from the Mathloger video closely |
||
< |
<syntaxhighlight lang="python">from itertools import islice |
||
def posd(): |
def posd(): |
||
Line 376: | Line 1,345: | ||
print(" pos_gen:", list(islice(pos_gen(), 10))) |
print(" pos_gen:", list(islice(pos_gen(), 10))) |
||
print(" plus_minus:", list(islice(plus_minus(), 15))) |
print(" plus_minus:", list(islice(plus_minus(), 15))) |
||
print("\nPartitions:", [part(x) for x in range(15)])</ |
print("\nPartitions:", [part(x) for x in range(15)])</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 388: | Line 1,357: | ||
;Stretch goal: |
;Stretch goal: |
||
From command line after running the above |
From command line after running the above |
||
<pre>In [ |
<pre>In [1]: part(6666) |
||
Out[1]: 193655306161707661080005073394486091998480950338405932486880600467114423441282418165863 |
|||
Wall time: 243 ms |
|||
Out[2]: 193655306161707661080005073394486091998480950338405932486880600467114423441282418165863 |
|||
In [2]: %timeit part(6666) |
|||
103 ms ± 477 µs per loop (mean ± std. dev. of 7 runs, 10 loops each) |
|||
</pre> |
</pre> |
||
===Python: Mathloger video prime generator=== |
====Python: Mathloger video prime generator==== |
||
Add the following code after that above |
Add the following code after that above |
||
< |
<syntaxhighlight lang="python">def par_primes(): |
||
"Prime number generator from the partition machine" |
"Prime number generator from the partition machine" |
||
p = [1] |
p = [1] |
||
Line 412: | Line 1,383: | ||
yield p |
yield p |
||
print("\nPrimes:", list(islice(par_primes(), 15)))</ |
print("\nPrimes:", list(islice(par_primes(), 15)))</syntaxhighlight> |
||
{{Out}} |
{{Out}} |
||
<pre>Primes: [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]</pre> |
<pre>Primes: [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]</pre> |
||
===Python: Alternative=== |
|||
{{trans|C++}} |
|||
<syntaxhighlight lang="python">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)])</syntaxhighlight> |
|||
{{out}} |
|||
<pre>Partitions: [1, 1, 2, 3, 5, 7, 11, 15, 22, 30, 42, 56, 77, 101, 135]</pre> |
|||
;Stretch goal: |
|||
From command line after running the above |
|||
<pre>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) |
|||
</pre> |
|||
=={{header|Quackery}}== |
|||
<code>0 partitions</code> returns <code>1</code> as per [https://oeis.org/A000041 oeis.org/A000041 (Partitions of n)]. |
|||
This is a naive recursive solution, so computing the partitions of 6666 would take a hideously long time. |
|||
<syntaxhighlight lang="Quackery"> [ 1 swap |
|||
dup 0 = iff drop done |
|||
[ 2dup = iff [ 2drop 1 ] done |
|||
2dup > iff [ 2drop 0 ] done |
|||
2dup dip 1+ recurse |
|||
unrot over - recurse + ] ] is partitions ( n --> n ) |
|||
say "Partitions of 0 to 29" cr |
|||
30 times [ i^ partitions echo sp ] |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre>Partitions of 0 to 29 |
|||
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 2436 3010 3718 4565 |
|||
</pre> |
|||
=={{header|Racket}}== |
|||
Backwards range was used to get responsive feedback for progress. |
|||
<syntaxhighlight lang="racket">#lang racket |
|||
(require math/number-theory) |
|||
(define σ |
|||
(let ((memo (make-hash))) |
|||
(λ (z) |
|||
(hash-ref! memo z |
|||
(λ () (apply + (divisors z))))))) |
|||
(define p |
|||
(let ((memo (make-hash '((0 . 1))))) |
|||
(λ (n) |
|||
(hash-ref! |
|||
memo n |
|||
(λ () |
|||
(let ((r (if (zero? n) 1 |
|||
(/ (for/sum ((k (in-range (sub1 n) -1 -1))) |
|||
(* (σ (- n k)) |
|||
(p k))) |
|||
n)))) |
|||
(when (zero? (modulo n 1000)) (displayln (cons n r) (current-error-port))) |
|||
r)))))) |
|||
(map p (range 1 30)) |
|||
(p 666) |
|||
(p 1000) |
|||
(p 10000)</syntaxhighlight> |
|||
{{out}} |
|||
<pre>'(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 2436 3010 3718 4565) |
|||
11956824258286445517629485 |
|||
(1000 . 24061467864032622473692149727991) |
|||
24061467864032622473692149727991 |
|||
(2000 . 4720819175619413888601432406799959512200344166) |
|||
(3000 . 496025142797537184410324879054927095334462742231683423624) |
|||
(4000 . 1024150064776551375119256307915896842122498030313150910234889093895) |
|||
(5000 . 169820168825442121851975101689306431361757683049829233322203824652329144349) |
|||
(6000 . 4671727531970209092971024643973690643364629153270037033856605528925072405349246129) |
|||
(7000 . 32856930803440615786280925635924166861950151574532240659699032157432236394374450791229199) |
|||
(8000 . 78360264351568349490593145013364599719010769352985864331118600209417827764524450990388402844164) |
|||
(9000 . 77133638117808884907320791427403134961639798322072034262647713694605367979684296948790335590435626459) |
|||
(10000 . 36167251325636293988820471890953695495016030339315650422081868605887952568754066420592310556052906916435144) |
|||
36167251325636293988820471890953695495016030339315650422081868605887952568754066420592310556052906916435144</pre> |
|||
=={{header|Raku}}== |
=={{header|Raku}}== |
||
{{works with|Rakudo|2020.09}} |
{{works with|Rakudo|2020.09}} |
||
Not the fastest, but it gets the job done. |
Not the fastest, but it gets the job done. |
||
<lang |
<syntaxhighlight lang="raku" line>my @P = 1, { p(++$) } … *; |
||
my @i = lazy [\+] flat 1, ( |
my @i = lazy [\+] flat 1, ( 1..* Z (1..*).map: * × 2 + 1 ); |
||
sub p ($n) { sum @P[$n X- @i |
sub p ($n) { sum @P[$n X- @i] Z× (flat (1, 1, -1, -1) xx *) } |
||
put @P[^26]; |
put @P[^26]; |
||
put @P[6666];</ |
put @P[6666];</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre>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 |
<pre>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 |
||
Line 431: | Line 1,520: | ||
=={{header|REXX}}== |
=={{header|REXX}}== |
||
These three REXX versions are recursive. |
|||
=== version 1 === |
|||
<lang rexx>/*REXX program calculates and displays a specific value (or a range of) partitionsP(N).*/ |
|||
<syntaxhighlight 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*/ |
numeric digits 1000 /*able to handle some ginormous numbers*/ |
||
parse arg lo hi . /*obtain optional arguments from the CL*/ |
parse arg lo hi . /*obtain optional arguments from the CL*/ |
||
Line 438: | Line 1,528: | ||
if hi=='' | hi=="," then hi= lo /* " " " " " " */ |
if hi=='' | hi=="," then hi= lo /* " " " " " " */ |
||
@.= 0; @.0= 1; @.1= 1; @.2= 2; @.3= 3; @.4= 5 /*assign default value and low values. */ |
@.= 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. */ |
w= length( commas(hi) ) /*W: is used for aligning the index. */ |
||
Line 447: | Line 1,538: | ||
commas: parse arg ?; do jc=length(?)-3 to 1 by -3; ?=insert(',', ?, jc); end; return ? |
commas: parse arg ?; do jc=length(?)-3 to 1 by -3; ?=insert(',', ?, jc); end; return ? |
||
/*──────────────────────────────────────────────────────────────────────────────────────*/ |
/*──────────────────────────────────────────────────────────────────────────────────────*/ |
||
partP: procedure expose @.; parse arg n |
partP: procedure expose @. !.; parse arg n /*obtain number (index) for computation*/ |
||
if @.n\==0 then return @.n /*Is it already computed? Return it. */ |
if @.n\==0 then return @.n /*Is it already computed? Return it. */ |
||
#= 0 /*initialize part P number.*/ |
#= 0 /*initialize part P number.*/ |
||
do k=1 for n; z= n - (k*3 - 1) * k % 2 |
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 leave /*Is Z negative? Then leave.*/ |
||
if @.z==0 then x= partP(z) /*use recursion if not known.*/ |
if @.z==0 then x= partP(z) /*use recursion if not known.*/ |
||
Line 458: | Line 1,549: | ||
else if @.z==0 then y= partP(z) /*use recursion if not known.*/ |
else if @.z==0 then y= partP(z) /*use recursion if not known.*/ |
||
else y= @.z /*use the pre─computed number*/ |
else y= @.z /*use the pre─computed number*/ |
||
if k//2 then #= # + x + y |
if k//2 then #= # + x + y /*Odd? Then sum X and Y.*/ |
||
else #= # - x |
else #= # - (x + y) /*Even? " subtract " " " */ |
||
end /*k*/ |
end /*k*/ |
||
@.n= #; |
@.n= #; return # /*define and return partitionsP of N. */</syntaxhighlight> |
||
{{out|output|text= when using the input of: <tt> 6666 </tt>}} |
{{out|output|text= when using the input of: <tt> 6666 </tt>}} |
||
<pre> |
<pre> |
||
Line 470: | Line 1,561: | ||
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 |
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 |
||
</pre> |
</pre> |
||
=== 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'''. |
|||
<syntaxhighlight 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 /*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. */</syntaxhighlight> |
|||
{{out|output|text= is identical to the 1<sup>st</sup> 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'''. |
|||
<syntaxhighlight 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 /*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. */</syntaxhighlight> |
|||
{{out|output|text= is identical to the 1<sup>st</sup> REXX version.}} <br><br> |
|||
=={{header|Rust}}== |
=={{header|Rust}}== |
||
< |
<syntaxhighlight lang="rust">// [dependencies] |
||
// rug = "1.11" |
// rug = "1.11" |
||
Line 517: | Line 1,686: | ||
println!("P({}) = {}", n, result); |
println!("P({}) = {}", n, result); |
||
println!("elapsed time: {} microseconds", time.as_micros()); |
println!("elapsed time: {} microseconds", time.as_micros()); |
||
}</ |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
Line 524: | Line 1,693: | ||
elapsed time: 8912 microseconds |
elapsed time: 8912 microseconds |
||
</pre> |
</pre> |
||
=={{header|Sidef}}== |
|||
Built-in: |
|||
<syntaxhighlight lang="ruby">say partitions(6666) # very fast</syntaxhighlight> |
|||
User-defined: |
|||
<syntaxhighlight lang="ruby">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)</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
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 |
|||
</pre> |
|||
=={{header|Swift}}== |
|||
{{trans|Rust}} |
|||
Using AttaSwift's BigInt library. |
|||
<syntaxhighlight lang="swift">import BigInt |
|||
func partitions(n: Int) -> BigInt { |
|||
var p = [BigInt(1)] |
|||
for i in 1...n { |
|||
var num = BigInt(0) |
|||
var k = 1 |
|||
while true { |
|||
var 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.append(num) |
|||
} |
|||
return p[n] |
|||
} |
|||
print("partitions(6666) = \(partitions(n: 6666))")</syntaxhighlight> |
|||
{{out}} |
|||
<pre>partitions(6666) = 193655306161707661080005073394486091998480950338405932486880600467114423441282418165863</pre> |
|||
=={{header|Wren}}== |
=={{header|Wren}}== |
||
Line 529: | Line 1,791: | ||
{{libheader|Wren-big}} |
{{libheader|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. |
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. |
||
< |
<syntaxhighlight lang="wren">import "./big" for BigInt |
||
var p = [] |
var p = [] |
||
Line 564: | Line 1,826: | ||
for (n in 2..N) partitionsP.call(n) |
for (n in 2..N) partitionsP.call(n) |
||
System.print("p[%(N)] = %(p[N])") |
System.print("p[%(N)] = %(p[N])") |
||
System.print("Took %(System.clock - start) seconds")</ |
System.print("Took %(System.clock - start) seconds")</syntaxhighlight> |
||
{{out}} |
{{out}} |
Latest revision as of 10:57, 24 January 2024
You are encouraged to solve this task according to the task description, using any language you may know.
The Partition Function P is the function P(n), where n∈ℤ, defined as the number of distinct ways in which n can be expressed as the sum of non-increasing 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 was viewed more than 100,000 times in the first couple of weeks after 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
11l
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
#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#
This can also be done using BigIntegers, but will take around five times longer. Since only adding and subtracting are required, some simple routines can be created to handle the tabulations. Also note that detecting odd and even numbers on each loop iteration is avoided by coding four increments per loop.
using System;
class Program {
const long Lm = (long)1e18;
const string Fm = "D18";
// provides 5 x 18 = 90 decimal digits
struct LI { public long lo, ml, mh, hi, tp; }
static void inc(ref LI d, LI s) { // d += s
if ((d.lo += s.lo) >= Lm) { d.ml++; d.lo -= Lm; }
if ((d.ml += s.ml) >= Lm) { d.mh++; d.ml -= Lm; }
if ((d.mh += s.mh) >= Lm) { d.hi++; d.mh -= Lm; }
if ((d.hi += s.hi) >= Lm) { d.tp++; d.hi -= Lm; }
d.tp += s.tp;
}
static void dec(ref LI d, LI s) { // d -= s
if ((d.lo -= s.lo) < 0) { d.ml--; d.lo += Lm; }
if ((d.ml -= s.ml) < 0) { d.mh--; d.ml += Lm; }
if ((d.mh -= s.mh) < 0) { d.hi--; d.mh += Lm; }
if ((d.hi -= s.hi) < 0) { d.tp--; d.hi += Lm; }
d.tp -= s.tp;
}
static LI set(long s) { LI d;
d.lo = s; d.ml = d.mh = d.hi = d.tp = 0; return d; }
static string fmt(LI x) { // returns formatted string value of x
if (x.tp > 0) return x.tp.ToString() + x.hi.ToString(Fm) + x.mh.ToString(Fm) + x.ml.ToString(Fm) + x.lo.ToString(Fm);
if (x.hi > 0) return x.hi.ToString() + x.mh.ToString(Fm) + x.ml.ToString(Fm) + x.lo.ToString(Fm);
if (x.mh > 0) return x.mh.ToString() + x.ml.ToString(Fm) + x.lo.ToString(Fm);
if (x.ml > 0) return x.ml.ToString() + x.lo.ToString(Fm);
return x.lo.ToString();
}
static LI partcount(int n) {
var P = new LI[n + 1]; P[0] = set(1);
for (int i = 1; i <= n; i++) {
int k = 0, d = -2, j = i;
LI x = set(0);
while (true) {
if ((j -= (d += 3) -k) >= 0) inc(ref x, P[j]); else break;
if ((j -= ++k) >= 0) inc(ref x, P[j]); else break;
if ((j -= (d += 3) -k) >= 0) dec(ref x, P[j]); else break;
if ((j -= ++k) >= 0) dec(ref x, P[j]); else break;
}
P[i] = x;
}
return P[n];
}
static void Main(string[] args) {
var sw = System.Diagnostics.Stopwatch.StartNew ();
var res = partcount(6666); sw.Stop();
Console.Write("{0} {1} ms", fmt(res), sw.Elapsed.TotalMilliseconds);
}
}
- Output:
193655306161707661080005073394486091998480950338405932486880600467114423441282418165863 12.9365 ms
Timing from Tio.run.
C++
GMP version
#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
Non GMP version
#include <chrono>
#include <iostream>
using namespace std;
using namespace chrono;
const long long Lm = (long)1e18;
const int Fm = 18;
struct LI { long long lo, ml, mh, hi, tp; };
LI set(long long s) { LI d;
d.lo = s; d.ml = d.mh = d.hi = d.tp = 0; return d; }
void inc(LI& d, LI s) { // d += s
if ((d.lo += s.lo) >= Lm) { d.ml++; d.lo -= Lm; }
if ((d.ml += s.ml) >= Lm) { d.mh++; d.ml -= Lm; }
if ((d.mh += s.mh) >= Lm) { d.hi++; d.mh -= Lm; }
if ((d.hi += s.hi) >= Lm) { d.tp++; d.hi -= Lm; }
d.tp += s.tp;
}
void dec(LI& d, LI s) { // d -= s
if ((d.lo -= s.lo) < 0) { d.ml--; d.lo += Lm; }
if ((d.ml -= s.ml) < 0) { d.mh--; d.ml += Lm; }
if ((d.mh -= s.mh) < 0) { d.hi--; d.mh += Lm; }
if ((d.hi -= s.hi) < 0) { d.tp--; d.hi += Lm; }
d.tp -= s.tp;
}
inline string sf(long long n) {
int len = Fm;
string result(len--, '0');
for (long long i = n; len >= 0 && i > 0; --len, i /= 10)
result[len] = '0' + i % 10;
return result;
}
string fmt(LI x) { // returns formatted string value of x
if (x.tp > 0) return to_string(x.tp) + sf(x.hi) + sf(x.mh) + sf(x.ml) + sf(x.lo);
if (x.hi > 0) return to_string(x.hi) + sf(x.mh) + sf(x.ml) + sf(x.lo);
if (x.mh > 0) return to_string(x.mh) + sf(x.ml) + sf(x.lo);
if (x.ml > 0) return to_string(x.ml) + sf(x.lo);
return to_string(x.lo);
}
LI partcount(int n) { LI p[n + 1]; p[0] = set(1);
for (int i = 1; i <= n; i++) {
int k = 0, d = -2, j = i; LI x = set(0); while (true) {
if ((j -= (d += 3) - k) >= 0) inc(x, p[j]); else break;
if ((j -= ++k) >= 0) inc(x, p[j]); else break;
if ((j -= (d += 3) - k) >= 0) dec(x, p[j]); else break;
if ((j -= ++k) >= 0) dec(x, p[j]); else break;
} p[i] = x; } return p[n]; }
int main() {
auto start = steady_clock::now();
auto result = partcount(6666);
auto end = steady_clock::now();
duration<double, milli> ms(end - start);
cout << fmt(result) << " " << ms.count() << " ms";
}
- Output:
Timing from Tio.run, but execution time can't be directly compared to the GMP version, as GMP isn't accessible at Tio.run.
193655306161707661080005073394486091998480950338405932486880600467114423441282418165863 7.32706 ms
Delphi
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
end
end
Partition.init
IO.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., 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
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
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
From the 9_billion_names_of_God_the_integer entry
' 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 n
Dim As ZString Ptr ans
Dim As Double t = Timer
ReDim big_p(max) As Mpz_ptr
For 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 buffer
While InKey <> "" : Wend
Print : Print "hit any key to end program"
Sleep
End
- Output:
PartitionsP(6666) = 193655306161707661080005073394486091998480950338405932486880600467114423441282418165863 time = 32.97 ms
Frink
Frink has a built-in function for counting partitions that uses Euler's pentagonal method. It works for arbitrarily-large integers and caches results.
println[partitionCount[6666]]
- Output:
193655306161707661080005073394486091998480950338405932486880600467114423441282418165863
Go
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
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
Haskell
{-# LANGUAGE DeriveFunctor #-}
------------------------------------------------------------
-- memoization utilities
data Memo a = Node a (Memo a) (Memo a)
deriving (Functor)
memo :: Integral a => Memo p -> a -> p
memo (Node a l r) n
| n == 0 = a
| odd n = memo l (n `div` 2)
| otherwise = memo r (n `div` 2 - 1)
nats :: Memo Int
nats =
Node
0
((+ 1) . (* 2) <$> nats)
((* 2) . (+ 1) <$> nats)
------------------------------------------------------------
-- calculating partitions
partitions :: Memo Integer
partitions = partitionP <$> nats
partitionP :: Int -> Integer
partitionP 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 :: [Int]
ofsets = scanl1 (+) $ mix [1, 3 ..] [1, 2 ..]
where
mix a b = concat $ zipWith (\x y -> [x, y]) a b
main :: IO ()
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 * ]) @ (>:@i.@>.@%:@((2%3)&*))
- Output:
pn 6 11 pn 66 2323520 pn 666 11956824258286445517629485 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]
Using gojq 0.12.11, `partitions(6666)` yields (in about 12 minutes (u+s) on a 3GHz machine):
193655306161707661080005073394486091998480950338405932486880600467114423441282418165863
The integer precision of the C implementation of jq is insufficient for computing ``partitions(6666)``, but as a test of the BigInt.jq library for jq, 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"
Curiously, the u+s time is 7m3s, which is significantly less than the above-mentioned gojq time, even though the BigInt library is written in jq.
Recursive
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 memoization
def 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 memoization
def 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+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)
Lingo
Lingo natively only supports 32 bit integers, so P(6666) would be way too big.
-- returns number of partitions of n
on partitions(n, res_table)
if n < 2 then return 1
if voidP(res_table) then
res_table = []
res_table[n] = 0
else if res_table[n] then
return res_table[n]
end if
res = 0
i = 0
param = 1
repeat while param <= n
if i mod 4 < 2 then
res = res + partitions(n - param, res_table)
else
res = res - partitions(n - param, res_table)
end if
if i mod 2 then
param = param + i + 2
else
param = param + i / 2 + 1
end if
i = i + 1
end repeat
res_table[n] = res
return res
end
- Output:
ms = _system.milliseconds put "P(121):", partitions(121) put "Time (ms):", _system.milliseconds - ms -- P(121): 2056148051 -- Time (ms): 3
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;
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
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
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
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
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 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
/* 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.
PicoLisp
Based on the Erlang implementation.
(de gpentagonals (Max)
(make
(let (N 0 M 1)
(loop
(inc 'N (if (=0 (& M 1)) (>> 1 M) M))
(T (> N Max))
(link N)
(inc 'M)))))
(de p (N)
(cache '(NIL) N
(if (=0 N)
1
(let (Sum 0 Sgn 0)
(for G (gpentagonals N)
((if (< Sgn 2) 'inc 'dec) 'Sum (p (- N G)))
(setq Sgn (& 3 (inc Sgn))))
Sum))))
- Output:
: (bench (p 6666)) 0.959 sec -> 193655306161707661080005073394486091998480950338405932486880600467114423441282418165863
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 = 1936553061617076610800050733944860919984809503384
05932486880600467114423441282418165863.
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
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)
Quackery
0 partitions
returns 1
as per oeis.org/A000041 (Partitions of n).
This is a naive recursive solution, so computing the partitions of 6666 would take a hideously long time.
[ 1 swap
dup 0 = iff drop done
[ 2dup = iff [ 2drop 1 ] done
2dup > iff [ 2drop 0 ] done
2dup dip 1+ recurse
unrot over - recurse + ] ] is partitions ( n --> n )
say "Partitions of 0 to 29" cr
30 times [ i^ partitions echo sp ]
- Output:
Partitions of 0 to 29 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 2436 3010 3718 4565
Racket
Backwards range was used to get responsive feedback for progress.
#lang racket
(require math/number-theory)
(define σ
(let ((memo (make-hash)))
(λ (z)
(hash-ref! memo z
(λ () (apply + (divisors z)))))))
(define p
(let ((memo (make-hash '((0 . 1)))))
(λ (n)
(hash-ref!
memo n
(λ ()
(let ((r (if (zero? n) 1
(/ (for/sum ((k (in-range (sub1 n) -1 -1)))
(* (σ (- n k))
(p k)))
n))))
(when (zero? (modulo n 1000)) (displayln (cons n r) (current-error-port)))
r))))))
(map p (range 1 30))
(p 666)
(p 1000)
(p 10000)
- Output:
'(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 2436 3010 3718 4565) 11956824258286445517629485 (1000 . 24061467864032622473692149727991) 24061467864032622473692149727991 (2000 . 4720819175619413888601432406799959512200344166) (3000 . 496025142797537184410324879054927095334462742231683423624) (4000 . 1024150064776551375119256307915896842122498030313150910234889093895) (5000 . 169820168825442121851975101689306431361757683049829233322203824652329144349) (6000 . 4671727531970209092971024643973690643364629153270037033856605528925072405349246129) (7000 . 32856930803440615786280925635924166861950151574532240659699032157432236394374450791229199) (8000 . 78360264351568349490593145013364599719010769352985864331118600209417827764524450990388402844164) (9000 . 77133638117808884907320791427403134961639798322072034262647713694605367979684296948790335590435626459) (10000 . 36167251325636293988820471890953695495016030339315650422081868605887952568754066420592310556052906916435144) 36167251325636293988820471890953695495016030339315650422081868605887952568754066420592310556052906916435144
Raku
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] 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
Swift
Using AttaSwift's BigInt library.
import BigInt
func partitions(n: Int) -> BigInt {
var p = [BigInt(1)]
for i in 1...n {
var num = BigInt(0)
var k = 1
while true {
var 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.append(num)
}
return p[n]
}
print("partitions(6666) = \(partitions(n: 6666))")
- Output:
partitions(6666) = 193655306161707661080005073394486091998480950338405932486880600467114423441282418165863
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.
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
- Programming Tasks
- Recursion
- Memoization
- 11l
- C
- GMP
- C sharp
- C++
- Delphi
- System.SysUtils
- Velthuis.BigIntegers
- System.Diagnostics
- Elixir
- Erlang
- F Sharp
- Factor
- FreeBASIC
- Frink
- Go
- Java
- JavaScript
- Haskell
- J
- Jq
- Julia
- Lingo
- Maple
- Mathematica
- Wolfram Language
- Nim
- Bignum
- Perl
- Phix
- Phix/mpfr
- Picat
- PicoLisp
- Prolog
- Python
- Quackery
- Racket
- Raku
- REXX
- Rust
- Sidef
- Swift
- Wren
- Wren-big