(* :Name: DiscreteMath`IntegerPartitions` *)

(* :Title: Integer Partitions. *)

(* :Author: Mark Sofroniou *)

(* :Summary:
 This package provides functions for the construction and enumeration of
 restricted and unrestricted integer partitions.
*)

(* :Context: DiscreteMath`IntegerPartitions` *)

(* :Package Version: 1.0 *)

(* :Copyright: Copyright 1995-1999, Wolfram Research, Inc. *)

(* :History:
 Original Version by Mark Sofroniou, May, 1995.
 Updated by Mark Sofroniou, January 1998. Removed some top-level
 implementations and replaced them with an interface to kernel functions.
*)

(* :Keywords:
 Partitions, Integer Partitions, Unrestricted Partitions,
 Restricted Partitions.
*)

(* :Discussion:
 Integer partitions play a key role in combinatorial number theory.
 There also exist numerous applications such as the construction of rooted
 trees (and hence Frechet derivatives) via Faa di Bruno's formula.

 Ferrers diagrams are a useful way of representing partitions pictorially.
 Some interesting properties include the relationship between transpositions
 of partitions, where a partition of length k is in one-to-one correspondence
 with a partition of size k.

 The multiset representation (a set with multiplicities) is the key
 to efficient running time.
*)

(* :Source:
 A. Nijenhuis, H. S. Wilf, Combinatorial Algorithms. Second Edition.
 Academic Press, New York, 1978.
*)

(* :Mathematica Version: 4.0 *)

(* :Limitations:
 The limiting factor here is primarily memory consumption. For example there
 are 190569292 partitions of 100. *)

BeginPackage["DiscreteMath`IntegerPartitions`","Utilities`FilterOptions`"]

IntegerPartitions::usage =
"IntegerPartitions[n] constructs the unrestricted partitions of the
integer n in reverse lexicographic order. The result is a list of lists of
integers i grouped together with their multiplicities m as {i,m} such that
the sum of i*m is n. IntegerPartitions[n,m] gives the partitions of n with
size at most m."

NextIntegerPartition::usage =
"NextIntegerPartition[p] returns the integer partition following p in reverse
lexicographic order."

PreviousIntegerPartition::usage =
"PreviousIntegerPartition[p] returns the integer partition preceding p in
reverse lexicographic order."

ConstrainedIntegerPartitionsP::usage =
"ConstrainedIntegerPartitionsP[n,k] gives the number of integer partitions
of n into k parts."

FerrersPlot::usage =
"FerrersPlot[p, opts] plots the Ferrers diagram of the integer partitions p
of n, where opts can be FerrersPlotColumns, FerrersPlotLabel,
FerrersPlotPointSize or any Graphics option."

FerrersPlotColumns::usage =
"FerrersPlotColumns is an option to FerrersPlot
that specifies the number of columns in the GraphicsArray plot of the
partitions of an integer."

FerrersPlotLabel::usage =
"FerrersPlotLabel is an option to FerrersPlot that
specifies a list of plot labels, which should be the same length as the
list of integer partitions to be plotted."

FerrersPlotPointSize::usage =
"FerrersPlotPointSize is an option to FerrersPlot
that specifies a scaling factor to be applied to the default size for
point data in the Ferrers plot."

IntegerPartitionQ::usage =
"IntegerPartitionQ[p] returns True if p is an integer partition and False
otherwise."

MultiplicityFromRepetition::usage =
"MultiplicityFromRepetition[p] transforms the integer partition p from
a representation with repeated elements i,i,...,i to a multiset representation
with multiplicites {i,m}."

RandomIntegerPartition::usage =
"RandomIntegerPartition[n] constructs a random partition of the integer n."

RepetitionFromMultiplicity::usage =
"RepetitionFromMultiplicity[p] transforms the integer partition p from
a multiset representation with multiplicites {i,m} to a repeated representation
i,i,...,i."

RestrictedIntegerPartitionsP::usage =
"RestrictedIntegerPartitionsP[n,k] gives the total number of restricted
integer partitions of n into parts of size at most k. This is the sum of the
integer partitions constrained to j parts with 1<=j<=k."

TransposeIntegerPartition::usage =
"TransposeIntegerPartition[p] creates a partition with maximum part k
from a partition p with k parts."

Unprotect[ ConstrainedIntegerPartitionsP, FerrersPlot,
FerrersPlotColumns, FerrersPlotLabel, FerrersPlotPointSize,
IntegerPartitions, IntegerPartitionQ, MultiplicityFromRepetition,
NextIntegerPartition, PreviousIntegerPartition,
RepetitionFromMultiplicity, RandomIntegerPartition,
RestrictedIntegerPartitionsP, TransposeIntegerPartition];

Begin["`Private`"];

(**** Interfaces to kernel code. ****)

IntegerPartitions[_Integer?Positive,0]:= {};
IntegerPartitions[n_Integer,m_Integer]:=
  System`Private`IntegerPartitions[n,m] /; 1<=m<=n;
IntegerPartitions[n_Integer?Positive]:=
  System`Private`IntegerPartitions[ n ];

(* Transform to next adjacent partition. *)

NextIntegerPartition[p_?IntegerPartitionQ]:=
   System`Private`NextIntegerPartition[p];

(* Transform to previous adjacent partition. *)

PreviousIntegerPartition[p_?IntegerPartitionQ]:=
  System`Private`PreviousIntegerPartition[p];


(**** Syntax test for a valid partition. ****)

IntegerPartitionQ[p_]:=
  Module[{res},
    res = System`Private`IntegerPartitionQ[p];
	res /; Head[res]=!=System`Private`IntegerPartitionQ
  ];


(**** Transpose partition ****)

TransposeIntegerPartition[p_?IntegerPartitionQ]:=
  System`Private`TransposeIntegerPartition[p];

(* Total number of restricted partitions of n into parts of size
 at most k. *)

(* Number of restricted partitions *)

RestrictedIntegerPartitionsP[n_Integer, k_Integer]:=
  Block[{$RecursionLimit=Infinity}, rip[n,k] ] /; 0<=k<=n;
RestrictedIntegerPartitionsP[_Integer, _Integer] = 0;

rip[n_,k_]:= rip[n,k] = Sum[cip[n,i],{i,k}];


(**** Number of partitions of n into k parts ****)

ConstrainedIntegerPartitionsP[n_Integer, k_Integer]:=
  Block[{$RecursionLimit=Infinity}, cip[n,k] ] /; 0<=k<=n;
ConstrainedIntegerPartitionsP[_Integer,_Integer] = 0;

cip[n_,k_]:= 0 /; n<k;
cip[_,1] = 1;
cip[n_,k_]:= cip[n,k]= cip[n-k,k] + cip[n-1,k-1];


(**** Random Partition ****)

RandomIntegerPartition[(n_Integer)?Positive] := 
  Module[{mult = Table[0,{n}], j, d, r = n, z}, 
    While[r > 0, d = 1; j = 0; 
      z = Random[]*r*PartitionsP[r]; 
      While[z >= 0,
        j++;
        If[r - j*d < 0, {j = 1; d++; }];
        z -= j*PartitionsP[r - j*d]; ];
        r -= j*d;
        mult[[j]] += d;
      ];
      Table[randomElement[Part[mult,j],j], {j,n,1,-1}]
  ];

randomElement[0,_]:= Sequence[];
randomElement[m_,i_]:= {i,m};


(**** Plot Ferrers diagram. ****)

Options[FerrersPlot] =
  {FerrersPlotColumns -> Automatic, FerrersPlotLabel->{},
     FerrersPlotPointSize -> 1};

FerrersPlot::opts =
"The option `1` in FerrersPlot did not evaluate to `2`"

FPmessages = {
  {"FerrersPlotColumns","a positive integer or Automatic"},
  {"FerrersPlotLabel","a list of plot labels"},
  {"FerrersPlotPointSize","a positive number"}};

FPtest =
  If[#1,
    True,
    Message[FerrersPlot::opts,Apply[Sequence,#2]]; False
  ]&;

FPOptionsTest[lenp_,opts___]:=
  Module[{col, lbl, ps, datatypes},

    {col,lbl,ps} =
      {FerrersPlotColumns, FerrersPlotLabel, FerrersPlotPointSize} /.
        {opts} /. Options[FerrersPlot];

    lbl = Flatten[lbl];

    datatypes = {
      MatchQ[col, Automatic | (_Integer?Positive)],
      MatchQ[lbl,{}|_?(VectorQ[#]&&Length[#]===Length[lenp]&)],
      TrueQ[Positive[ps]]};

    If[Apply[And, MapThread[FPtest, {datatypes, FPmessages}]],
      {col,lbl,ps},
      $Failed
    ]
  ];

FerrersPlot[p_?IntegerPartitionQ, opts___?OptionQ]:=
  Module[{ans},
    ans /; ((ans = FPOptionsTest[1,opts]) =!= $Failed &&
              (ans = fplot[{p}, ans, FilterOptions[Graphics,opts]];
                True))
  ];

FerrersPlot[p:{__?IntegerPartitionQ}, opts___?OptionQ]:=
  Module[{ans, lenp = Length[p]},
    ans /; ((ans = FPOptionsTest[lenp,opts]) =!= $Failed &&
              (ans = fplot[p, ans, FilterOptions[Graphics,opts]];
                True))
  ];

(* Construct the expanded form of a partition from the multiset
 representation. This converts to the form used in Combinatorica.m. *)

RepetitionFromMultiplicity[p_?IntegerPartitionQ]:=
  System`Private`RepetitionFromMultiplicity[p];

(* Construct the multiset form of a partition from the expanded
 representation. This converts from the form used in Combinatorica.m. *)

MultiplicityFromRepetition[p:{n__Integer?Positive}]:=
  mfr[p] /; GreaterEqual[n];

mfr[p_]:= Map[ {First[#], Length[#]}&, Split[p] ];

makeArray[p_, n_, m_]:=
  Module[{pp, rem},
    pp = Partition[p, m];
    rem = Mod[n,m];
    If[rem!=0, AppendTo[pp, Take[p, -rem]]];
    pp
  ];

fplot[p_, {col_, lbl_, ps_},gropts___]:=
  Module[{lenp, maxp, n, pr, prims, ptsz, scale, safety},

    lenp = Length[p];
    safety = 1.0; (* Border around each Graphics cell *)
    scale = 0.05 ps; (* Scale for point size *)

(* Find the maximum number that was used to generate any partition. This is used
 to scale the elements of the GraphicsArray. *)

    maxp = Max[ Apply[Plus,Apply[Times,p,{2}],{1}] ];

    If[lenp==1, (* Special case for one partition *)
      ptsz = scale (maxp)^-0.25;
      pr = All, (* PlotRange->All for one partition *)
      ptsz = scale (maxp)^0.25; (* Scales as Sqrt[Sqrt[maxp]] *)
      pr = {{1-safety, maxp+safety},{1-safety, maxp+safety}};
    ];

    n = If[col === Automatic, Ceiling[Sqrt[lenp]], col];

    prims =
      MapThread[
        fpprimitive[Reverse[rfm[#1]], pr, ptsz, #2]&,
        {p, If[lbl==={}, Table[None,{lenp}], lbl]}
      ];

    Show[
      GraphicsArray[makeArray[prims, lenp, n]], gropts, PlotRange->All
    ]
  ];

fpprimitive[p_,pr_,ptsz_,lbl_]:=
  Graphics[
    Flatten[{
      PointSize[ ptsz ],
      Table[Point[{i,j}], {j,Length[p]}, {i,Part[p,j]}]
    }], PlotLabel -> lbl, AspectRatio -> 1, PlotRange -> pr
  ];


End[];

SetAttributes[{ConstrainedIntegerPartitionsP, FerrersPlot,
IntegerPartitions, IntegerPartitionQ, MultiplicityFromRepetition,
NextIntegerPartition, PreviousIntegerPartition,
RepetitionFromMultiplicity, RandomIntegerPartition,
RestrictedIntegerPartitionsP, TransposeIntegerPartition},
{ReadProtected}];

SetAttributes[{ConstrainedIntegerPartitionsP, IntegerPartitions,
RandomIntegerPartition, RestrictedIntegerPartitionsP},
Listable];

Protect[ConstrainedIntegerPartitionsP, FerrersPlot,
FerrersPlotColumns, FerrersPlotLabel, FerrersPlotPointSize,
IntegerPartitions, IntegerPartitionQ, MultiplicityFromRepetition,
NextIntegerPartition, PreviousIntegerPartition,
RepetitionFromMultiplicity, RandomIntegerPartition,
RestrictedIntegerPartitionsP, TransposeIntegerPartition];

EndPackage[];

