home *** CD-ROM | disk | FTP | other *** search
/ Math Solutions 1995 October / Math_Solutions_CD-ROM_Walnut_Creek_October_1995.iso / pc / mac / discrete / lib / combinat.g < prev    next >
Encoding:
Text File  |  1993-05-05  |  48.8 KB  |  1,564 lines
  1. #############################################################################
  2. ##
  3. #A  combinat.g                  GAP library                  Martin Schoenert
  4. ##
  5. #A  @(#)$Id: combinat.g,v 3.21 1992/07/07 13:30:02 goetz Rel $
  6. ##
  7. #Y  Copyright 1990-1992,  Lehrstuhl D fuer Mathematik,  RWTH Aachen,  Germany
  8. ##
  9. ##  This file contains  the functions that mainly  deal  with  combinatorics.
  10. ##
  11. #H  $Log: combinat.g,v $
  12. #H  Revision 3.21  1992/07/07  13:30:02  goetz
  13. #H  made 'PowerPartition' faster.
  14. #H
  15. #H  Revision 3.20  1992/06/01  15:13:21  goetz
  16. #H  a more readable version of 'PartitionTuples' which accepts 0.
  17. #H
  18. #H  Revision 3.19  1992/03/13  13:19:49  goetz
  19. #H  added 'PartitionTuples' from 'ctsymmet.g'.
  20. #H
  21. #H  Revision 3.18  1992/03/13  13:13:20  goetz
  22. #H  added 'SignPartition', 'AssociatedPartition' and
  23. #H  'PowerPartition' from 'ctsymmet.g'.
  24. #H
  25. #H  Revision 3.17  1991/10/08  16:10:50  martin
  26. #H  'Position' now returns 'false'
  27. #H
  28. #H  Revision 3.16  1991/08/08  16:12:33  martin
  29. #H  fixed functions that used 'Size' instead of 'Length' for multisets
  30. #H
  31. #H  Revision 3.15  1991/07/30  15:13:58  martin
  32. #H  improved 'NrPartitions' to use Eulers theorem
  33. #H
  34. #H  Revision 3.14  1991/07/22  17:40:24  martin
  35. #H  added 'RestrictedPartitions'
  36. #H
  37. #H  Revision 3.13  1991/07/21  17:56:24  martin
  38. #H  added lots of comments
  39. #H
  40. #H  Revision 3.12  1991/07/21  13:04:14  martin
  41. #H  removed avoidable 'return'-s
  42. #H
  43. #H  Revision 3.11  1991/07/21  12:00:56  martin
  44. #H  changed 'Partitions' to return partitions in standard form
  45. #H
  46. #H  Revision 3.10  1991/07/21  11:39:38  martin
  47. #H  fixed 'Partitions' for 0
  48. #H
  49. #H  Revision 3.9  1991/07/10  12:00:00  martin
  50. #H  fixed 'PartitionsSet' for the empty set
  51. #H
  52. #H  Revision 3.8  1991/07/09  12:00:00  martin
  53. #H  fixed 'OrderedPartitions' for the empty set
  54. #H
  55. #H  Revision 3.7  1991/07/05  12:00:00  martin
  56. #H  fixed 'NrArrangements' for the empty set
  57. #H
  58. #H  Revision 3.6  1991/07/01  13:00:00  martin
  59. #H  added 'Bell'
  60. #H
  61. #H  Revision 3.5  1991/07/01  12:00:00  martin
  62. #H  fixed 'NrDerangements' for the empty set
  63. #H
  64. #H  Revision 3.4  1990/11/20  13:00:00  martin
  65. #H  added 'Bernoulli'
  66. #H
  67. #H  Revision 3.3  1990/11/20  12:00:00  martin
  68. #H  moved functions from numtheor to combinatorics package
  69. #H
  70. #H  Revision 3.2  1990/11/16  12:00:00  martin
  71. #H  moved functions from integer to combinatorics package
  72. #H
  73. #H  Revision 3.1  1990/11/11  12:00:00  martin
  74. #H  added 'Permanent'
  75. #H
  76. #H  Revision 3.0  1990/11/09  12:00:00  martin
  77. #H  initial revision under RCS
  78. #H
  79. ##
  80.  
  81.  
  82. #############################################################################
  83. ##
  84. #F  Factorial( <n> )  . . . . . . . . . . . . . . . . factorial of an integer
  85. ##
  86. Factorial := function ( n )
  87.     if n < 0  then Error("<n> must be nonnegative");  fi;
  88.     return Product( [1..n] );
  89. end;
  90.  
  91.  
  92. #############################################################################
  93. ##
  94. #F  Binomial( <n>, <k> )  . . . . . . . . .  binomial coefficient of integers
  95. ##
  96. Binomial := function ( n, k )
  97.     local   bin, i, j;
  98.     if   k < 0  then
  99.         bin := 0;
  100.     elif k = 0  then
  101.         bin := 1;
  102.     elif n < 0  then
  103.         bin := (-1)^k * Binomial( -n+k-1, k );
  104.     elif n < k  then
  105.         bin := 0;
  106.     elif n = k  then
  107.         bin := 1;
  108.     elif n-k < k  then
  109.         bin := Binomial( n, n-k );
  110.     else
  111.         bin := 1;  j := 1;
  112.         for i  in [n-k+1..n]  do
  113.             bin := bin * i;
  114.             while j <= k  and bin mod j = 0  do
  115.                 bin := bin / j;
  116.                 j := j + 1;
  117.             od;
  118.         od;
  119.     fi;
  120.     return bin;
  121. end;
  122.  
  123.  
  124. #############################################################################
  125. ##
  126. #F  Bell( <n> ) . . . . . . . . . . . . . . . . .  value of the Bell sequence
  127. ##
  128. Bell := function ( n )
  129.     local   bell, k, i;
  130.     bell := [ 1 ];
  131.     for i  in [1..n-1]  do
  132.         bell[i+1] := bell[1];
  133.         for k  in [0..i-1]  do
  134.             bell[i-k] := bell[i-k] + bell[i-k+1];
  135.         od;
  136.     od;
  137.     return bell[1];
  138. end;
  139.  
  140.  
  141. #############################################################################
  142. ##
  143. #F  Stirling1( <n>, <k> ) . . . . . . . . . Stirling number of the first kind
  144. ##
  145. Stirling1 := function ( n, k )
  146.     local   sti, i, j;
  147.     if   n < k  then
  148.         sti := 0;
  149.     elif n = k  then
  150.         sti := 1;
  151.     elif n < 0 and k < 0  then
  152.         sti := Stirling2( -k, -n );
  153.     elif k <= 0  then
  154.         sti := 0;
  155.     else
  156.         sti := [ 1 ];
  157.         for j  in [2..n-k+1]  do
  158.             sti[j] := 0;
  159.         od;
  160.         for i  in [1..k]  do
  161.             sti[1] := 1;
  162.             for j  in [2..n-k+1]  do
  163.                 sti[j] := (i+j-2) * sti[j-1] + sti[j];
  164.             od;
  165.         od;
  166.         sti := sti[n-k+1];
  167.     fi;
  168.     return sti;
  169. end;
  170.  
  171.  
  172. #############################################################################
  173. ##
  174. #F  Stirling2( <n>, <k> ) . . . . . . . .  Stirling number of the second kind
  175. ##
  176. ##  Uses $S_2(n,k) = (-1)^k \sum_{i=1}^{k}{(-1)^i {k \choose i} i^k} / k!$.
  177. ##
  178. Stirling2 := function ( n, k )
  179.     local   sti, bin, fib, i;
  180.     if   n < k  then
  181.         sti := 0;
  182.     elif n = k  then
  183.         sti := 1;
  184.     elif n < 0 and k < 0  then
  185.         sti := Stirling1( -k, -n );
  186.     elif k <= 0  then
  187.         sti := 0;
  188.     else
  189.         bin := 1;                       # (k 0)
  190.         sti := 0;                       # (-1)^0 (k 0) 0^k
  191.         fib := 1;                       # 0!
  192.         for i  in [1..k]  do
  193.             bin := (k-i+1)/i * bin;     # (k i) = (k-(i-1))/i (k i-1)
  194.             sti := bin * i^n - sti;     # (-1)^i sum (-1)^j (k j) j^k
  195.             fib := fib * i;             # i!
  196.         od;
  197.         sti := sti / fib;
  198.     fi;
  199.     return sti;
  200. end;
  201.  
  202.  
  203. #############################################################################
  204. ##
  205. #F  Combinations( <mset>, <k> ) . . . .  set of sorted sublists of a multiset
  206. ##
  207. ##  'CombinationsA( <mset>, <m>,  <n>, <comb>, <i> )' returns  the set of all
  208. ##  combinations of the multiset <mset>, which has size  <n>, that begin with
  209. ##  '<comb>[[1..<i>-1]]'.  To do  this it finds  all elements of <mset>  that
  210. ##  can go at '<comb>[<i>]' and calls itself  recursively for each candidate.
  211. ##  <m>-1 is the position of '<comb>[<i>-1]' in <mset>, so the candidates for
  212. ##  '<comb>[<i>]' are exactely the elements 'Set( <mset>[[<m>..<n>]] )'.
  213. ##
  214. ##  'CombinationsK( <mset>, <m>, <n>, <k>, <comb>, <i>  )' returns the set of
  215. ##  all combinations  of the multiset <mset>,  which has size  <n>, that have
  216. ##  length '<i>+<k>-1', and that begin with '<comb>[[1..<i>-1]]'.  To do this
  217. ##  it finds  all elements of  <mset> that can go  at '<comb>[<i>]' and calls
  218. ##  itself recursively  for   each candidate.    <m>-1 is   the  position  of
  219. ##  '<comb>[<i>-1]'  in <mset>,  so  the  candidates  for '<comb>[<i>]'   are
  220. ##  exactely the elements 'Set( <mset>[<m>..<n>-<k>+1] )'.
  221. ##
  222. ##  'Combinations' only calls 'CombinationsA' or 'CombinationsK' with initial
  223. ##  arguments.
  224. ##
  225. CombinationsA := function ( mset, m, n, comb, i )
  226.     local   combs, l;
  227.     if m = n+1  then
  228.         comb  := ShallowCopy(comb);
  229.         combs := [ comb ];
  230.     else
  231.         comb  := ShallowCopy(comb);
  232.         combs := [ ShallowCopy(comb) ];
  233.         for l  in [m..n]  do
  234.             if l = m or mset[l] <> mset[l-1]  then
  235.                 comb[i] := mset[l];
  236.                 Append( combs, CombinationsA(mset,l+1,n,comb,i+1) );
  237.             fi;
  238.         od;
  239.     fi;
  240.     return combs;
  241. end;
  242.  
  243. CombinationsK := function ( mset, m, n, k, comb, i )
  244.     local   combs, l;
  245.     if k = 0  then
  246.         comb  := ShallowCopy(comb);
  247.         combs := [ comb ];
  248.     else
  249.         combs := [];
  250.         for l  in [m..n-k+1]  do
  251.             if l = m or mset[l] <> mset[l-1]  then
  252.                 comb[i] := mset[l];
  253.                 Append( combs, CombinationsK(mset,l+1,n,k-1,comb,i+1) );
  254.             fi;
  255.         od;
  256.     fi;
  257.     return combs;
  258. end;
  259.  
  260. Combinations := function ( arg )
  261.     local   combs, mset;
  262.     if Length(arg) = 1  then
  263.         mset := ShallowCopy(arg[1]);  Sort( mset );
  264.         combs := CombinationsA( mset, 1, Length(mset), [], 1 );
  265.     elif Length(arg) = 2  then
  266.         mset := ShallowCopy(arg[1]);  Sort( mset );
  267.         combs := CombinationsK( mset, 1, Length(mset), arg[2], [], 1 );
  268.     else
  269.         Error("usage: Combinations( <mset> [, <k>] )");
  270.     fi;
  271.     return combs;
  272. end;
  273.  
  274.  
  275. #############################################################################
  276. ##
  277. #F  NrCombinations( <mset>, <k> ) . . number of sorted sublists of a multiset
  278. ##
  279. ##  'NrCombinations' uses well known identities if <mset> is a proper set, or
  280. ##  if no <k> is given.  For a multiset it uses 'NrCombinationsK' which works
  281. ##  just like 'CombinationsK'.
  282. ##
  283. #N  26-Oct-90 martin there should be a better way to do this for multisets
  284. ##
  285. NrCombinationsK := function ( mset, m, n, k )
  286.     local   combs, l;
  287.     if k = 0  then
  288.         combs := 1;
  289.     else
  290.         combs := 0;
  291.         for l  in [m..n-k+1]  do
  292.             if l = m  or mset[l] <> mset[l-1]  then
  293.                 combs := combs + NrCombinationsK( mset, l+1, n, k-1 );
  294.             fi;
  295.         od;
  296.     fi;
  297.     return combs;
  298. end;
  299.  
  300. NrCombinations := function ( arg )
  301.     local   nr, mset, i, j;
  302.     if Length(arg) = 1  then
  303.         mset := ShallowCopy(arg[1]);  Sort( mset );
  304.         if IsSet(mset)  then
  305.             nr := 2^Size(mset);
  306.         else
  307.             nr := Product( Set(mset), i->Number(mset,j->i=j)+1 );
  308.         fi;
  309.     elif Length(arg) = 2  then
  310.         mset := ShallowCopy(arg[1]);  Sort( mset );
  311.         if IsSet(mset)  then
  312.             nr := Binomial( Size(mset), arg[2] );
  313.         else
  314.             nr := NrCombinationsK( mset, 1, Length(mset), arg[2] );
  315.         fi;
  316.     else
  317.         Error("usage: NrCombinations( <mset> [, <k>] )");
  318.     fi;
  319.     return nr;
  320. end;
  321.  
  322.  
  323. #############################################################################
  324. ##
  325. #F  Arrangements( <mset> )  . . . . set of ordered combinations of a multiset
  326. ##
  327. ##  'ArrangementsA( <mset>,  <m>, <n>, <comb>, <i>  )' returns the set of all
  328. ##  arrangements of the multiset <mset>, which has  size <n>, that begin with
  329. ##  '<comb>[[1..<i>-1]]'.   To do this it  finds all  elements of <mset> that
  330. ##  can go at '<comb>[<i>]' and calls itself  recursively for each candidate.
  331. ##  <m> is a boolean list of size <n> that contains  'true' for every element
  332. ##  of <mset> that we have not yet taken, so the candidates for '<comb>[<i>]'
  333. ##  are exactely the elements '<mset>[<l>]' such that '<m>[<l>]'  is  'true'.
  334. ##  Some care must be taken to take a candidate only once if it appears  more
  335. ##  than once in <mset>.
  336. ##
  337. ##  'ArrangementsK( <mset>, <m>, <n>, <k>, <comb>, <i> )'  returns the set of
  338. ##  all arrangements  of the multiset <mset>,  which has size  <n>, that have
  339. ##  length '<i>+<k>-1', and that begin with '<comb>[[1..<i>-1]]'.  To do this
  340. ##  it finds  all elements of  <mset> that can  go at '<comb>[<i>]' and calls
  341. ##  itself recursively for each candidate.  <m> is a boolean list of size <n>
  342. ##  that contains 'true' for every element  of <mset>  that we  have  not yet
  343. ##  taken,  so  the candidates for   '<comb>[<i>]' are  exactely the elements
  344. ##  '<mset>[<l>]' such that '<m>[<l>]' is 'true'.  Some care must be taken to
  345. ##  take a candidate only once if it appears more than once in <mset>.
  346. ##
  347. ##  'Arrangements' only calls 'ArrangementsA' or 'ArrangementsK' with initial
  348. ##  arguments.
  349. ##
  350. ArrangementsA := function ( mset, m, n, comb, i )
  351.     local   combs, l;
  352.     if i = n+1  then
  353.         comb  := ShallowCopy(comb);
  354.         combs := [ comb ];
  355.     else
  356.         comb  := ShallowCopy(comb);
  357.         combs := [ ShallowCopy(comb) ];
  358.         for l  in [1..n]  do
  359.             if m[l] and (l=1 or m[l-1]=false or mset[l]<>mset[l-1])  then
  360.                 comb[i] := mset[l];
  361.                 m[l] := false;
  362.                 Append( combs, ArrangementsA( mset, m, n, comb, i+1 ) );
  363.                 m[l] := true;
  364.             fi;
  365.         od;
  366.     fi;
  367.     return combs;
  368. end;
  369.  
  370. ArrangementsK := function ( mset, m, n, k, comb, i )
  371.     local   combs, l;
  372.     if k = 0  then
  373.         comb := ShallowCopy(comb);
  374.         combs := [ comb ];
  375.     else
  376.         combs := [];
  377.         for l  in [1..n]  do
  378.             if m[l] and (l=1 or m[l-1]=false or mset[l]<>mset[l-1])  then
  379.                 comb[i] := mset[l];
  380.                 m[l] := false;
  381.                 Append( combs, ArrangementsK( mset, m, n, k-1, comb, i+1 ) );
  382.                 m[l] := true;
  383.             fi;
  384.         od;
  385.     fi;
  386.     return combs;
  387. end;
  388.  
  389. Arrangements := function ( arg )
  390.     local   combs, mset, m, i;
  391.     if Length(arg) = 1  then
  392.         mset := ShallowCopy(arg[1]);  Sort( mset );
  393.         m := List( mset, i->true );
  394.         combs := ArrangementsA( mset, m, Length(mset), [], 1 );
  395.     elif Length(arg) = 2  then
  396.         mset := ShallowCopy(arg[1]);  Sort( mset );
  397.         m := List( mset, i->true );
  398.         combs := ArrangementsK( mset, m, Length(mset), arg[2], [], 1 );
  399.     else
  400.         Error("usage: Arrangements( <mset> [, <k>] )");
  401.     fi;
  402.     return combs;
  403. end;
  404.  
  405.  
  406. #############################################################################
  407. ##
  408. #F  NrArrangements( <mset> )  .  number of ordered combinations of a multiset
  409. ##
  410. ##  'NrArrangements' uses well known identities if <mset>  is a  proper  set.
  411. ##  If <mset> is a multiset  it uses 'NrArrangementsA' and 'NrArrangementsK',
  412. ##  which work just like 'ArrangementsA' and 'ArrangementsK'.
  413. ##
  414. #N  26-Oct-90 martin there should be a better way to do this for multisets
  415. ##
  416. NrArrangementsA := function ( mset, m, n, i )
  417.     local   combs, l;
  418.     if i = n+1  then
  419.         combs := 1;
  420.     else
  421.         combs := 1;
  422.         for l  in [1..n]  do
  423.             if m[l] and (l=1 or m[l-1]=false or mset[l]<>mset[l-1])  then
  424.                 m[l] := false;
  425.                 combs := combs + NrArrangementsA( mset, m, n, i+1 );
  426.                 m[l] := true;
  427.             fi;
  428.         od;
  429.     fi;
  430.     return combs;
  431. end;
  432.  
  433. NrArrangementsK := function ( mset, m, n, k )
  434.     local   combs, l;
  435.     if k = 0  then
  436.         combs := 1;
  437.     else
  438.         combs := 0;
  439.         for l  in [1..n]  do
  440.             if m[l] and (l=1 or m[l-1]=false or mset[l]<>mset[l-1])  then
  441.                 m[l] := false;
  442.                 combs := combs + NrArrangementsK( mset, m, n, k-1 );
  443.                 m[l] := true;
  444.             fi;
  445.         od;
  446.     fi;
  447.     return combs;
  448. end;
  449.  
  450. NrArrangements := function ( arg )
  451.     local   nr, i, mset, m;
  452.     if Length(arg) = 1  then
  453.         mset := ShallowCopy(arg[1]);  Sort( mset );
  454.         if IsSet( mset )  then
  455.             nr := 0;
  456.             for i  in [0..Size(mset)]  do
  457.                 nr := nr + Factorial(Size(mset)) / Factorial(Size(mset)-i);
  458.             od;
  459.         else
  460.             m := List( mset, i->true );
  461.             nr := NrArrangementsA( mset, m, Length(mset), 1 );
  462.         fi;
  463.     elif Length(arg) = 2  then
  464.         mset := ShallowCopy(arg[1]);  Sort( mset );
  465.         if IsSet( mset )  then
  466.             if arg[2] <= Size(mset)  then
  467.                 nr := Factorial(Size(mset)) / Factorial(Size(mset)-arg[2]);
  468.             else
  469.                 nr := 0;
  470.             fi;
  471.         else
  472.             m := List( mset, i->true );
  473.             nr := NrArrangementsK( mset, m, Length(mset), arg[2] );
  474.         fi;
  475.     else
  476.         Error("usage: NrArrangements( <mset> [, <k>] )");
  477.     fi;
  478.     return nr;
  479. end;
  480.  
  481.  
  482. #############################################################################
  483. ##
  484. #F  UnorderedTuples( <set>, <k> ) . . . .  set of unordered tuples from a set
  485. ##
  486. ##  'UnorderedTuplesK( <set>, <n>, <m>, <k>, <tup>, <i> )' returns the set of
  487. ##  all unordered tuples of  the  set <set>,  which  has size <n>, that  have
  488. ##  length '<i>+<k>-1', and that begin with  '<tup>[[1..<i>-1]]'.  To do this
  489. ##  it  finds all elements  of <set>  that  can go  at '<tup>[<i>]' and calls
  490. ##  itself   recursively  for   each  candidate.  <m>  is    the  position of
  491. ##  '<tup>[<i>-1]' in <set>, so the  candidates for '<tup>[<i>]' are exactely
  492. ##  the elements '<set>[[<m>..<n>]]',  since we require that unordered tuples
  493. ##  be sorted.
  494. ##
  495. ##  'UnorderedTuples' only calls 'UnorderedTuplesK' with initial arguments.
  496. ##
  497. UnorderedTuplesK := function ( set, n, m, k, tup, i )
  498.     local   tups, l;
  499.     if k = 0  then
  500.         tup := ShallowCopy(tup);
  501.         tups := [ tup ];
  502.     else
  503.         tups := [];
  504.         for l  in [m..n]  do
  505.             tup[i] := set[l];
  506.             Append( tups, UnorderedTuplesK( set, n, l, k-1, tup, i+1 ) );
  507.         od;
  508.     fi;
  509.     return tups;
  510. end;
  511.  
  512. UnorderedTuples := function ( set, k )
  513.     set := Set(set);
  514.     return UnorderedTuplesK( set, Size(set), 1, k, [], 1 );
  515. end;
  516.  
  517.  
  518. #############################################################################
  519. ##
  520. #F  NrUnorderedTuples( <set>, <k> ) . . number unordered of tuples from a set
  521. ##
  522. NrUnorderedTuples := function ( set, k )
  523.     return Binomial( Size(Set(set))+k-1, k );
  524. end;
  525.  
  526.  
  527. #############################################################################
  528. ##
  529. #F  Tuples( <set>, <k> )  . . . . . . . . .  set of ordered tuples from a set
  530. ##
  531. ##  'TuplesK( <set>, <k>, <tup>, <i> )' returns the set  of all tuples of the
  532. ##  set   <set>   that have   length     '<i>+<k>-1',  and  that  begin  with
  533. ##  '<tup>[[1..<i>-1]]'.  To do this  it loops  over  all elements of  <set>,
  534. ##  puts them at '<tup>[<i>]' and calls itself recursively.
  535. ##
  536. ##  'Tuples' only calls 'TuplesK' with initial arguments.
  537. ##
  538. TuplesK := function ( set, k, tup, i )
  539.     local   tups, l;
  540.     if k = 0  then
  541.         tup := ShallowCopy(tup);
  542.         tups := [ tup ];
  543.     else
  544.         tups := [];
  545.         for l  in set  do
  546.             tup[i] := l;
  547.             Append( tups, TuplesK( set, k-1, tup, i+1 ) );
  548.         od;
  549.     fi;
  550.     return tups;
  551. end;
  552.  
  553. Tuples := function ( set, k )
  554.     set := Set(set);
  555.     return TuplesK( set, k, [], 1 );
  556. end;
  557.  
  558.  
  559. #############################################################################
  560. ##
  561. #F  NrTuples( <set>, <k> )  . . . . . . . number of ordered tuples from a set
  562. ##
  563. NrTuples := function ( set, k )
  564.     return Size(Set(set)) ^ k;
  565. end;
  566.  
  567.  
  568. #############################################################################
  569. ##
  570. #F  PermutationsList( <mset> )  . . . . . . set of permutations of a multiset
  571. ##
  572. ##  'PermutationsListK( <mset>, <m>, <n>, <k>, <perm>, <i> )' returns the set
  573. ##  of all  permutations  of the multiset  <mset>, which  has  size <n>, that
  574. ##  begin  with '<perm>[[1..<i>-1]]'.  To do   this it finds all elements  of
  575. ##  <mset> that can go at '<perm>[<i>]' and calls itself recursively for each
  576. ##  candidate.  <m> is a  boolean  list of size  <n> that contains 'true' for
  577. ##  every element of <mset> that we have not yet taken, so the candidates for
  578. ##  '<perm>[<i>]'  are    exactely the   elements   '<mset>[<l>]' such   that
  579. ##  '<m>[<l>]' is 'true'.  Some care must be taken to take a  candidate  only
  580. ##  once if it apears more than once in <mset>.
  581. ##
  582. ##  'Permutations' only calls 'PermutationsListK' with initial arguments.
  583. ##
  584. PermutationsListK := function ( mset, m, n, k, perm, i )
  585.     local   perms, l;
  586.     if k = 0  then
  587.         perm := ShallowCopy(perm);
  588.         perms := [ perm ];
  589.     else
  590.         perms := [];
  591.         for l  in [1..n]  do
  592.             if m[l] and (l=1 or m[l-1]=false or mset[l]<>mset[l-1])  then
  593.                 perm[i] := mset[l];
  594.                 m[l] := false;
  595.                 Append( perms, PermutationsListK(mset,m,n,k-1,perm,i+1) );
  596.                 m[l] := true;
  597.             fi;
  598.         od;
  599.     fi;
  600.     return perms;
  601. end;
  602.  
  603. PermutationsList := function ( mset )
  604.     local   m, i;
  605.     mset := ShallowCopy(mset);  Sort( mset );
  606.     m := List( mset, i->true );
  607.     return PermutationsListK(mset,m,Length(mset),Length(mset),[],1);
  608. end;
  609.  
  610.  
  611. #############################################################################
  612. ##
  613. #F  NrPermutationsList( <mset> )  . . .  number of permutations of a multiset
  614. ##
  615. ##  'NrPermutationsList' uses the well known multinomial coefficient formula.
  616. ##
  617. NrPermutationsList := function ( mset )
  618.     local   nr, m, i;
  619.     nr := Factorial( Length(mset) );
  620.     for m  in Set(mset)  do
  621.         nr := nr / Factorial( Number( mset, i->i = m ) );
  622.     od;
  623.     return nr;
  624. end;
  625.  
  626.  
  627. #############################################################################
  628. ##
  629. #F  Derangements( <list> ) . . . . set of fixpointfree permutations of a list
  630. ##
  631. ##  'DerangementsK( <mset>, <m>, <n>, <list>, <k>, <perm>, <i> )' returns the
  632. ##  set of all permutations of the multiset <mset>, which  has size <n>, that
  633. ##  have no element  at the  same position  as <list>,  and that begin   with
  634. ##  '<perm>[[1..<i>-1]]'.   To do this it finds  all elements  of <mset> that
  635. ##  can go at '<perm>[<i>]' and calls itself  recursively for each candidate.
  636. ##  <m> is a boolean list of size <n> that contains  'true' for every element
  637. ##  that we have not  yet taken, so the  candidates for '<perm>[<i>]' are the
  638. ##  elements '<mset>[<l>]' such that '<m>[<l>]' is 'true'.  Some care must be
  639. ##  taken  to  take a candidate   only once if  it  append more  than once in
  640. ##  <mset>.
  641. ##
  642. DerangementsK := function ( mset, m, n, list, k, perm, i )
  643.     local   perms, l;
  644.     if k = 0  then
  645.         perm := ShallowCopy(perm);
  646.         perms := [ perm ];
  647.     else
  648.         perms := [];
  649.         for l  in [1..n]  do
  650.             if m[l] and (l=1 or m[l-1]=false or mset[l]<>mset[l-1])
  651.               and mset[l] <> list[i]  then
  652.                 perm[i] := mset[l];
  653.                 m[l] := false;
  654.                 Append( perms, DerangementsK(mset,m,n,list,k-1,perm,i+1) );
  655.                 m[l] := true;
  656.             fi;
  657.         od;
  658.     fi;
  659.     return perms;
  660. end;
  661.  
  662. Derangements := function ( list )
  663.     local   mset, m, i;
  664.     mset := ShallowCopy(list);  Sort( mset );
  665.     m := List( mset, i->true );
  666.     return DerangementsK(mset,m,Length(mset),list,Length(mset),[],1);
  667. end;
  668.  
  669.  
  670. #############################################################################
  671. ##
  672. #F  NrDerangements( <list> ) .  number of fixpointfree permutations of a list
  673. ##
  674. ##  'NrDerangements' uses well known  identities if <mset>  is a proper  set.
  675. ##  If <mset> is a multiset it  uses 'NrDerangementsK', which works just like
  676. ##  'DerangementsK'.
  677. ##
  678. NrDerangementsK := function ( mset, m, n, list, k, i )
  679.     local   perms, l;
  680.     if k = 0  then
  681.         perms := 1;
  682.     else
  683.         perms := 0;
  684.         for l  in [1..n]  do
  685.             if m[l] and (l=1 or m[l-1]=false or mset[l]<>mset[l-1])
  686.               and mset[l] <> list[i]  then
  687.                 m[l] := false;
  688.                 perms := perms + NrDerangementsK(mset,m,n,list,k-1,i+1);
  689.                 m[l] := true;
  690.             fi;
  691.         od;
  692.     fi;
  693.     return perms;
  694. end;
  695.  
  696. NrDerangements := function ( list )
  697.     local   nr, mset, m, i;
  698.     mset := ShallowCopy(list);  Sort( mset );
  699.     if IsSet(mset)  then
  700.         if Size(mset) = 0  then
  701.             nr := 1;
  702.         elif Size(mset) = 1  then
  703.             nr := 0;
  704.         else
  705.             m := - Factorial(Size(mset));
  706.             nr := 0;
  707.             for i  in [2..Size(mset)]  do
  708.                 m := - m / i;
  709.                 nr := nr + m;
  710.             od;
  711.         fi;
  712.     else
  713.         m := List( mset, i->true );
  714.         nr := NrDerangementsK(mset,m,Length(mset),list,Length(mset),1);
  715.     fi;
  716.     return nr;
  717. end;
  718.  
  719.  
  720. #############################################################################
  721. ##
  722. #F  Permanent( <mat> )  . . . . . . . . . . . . . . . . permanent of a matrix
  723. ##
  724. Permanent2 := function ( mat, n, i, sum )
  725.     local   p,  k;
  726.     if i = n+1  then
  727.         p := 1;
  728.         for k  in sum  do p := p * k;  od;
  729.     else
  730.         p := Permanent2( mat, n, i+1, sum )
  731.            - Permanent2( mat, n, i+1, sum+mat[i] );
  732.     fi;
  733.     return p;
  734. end;
  735.  
  736. Permanent := function ( mat )
  737.     return (-1)^Length(mat) * Permanent2( mat, Length(mat), 1, 0*mat[1] );
  738. end;
  739.  
  740.  
  741. #############################################################################
  742. ##
  743. #F  PartitionsSet( <set> )  . . . . . . . . . . .  set of partitions of a set
  744. ##
  745. ##  'PartitionsSetA( <set>,  <n>, <m>, <o>, <part>,  <i>, <j> )'  returns the
  746. ##  set  of all partitions of  the set <set>, which  has size <n>, that begin
  747. ##  with  '<part>[[1..<i>-1]]'  and   where the    <i>-th   set begins   with
  748. ##  '<part>[<i>][[1..<j>]]'.    To do so  it  does two things.   It finds all
  749. ##  elements of  <mset> that can  go at '<part>[<i>][<j>+1]' and calls itself
  750. ##  recursively for  each candidate.  And it  considers the set '<part>[<i>]'
  751. ##  to be complete  and starts a  new  set '<part>[<i>+1]', which must  start
  752. ##  with the smallest element of <mset> not yet taken, because we require the
  753. ##  returned partitions to be  sorted lexicographically.  <mset> is a boolean
  754. ##  list that contains 'true' for every element of <mset> not yet taken.  <o>
  755. ##  is the position  of '<part>[<i>][<j>]' in  <mset>, so the candidates  for
  756. ##  '<part>[<i>][<j>+1]' are  those elements '<mset>[<l>]'  for  which '<o> <
  757. ##  <l>' and '<m>[<l>]' is 'true'.
  758. ##
  759. ##  'PartitionsSetK( <set>, <n>,  <m>, <o>, <k>,  <part>, <i>, <j> )' returns
  760. ##  the set of all partitions of the set <set>, which has size <n>, that have
  761. ##  '<k>+<i>-1' subsets, and  that begin with '<part>[[1..<i>-1]]'  and where
  762. ##  the <i>-th set begins with '<part>[<i>][[1..<j>]]'.  To do so it does two
  763. ##  things.  It   finds    all  elements   of   <mset>    that  can   go   at
  764. ##  '<part>[<i>][<j>+1]'  and calls  itself  recursively for  each candidate.
  765. ##  And,  if <k> is  larger than 1, it considers  the set '<part>[<i>]' to be
  766. ##  complete and starts a new set '<part>[<i>+1]',  which must start with the
  767. ##  smallest element of <mset> not yet taken, because we require the returned
  768. ##  partitions to be sorted lexicographically.  <mset> is a boolean list that
  769. ##  contains 'true' for every element  of <mset> not yet  taken.  <o> is  the
  770. ##  position    of '<part>[<i>][<j>]'  in    <mset>,  so  the  candidates for
  771. ##  '<part>[<i>][<j>+1]' are those elements  '<mset>[<l>]'  for which '<o>  <
  772. ##  <l>' and '<m>[<l>]' is 'true'.
  773. ##
  774. ##  'PartitionsSet' only  calls   'PartitionsSetA' or  'PartitionsSetK'  with
  775. ##  initial arguments.
  776. ##
  777. PartitionsSetA := function ( set, n, m, o, part, i, j )
  778.     local   parts, npart, l;
  779.     l := Position(m,true);
  780.     if l = false  then
  781.         part := List(part,ShallowCopy);
  782.         parts := [ part ];
  783.     else
  784.         npart := ShallowCopy(part);
  785.         m[l] := false;
  786.         npart[i+1] := [ set[l] ];
  787.         parts := PartitionsSetA(set,n,m,l+1,npart,i+1,1);
  788.         m[l] := true;
  789.         part := ShallowCopy(part);
  790.         part[i] := ShallowCopy(part[i]);
  791.         for l  in [o..n]  do
  792.             if m[l]   then
  793.                 m[l] := false;
  794.                 part[i][j+1] := set[l];
  795.                 Append( parts, PartitionsSetA(set,n,m,l+1,part,i,j+1));
  796.                 m[l] := true;
  797.             fi;
  798.         od;
  799.     fi;
  800.     return parts;
  801. end;
  802.  
  803. PartitionsSetK := function ( set, n, m, o, k, part, i, j )
  804.     local   parts, npart, l;
  805.     l := Position(m,true);
  806.     parts := [];
  807.     if k = 1  then
  808.         part := List(part,ShallowCopy);
  809.         for l  in [k..n]  do
  810.             if m[l]  then
  811.                 Add( part[i], set[l] );
  812.             fi;
  813.         od;
  814.         parts := [ part ];
  815.     elif l <> false  then
  816.         npart := ShallowCopy(part);
  817.         m[l] := false;
  818.         npart[i+1] := [ set[l] ];
  819.         parts := PartitionsSetK(set,n,m,l+1,k-1,npart,i+1,1);
  820.         m[l] := true;
  821.         part := ShallowCopy(part);
  822.         part[i] := ShallowCopy(part[i]);
  823.         for l  in [o..n]  do
  824.             if m[l]  then
  825.                 m[l] := false;
  826.                 part[i][j+1] := set[l];
  827.                 Append( parts, PartitionsSetK(set,n,m,l+1,k,part,i,j+1));
  828.                 m[l] := true;
  829.             fi;
  830.         od;
  831.     fi;
  832.     return parts;
  833. end;
  834.  
  835. PartitionsSet := function ( arg )
  836.     local   parts, set, m, i;
  837.     if Length(arg) = 1  then
  838.         set := arg[1];
  839.         if not IsSet(arg[1])  then
  840.             Error("PartitionsSet: <set> must be a set");
  841.         fi;
  842.         if set = []  then
  843.             parts := [ [  ] ];
  844.         else
  845.             m := List( set, i->true );
  846.             m[1] := false;
  847.             parts := PartitionsSetA(set,Length(set),m,2,[[set[1]]],1,1);
  848.         fi;
  849.     elif Length(arg) = 2  then
  850.         set := arg[1];
  851.         if not IsSet(set)  then
  852.             Error("PartitionsSet: <set> must be a set");
  853.         fi;
  854.         if set = []  then
  855.             if arg[2] = 0  then
  856.                 parts := [ [ ] ];
  857.             else
  858.                 parts := [ ];
  859.             fi;
  860.         else
  861.             m := List( set, i->true );
  862.             m[1] := false;
  863.             parts := PartitionsSetK(
  864.                         set, Length(set), m, 2, arg[2], [[set[1]]], 1, 1 );
  865.         fi;
  866.     else
  867.         Error("usage: PartitionsSet( <n> [, <k>] )");
  868.     fi;
  869.     return parts;
  870. end;
  871.  
  872.  
  873. #############################################################################
  874. ##
  875. #F  NrPartitionsSet( <set> )  . . . . . . . . . number of partitions of a set
  876. ##
  877. NrPartitionsSet := function ( arg )
  878.     local   nr, set, k;
  879.     if Length(arg) = 1  then
  880.         set := arg[1];
  881.         if not IsSet(arg[1])  then
  882.             Error("NrPartitionsSet: <set> must be a set");
  883.         fi;
  884.         nr := Bell( Size(set) );
  885.     elif Length(arg) = 2  then
  886.         set := arg[1];
  887.         if not IsSet(set)  then
  888.             Error("NrPartitionsSet: <set> must be a set");
  889.         fi;
  890.         nr := Stirling2( Size(set), arg[2] );
  891.     else
  892.         Error("usage: NrPartitionsSet( <n> [, <k>] )");
  893.     fi;
  894.     return nr;
  895. end;
  896.  
  897.  
  898. #############################################################################
  899. ##
  900. #F  Partitions( <n> ) . . . . . . . . . . . . set of partitions of an integer
  901. ##
  902. ##  'PartitionsA( <n>, <m>, <part>, <i> )' returns the  set of all partitions
  903. ##  of '<n> +  Sum(<part>[[1..<i>-1]])' that begin with '<part>[[1..<i>-1]]'.
  904. ##  To do so  it finds  all values that  can go  at  '<part>[<i>]' and  calls
  905. ##  itself recursively for each   candidate.  <m> is '<part>[<i>-1]', so  the
  906. ##  candidates  for   '<part>[<i>]'  are '[1..Minimum(<m>,<n>)]',   since  we
  907. ##  require that partitions are nonincreasing.
  908. ##
  909. ##  There is one hack  that needs some comments.   Each call to 'PartitionsA'
  910. ##  contributes one  partition  without   going into recursion,    namely the
  911. ##  'Concatenation(  <part>[[1..<i>-1]], [1,1,...,1]  )'.  Of all  partitions
  912. ##  returned by 'PartitionsA'  this  is the smallest,  i.e.,  it will be  the
  913. ##  first  one in the result  set.   Therefor it is  put  into the result set
  914. ##  before anything else is done.  However it  is not immediately padded with
  915. ##  1, this is  the last  thing  'PartitionsA' does befor  returning.  In the
  916. ##  meantime the  list is  used as a   temporary that is passed  to recursive
  917. ##  invocations.  Note that the fact that each call contributes one partition
  918. ##  without going into recursion means that  the number of recursive calls to
  919. ##  'PartitionsA'  (and the number  of  calls to  'ShallowCopy') is equal  to
  920. ##  'NrPartitions(<n>)'.
  921. ##
  922. ##  'PartitionsK( <n>,  <m>,  <k>, <part>,  <i>  )' returns  the  set of  all
  923. ##  partitions    of  '<n>  +   Sum(<part>[[1..<i>-1]])'  that    have length
  924. ##  '<k>+<i>-1' and that begin with '<part>[[1..<i>-1]]'.   To do so it finds
  925. ##  all values that can go at '<part>[<i>]' and  calls itself recursively for
  926. ##  each  candidate.    <m>  is  '<part>[<i>-1]',   so  the   candidates  for
  927. ##  '<part>[<i>]' must be  less than or  equal to <m>,  since we require that
  928. ##  partitions    are  nonincreasing.   Also    '<part>[<i>]' must    be  \<=
  929. ##  '<n>+1-<k>', since  we need at   least  <k>-1  ones  to  fill the   <k>-1
  930. ##  positions of <part> remaining after  filling '<part>[<i>]'.  On the other
  931. ##  hand '<part>[<i>]'  must be  >=  '<n>/<k>', because otherwise we  can not
  932. ##  fill the <k>-1 remaining positions nonincreasingly.   It is not difficult
  933. ##  to show  that for  each  candidate satisfying these properties   there is
  934. ##  indeed a partition, i.e., we never run into a dead end.
  935. ##
  936. ##  'Partitions'  only  calls  'PartitionsA'  or  'PartitionsK'  with initial
  937. ##  arguments.
  938. ##
  939. PartitionsA := function ( n, m, part, i )
  940.     local   parts, l;
  941.     if n = 0  then
  942.         part := ShallowCopy(part);
  943.         parts := [ part ];
  944.     elif n <= m  then
  945.         part := ShallowCopy(part);
  946.         parts := [ part ];
  947.         for l  in [2..n]  do
  948.             part[i] := l;
  949.             Append( parts, PartitionsA( n-l, l, part, i+1 ) );
  950.         od;
  951.         for l  in [i..i+n-1]  do
  952.             part[l] := 1;
  953.         od;
  954.     else
  955.         part := ShallowCopy(part);
  956.         parts := [ part ];
  957.         for l  in [2..m]  do
  958.             part[i] := l;
  959.             Append( parts, PartitionsA( n-l, l, part, i+1 ) );
  960.         od;
  961.         for l  in [i..i+n-1]  do
  962.             part[l] := 1;
  963.         od;
  964.     fi;
  965.     return parts;
  966. end;
  967.  
  968. PartitionsK := function ( n, m, k, part, i )
  969.     local   parts, l;
  970.     if k = 1  then
  971.         part := ShallowCopy(part);
  972.         part[i] := n;
  973.         parts := [ part ];
  974.     elif n+1-k < m  then
  975.         parts := [];
  976.         for l  in [QuoInt(n+k-1,k)..n+1-k]  do
  977.             part[i] := l;
  978.             Append( parts, PartitionsK( n-l, l, k-1, part, i+1 ) );
  979.         od;
  980.     else
  981.         parts := [];
  982.         for l  in [QuoInt(n+k-1,k)..m]  do
  983.             part[i] := l;
  984.             Append( parts, PartitionsK( n-l, l, k-1, part, i+1 ) );
  985.         od;
  986.     fi;
  987.     return parts;
  988. end;
  989.  
  990. Partitions := function ( arg )
  991.     local   parts;
  992.     if Length(arg) = 1  then
  993.         parts := PartitionsA( arg[1], arg[1], [], 1 );
  994.     elif Length(arg) = 2  then
  995.         if arg[1] = 0  then
  996.             if arg[2] = 0  then
  997.                 parts := [ [  ] ];
  998.             else
  999.                 parts := [  ];
  1000.             fi;
  1001.         else
  1002.             if arg[2] = 0  then
  1003.                 parts := [  ];
  1004.             else
  1005.                 parts := PartitionsK( arg[1], arg[1], arg[2], [], 1 );
  1006.             fi;
  1007.         fi;
  1008.     else
  1009.         Error("usage: Partitions( <n> [, <k>] )");
  1010.     fi;
  1011.     return parts;
  1012. end;
  1013.  
  1014.  
  1015. #############################################################################
  1016. ##
  1017. #F  NrPartitions( <n> ) . . . . . . . . .  number of partitions of an integer
  1018. ##
  1019. ##  To compute $p(n) = NrPartitions(n)$ we use Euler\'s theorem, that asserts
  1020. ##  $p(n) = \sum_{k>0}{ (-1)^{k+1} (p(n-(3m^2-m)/2) + p(n-(3m^2+m)/2)) }$.
  1021. ##
  1022. ##  To compute $p(n,k)$ we use $p(m,1) = p(m,m) = 1$, $p(m,l) = 0$ if $m\<l$,
  1023. ##  and the recurrence  $p(m,l) = p(m-1,l-1) + p(m-l,l)$  if $1 \<   l \< m$.
  1024. ##  This recurrence can be proved by spliting the number of ways to write $m$
  1025. ##  as a  sum of $l$  summands in two subsets,  those  sums that have  1 as a
  1026. ##  summand and those that do not.  The number of ways  to write $m$ as a sum
  1027. ##  of $l$ summands that have 1 as a  summand is $p(m-1,l-1)$, because we can
  1028. ##  take away  the  1 and obtain  a new sums with   $l-1$ summands  and value
  1029. ##  $m-1$.  The number of ways to  write  $m$  as a sum  of $l$ summands such
  1030. ##  that no summand is 1 is $P(m-l,l)$, because we  can  subtract 1 from each
  1031. ##  summand and obtain new sums that still have $l$ summands but value $m-l$.
  1032. ##
  1033. NrPartitions := function ( arg )
  1034.     local   s, n, m, p, k, l;
  1035.  
  1036.     if Length(arg) = 1  then
  1037.         n := arg[1];
  1038.         s := 1;                             # p(0) = 1
  1039.         p := [ s ];
  1040.         for m  in [1..n]  do
  1041.             s := 0;
  1042.             k := 1;
  1043.             l := 1;                         # k*(3*k-1)/2
  1044.             while 0 <= m-(l+k)  do
  1045.                 s := s - (-1)^k * (p[m-l+1] + p[m-(l+k)+1]);
  1046.                 k := k + 1;
  1047.                 l := l + 3*k - 2;
  1048.             od;
  1049.             if 0 <= m-l  then
  1050.                 s := s - (-1)^k * p[m-l+1];
  1051.             fi;
  1052.             p[m+1] := s;
  1053.         od;
  1054.  
  1055.     elif Length(arg) = 2  then
  1056.         if arg[1] = arg[2]  then
  1057.             s := 1;
  1058.         elif arg[1] < arg[2]  or arg[2] = 0  then
  1059.             s := 0;
  1060.         else
  1061.             n := arg[1];  k := arg[2];
  1062.             p := [];
  1063.             for m  in [1..n]  do
  1064.                 p[m] := 1;                  # p(m,1) = 1
  1065.             od;
  1066.             for l  in [2..k]  do
  1067.                 for m  in [l+1..n-l+1]  do
  1068.                     p[m] := p[m] + p[m-l];  # p(m,l) = p(m,l-1) + p(m-l,l)
  1069.                 od;
  1070.             od;
  1071.             s := p[n-k+1];
  1072.         fi;
  1073.  
  1074.     else
  1075.         Error("usage: NrPartitions( <n> [, <k>] )");
  1076.     fi;
  1077.  
  1078.     return s;
  1079. end;
  1080.  
  1081.  
  1082. #############################################################################
  1083. ##
  1084. #F  OrderedPartitions( <n> ) . . . .  set of ordered partitions of an integer
  1085. ##
  1086. ##  'OrderedPartitionsA( <n>,  <part>, <i> )' returns  the set of all ordered
  1087. ##  partitions  of  '<n>  +    Sum(<part>[[1..<i>-1]])'   that  begin    with
  1088. ##  '<part>[[1..<i>-1]]'.   To do    so   it puts   all  possible  values  at
  1089. ##  '<part>[<i>]', which are of course exactely the elements of '[1..n]', and
  1090. ##  calls itself recursively.
  1091. ##
  1092. ##  'OrderedPartitionsK(  <n>,  <k>, <part>,  <i>  )' returns the set  of all
  1093. ##  ordered partitions of   '<n> + Sum(<part>[[1..<i>-1]])'  that have length
  1094. ##  '<k>+<i>-1', and that begin with '<part>[[1..<i>-1]]'.  To  do so it puts
  1095. ##  all possible  values at '<part>[<i>]', which are   of course exactely the
  1096. ##  elements of '[1..<n>-<k>+1]', and calls itself recursively.
  1097. ##
  1098. ##  'OrderedPartitions'      only       calls     'OrderedPartitionsA'     or
  1099. ##  'OrderedPartitionsK' with initial arguments.
  1100. ##
  1101. OrderedPartitionsA := function ( n, part, i )
  1102.     local   parts, l;
  1103.     if n = 0  then
  1104.         part := ShallowCopy(part);
  1105.         parts := [ part ];
  1106.     else
  1107.         part := ShallowCopy(part);
  1108.         parts := [];
  1109.         for l  in [1..n-1]  do
  1110.             part[i] := l;
  1111.             Append( parts, OrderedPartitionsA( n-l, part, i+1 ) );
  1112.         od;
  1113.         part[i] := n;
  1114.         Add( parts, part );
  1115.     fi;
  1116.     return parts;
  1117. end;
  1118.  
  1119. OrderedPartitionsK := function ( n, k, part, i )
  1120.     local   parts, l;
  1121.     if k = 1  then
  1122.         part := ShallowCopy(part);
  1123.         part[i] := n;
  1124.         parts := [ part ];
  1125.     else
  1126.         parts := [];
  1127.         for l  in [1..n-k+1]  do
  1128.             part[i] := l;
  1129.             Append( parts, OrderedPartitionsK( n-l, k-1, part, i+1 ) );
  1130.         od;
  1131.     fi;
  1132.     return parts;
  1133. end;
  1134.  
  1135. OrderedPartitions := function ( arg )
  1136.     local   parts;
  1137.     if Length(arg) = 1  then
  1138.         parts := OrderedPartitionsA( arg[1], [], 1 );
  1139.     elif Length(arg) = 2  then
  1140.         if arg[1] = 0  then
  1141.             if arg[2] = 0  then
  1142.                 parts := [ [  ] ];
  1143.             else
  1144.                 parts := [  ];
  1145.             fi;
  1146.         else
  1147.             if arg[2] = 0  then
  1148.                 parts := [  ];
  1149.             else
  1150.                 parts := OrderedPartitionsK( arg[1], arg[2], [], 1 );
  1151.             fi;
  1152.         fi;
  1153.     else
  1154.         Error("usage: OrderedPartitions( <n> [, <k>] )");
  1155.     fi;
  1156.     return parts;
  1157. end;
  1158.  
  1159.  
  1160. #############################################################################
  1161. ##
  1162. #F  NrOrderedPartitions( <n> ) . . number of ordered partitions of an integer
  1163. ##
  1164. ##  'NrOrderedPartitions' uses well known identities to compute the number of
  1165. ##  ordered partitions of <n>.
  1166. ##
  1167. NrOrderedPartitions := function ( arg )
  1168.     local   nr;
  1169.     if Length(arg) = 1  then
  1170.         if arg[1] = 0  then
  1171.             nr := 1;
  1172.         else
  1173.             nr := 2^(arg[1]-1);
  1174.         fi;
  1175.     elif Length(arg) = 2  then
  1176.         if arg[1] = 0  then
  1177.             if arg[2] = 0  then
  1178.                 nr := 1;
  1179.             else
  1180.                 nr := 0;
  1181.             fi;
  1182.         else
  1183.             nr := Binomial(arg[1]-1,arg[2]-1);
  1184.         fi;
  1185.     else
  1186.         Error("usage: NrOrderedPartitions( <n> [, <k>] )");
  1187.     fi;
  1188.     return nr;
  1189. end;
  1190.  
  1191.  
  1192. #############################################################################
  1193. ##
  1194. #F  RestrictedPartitions( <n>, <set> )  . restricted partitions of an integer
  1195. ##
  1196. ##  'RestrictedPartitionsA( <n>, <set>, <m>,  <part>, <i> )' returns the  set
  1197. ##  of  all partitions of  '<n> +  Sum(<part>[[1..<i>-1]])' that contain only
  1198. ##  elements of <set> and that begin with '<part>[[1..<i>-1]]'.   To do so it
  1199. ##  finds all elements of <set> that can go at '<part>[<i>]' and calls itself
  1200. ##  recursively for each candidate.  <m>  is the position of  '<part>[<i>-1]'
  1201. ##  in   <set>,  so the   candidates for  '<part>[<i>]'  are  the elements of
  1202. ##  '<set>[[1..<m>]]' that  are less    than  <n>,  since we   require   that
  1203. ##  partitions are nonincreasing.
  1204. ##
  1205. ##  'RestrictedPartitionsK( <n>, <set>, <m>, <k>,  <part>, <i> )' returns the
  1206. ##  set  of all  partitions of  '<n>  + Sum(<part>[[1..<i>-1]])' that contain
  1207. ##  only elements of <set>, that have length '<k>+<i>-1', and that begin with
  1208. ##  '<part>[[1..<i>-1]]'.   To do so it finds  all elements fo <set> that can
  1209. ##  go at '<part>[<i>]' and calls itself recursively for each candidate.  <m>
  1210. ##  is  the  position  of '<part>[<i>-1]'   in  <set>, so the  candidates for
  1211. ##  '<part>[<i>]' are the elements  of '<set>[[1..<m>]]'  that are less  than
  1212. ##  <n>, since we require that partitions are nonincreasing.
  1213. ##
  1214. RestrictedPartitionsA := function ( n, set, m, part, i )
  1215.     local   parts, l;
  1216.     if n = 0  then
  1217.         part := ShallowCopy(part);
  1218.         parts := [ part ];
  1219.     else
  1220.         part := ShallowCopy(part);
  1221.         if n mod set[1] = 0  then
  1222.             parts := [ part ];
  1223.         else
  1224.             parts := [ ];
  1225.         fi;
  1226.         for l  in [2..m]  do
  1227.             if set[l] <= n  then
  1228.                 part[i] := set[l];
  1229.                 Append(parts,RestrictedPartitionsA(n-set[l],set,l,part,i+1));
  1230.             fi;
  1231.         od;
  1232.         if n mod set[1] = 0  then
  1233.             for l  in [i..i+n/set[1]-1]  do
  1234.                 part[l] := set[1];
  1235.             od;
  1236.         fi;
  1237.     fi;
  1238.     return parts;
  1239. end;
  1240.  
  1241. RestrictedPartitionsK := function ( n, set, m, k, part, i )
  1242.     local   parts, l;
  1243.     if k = 1  then
  1244.         if n in set  then
  1245.             part := ShallowCopy(part);
  1246.             part[i] := n;
  1247.             parts := [ part ];
  1248.         else
  1249.             parts := [];
  1250.         fi;
  1251.     else
  1252.         part := ShallowCopy(part);
  1253.         parts := [ ];
  1254.         for l  in [1..m]  do
  1255.             if set[l]+(k-1)*set[1] <= n  and n <= k*set[l]  then
  1256.                 part[i] := set[l];
  1257.                 Append(parts,
  1258.                        RestrictedPartitionsK(n-set[l],set,l,k-1,part,i+1) );
  1259.             fi;
  1260.         od;
  1261.     fi;
  1262.     return parts;
  1263. end;
  1264.  
  1265. RestrictedPartitions := function ( arg )
  1266.     local   parts;
  1267.     if Length(arg) = 2  then
  1268.         parts := RestrictedPartitionsA(arg[1],arg[2],Length(arg[2]),[],1);
  1269.     elif Length(arg) = 3  then
  1270.         if arg[1] = 0  then
  1271.             if arg[3] = 0  then
  1272.                 parts := [ [  ] ];
  1273.             else
  1274.                 parts := [  ];
  1275.             fi;
  1276.         else
  1277.             if arg[2] = 0  then
  1278.                 parts := [  ];
  1279.             else
  1280.                 parts := RestrictedPartitionsK(
  1281.                              arg[1], arg[2], Length(arg[2]), arg[3], [], 1 );
  1282.             fi;
  1283.         fi;
  1284.     else
  1285.         Error("usage: RestrictedPartitions( <n>, <set> [, <k>] )");
  1286.     fi;
  1287.     return parts;
  1288. end;
  1289.  
  1290.  
  1291. #############################################################################
  1292. ##
  1293. #F  NrRestrictedPartitions(<n>,<set>) . . . . number of restricted partitions
  1294. ##
  1295. #N  22-Jul-91 martin there should be a better way to do this for given <k>
  1296. ##
  1297. NrRestrictedPartitionsK := function ( n, set, m, k, part, i )
  1298.     local   parts, l;
  1299.     if k = 1  then
  1300.         if n in set  then
  1301.             parts := 1;
  1302.         else
  1303.             parts := 0;
  1304.         fi;
  1305.     else
  1306.         part := ShallowCopy(part);
  1307.         parts := 0;
  1308.         for l  in [1..m]  do
  1309.             if set[l]+(k-1)*set[1] <= n  and n <= k*set[l]  then
  1310.                 part[i] := set[l];
  1311.                 parts := parts +
  1312.                         NrRestrictedPartitionsK(n-set[l],set,l,k-1,part,i+1);
  1313.             fi;
  1314.         od;
  1315.     fi;
  1316.     return parts;
  1317. end;
  1318.  
  1319. NrRestrictedPartitions := function ( arg )
  1320.     local  s, n, set, m, p, k, l;
  1321.  
  1322.     if Length(arg) = 2  then
  1323.         n := arg[1];
  1324.         set := arg[2];
  1325.         p := [];
  1326.         for m  in [1..n+1]  do
  1327.             if (m-1) mod set[1] = 0  then
  1328.                 p[m] := 1;
  1329.             else
  1330.                 p[m] := 0;
  1331.             fi;
  1332.         od;
  1333.         for l  in Sublist( set, [2..Length(set)] )  do
  1334.             for m  in [l+1..n+1]  do
  1335.                 p[m] := p[m] + p[m-l];
  1336.             od;
  1337.         od;
  1338.         s := p[n+1];
  1339.  
  1340.     elif Length(arg) = 3  then
  1341.         if arg[1] = 0  and arg[3] = 0  then
  1342.             s := 1;
  1343.         elif arg[1] < arg[3]  or arg[3] = 0  then
  1344.             s := 0;
  1345.         else
  1346.             s := NrRestrictedPartitionsK(
  1347.                         arg[1], arg[2], Length(arg[2]), arg[3], [], 1 );
  1348.         fi;
  1349.  
  1350.     else
  1351.         Error("usage: NrRestrictedPartitions( <n>, <set> [, <k>] )");
  1352.     fi;
  1353.  
  1354.     return s;
  1355. end;
  1356.  
  1357.  
  1358. #############################################################################
  1359. ##
  1360. #F  SignPartition( <pi> ) . . . . . . . . . . . . .  signum of partition <pi>
  1361. ##
  1362. SignPartition := function(pi)
  1363.    
  1364.    return (-1)^(Sum(pi) - Length(pi));
  1365.  
  1366. end;
  1367.  
  1368.  
  1369. #############################################################################
  1370. ##
  1371. #F  AssociatedPartition( <pi> ) . . . . . .  the associated partition of <pi>
  1372. ##
  1373. ##  'AssociatedPartition'  returns the associated partition  of the partition
  1374. ##  <pi> which is obtained by transposing the corresponding Young diagram.
  1375. ##
  1376. AssociatedPartition := function(pi)
  1377.  
  1378.    local i, j, mu;
  1379.  
  1380.    mu:= List([1..pi[1]], x->0);
  1381.    for i in pi do
  1382.       for j in [1..i] do
  1383.      mu[j]:= mu[j]+1;
  1384.       od;
  1385.    od;
  1386.  
  1387.    return(mu);
  1388.  
  1389. end;
  1390.  
  1391.  
  1392. #############################################################################
  1393. ##
  1394. #F  PowerPartition( <pi>, <k> ) . . . . . . . . . . . .  power of a partition
  1395. ##
  1396. ##  'PowerPartition'  returns the partition corresponding to the <k>-th power
  1397. ##  of a permutation with cycle structure <pi>.
  1398. ##
  1399. PowerPartition := function(pi, k)
  1400.  
  1401.    local res, i, d, part;
  1402.  
  1403.    res:= [];
  1404.    for part in pi do
  1405.       d:= GcdInt(part, k);
  1406.       for i in [1..d] do
  1407.          Add(res, part/d);
  1408.       od;
  1409.    od;
  1410.    Sort(res);
  1411.  
  1412.    return Reversed(res);
  1413.  
  1414. end;
  1415.  
  1416.  
  1417. #############################################################################
  1418. ##
  1419. #F  PartitionTuples( <n>, <r> ) . . . . . . . . . <r> partitions with sum <n>
  1420. ##
  1421. ##  'PartitionTuples'  returns the list of all <r>-tuples of partitions which
  1422. ##  together form a partition of <n>.
  1423. ##
  1424. PartitionTuples := function(n, r)
  1425.  
  1426.    local m, k, pm, i, t, t1, s, res, empty;
  1427.  
  1428.    empty:= rec(tup:= List([1..r], x->[]), pos:= List([1..n-1], x-> 1));
  1429.  
  1430.    if n = 0 then
  1431.       return [empty.tup];
  1432.    fi;
  1433.  
  1434.    pm:= List([1..n-1], x->[]);
  1435.    for m in [1..QuoInt(n,2)] do 
  1436.  
  1437.       # the m-cycle in all possible places.
  1438.       for i in [1..r] do
  1439.          s:= Copy(empty);
  1440.          s.tup[i]:= [m];
  1441.          s.pos[m]:= i;
  1442.          Add(pm[m], s);
  1443.       od;
  1444.  
  1445.       # add the m-cycle to everything you know.
  1446.       for k in [m+1..n-m] do
  1447.          for t in pm[k-m] do
  1448.             for i in [t.pos[m]..r] do
  1449.                t1:= Copy(t);
  1450.                s:= [m];
  1451.                Append(s, t.tup[i]);
  1452.                t1.tup[i]:= s;
  1453.                t1.pos[m]:= i;
  1454.                Add(pm[k], t1);
  1455.             od;
  1456.          od;
  1457.       od;
  1458.    od;
  1459.  
  1460.    #  collect.
  1461.    res:= [];
  1462.    for k in [1..n-1] do
  1463.       for t in pm[n-k] do
  1464.          for i in [t.pos[k]..r] do 
  1465.             t1:= Copy(t.tup);
  1466.             s:= [k];
  1467.             Append(s, t.tup[i]);
  1468.             t1[i]:= s;
  1469.             Add(res, t1);
  1470.          od;
  1471.       od;
  1472.    od;
  1473.    for i in [1..r] do
  1474.       s:= Copy(empty.tup);
  1475.       s[i]:= [n];
  1476.       Add(res, s);
  1477.    od;
  1478.  
  1479.    return res;
  1480.  
  1481. end;
  1482.  
  1483.  
  1484. #############################################################################
  1485. ##
  1486. #F  Lucas(<P>,<Q>,<k>)  . . . . . . . . . . . . . . value of a lucas sequence
  1487. ##
  1488. ##  'Lucas' uses the following relations to compute the result in $O(log(k))$
  1489. ##  $U_{2k} = U_k V_k,        U_{2k+1} = (P U_{2k} + V_{2k}) / 2$ and
  1490. ##  $V_{2k} = V_k^2 - 2 Q^k,  V_{2k+1} = ((P^2-4Q) U_{2k} + P V_{2k}) / 2$.
  1491. ##
  1492. Lucas := function ( P, Q, k )
  1493.     local   l;
  1494.     if k = 0  then
  1495.         l := [ 0, 2, 1 ];
  1496.     elif k < 0  then
  1497.         l := Lucas( P, Q, -k );
  1498.         l := [ -l[1]/l[3], l[2]/l[3], 1/l[3] ];
  1499.     elif k mod 2 = 0  then
  1500.         l := Lucas( P, Q, k/2 );
  1501.         l := [ l[1]*l[2], l[2]^2-2*l[3], l[3]^2 ];
  1502.     else
  1503.         l := Lucas( P, Q, k-1 );
  1504.         l := [ (P*l[1]+l[2])/2, ((P^2-4*Q)*l[1]+P*l[2])/2, Q*l[3] ];
  1505.     fi;
  1506.     return l;
  1507. end;
  1508.  
  1509.  
  1510. #############################################################################
  1511. ##
  1512. #F  Fibonacci( <n> )  . . . . . . . . . . . . value of the Fibonacci sequence
  1513. ##
  1514. ##  A recursive  Fibonacci needs $O( Fibonacci(n) ) = O(2^n)$ bit operations.
  1515. ##  An iterative version performs $n$ additions, the <i>th involving integers
  1516. ##  with $i$ bits,  so  we  need $\sum_{i=1}^{n}{i} = O(n^2)$ bit operations.
  1517. ##  The binary recursion of 'Lucas' reduces the number of calls to $log2(n)$.
  1518. ##  The number of bit operations now is $O(n)$, i.e., the size of the result.
  1519. ##
  1520. Fibonacci := function ( n )
  1521.     return Lucas( 1, -1, n )[ 1 ];
  1522. end;
  1523.  
  1524.  
  1525. #############################################################################
  1526. ##
  1527. #F  Bernoulli( <n> )  . . . . . . . . . . . . value of the Bernoulli sequence
  1528. ##
  1529. Bernoulli2 := [-1/2,1/6,0,-1/30,0,1/42,0,-1/30,0,5/66,0,-691/2730,0,7/6];
  1530.  
  1531. Bernoulli := function ( n )
  1532.     local   brn, bin, i, j;
  1533.     if   n < 0  then
  1534.         Error("Bernoulli: <n> must be nonnegative");
  1535.     elif n = 0  then
  1536.         brn := 1;
  1537.     elif n = 1  then
  1538.         brn := -1/2;
  1539.     elif n mod 2 = 1  then
  1540.         brn := 0;
  1541.     elif n <= Length(Bernoulli2)  then
  1542.         brn := Bernoulli2[n];
  1543.     else
  1544.         for i  in [Length(Bernoulli2)+1..n]  do
  1545.             if i mod 2 = 1  then
  1546.                 Bernoulli2[i] := 0;
  1547.             else
  1548.                 bin := 1;
  1549.                 brn := 1;
  1550.                 for j  in [1..i-1]  do
  1551.                     bin := (i+2-j)/j * bin;
  1552.                     brn := brn + bin * Bernoulli2[j];
  1553.                 od;
  1554.                 Bernoulli2[i] := - brn / (i+1);
  1555.             fi;
  1556.         od;
  1557.         brn := Bernoulli2[n];
  1558.     fi;
  1559.     return brn;
  1560. end;
  1561.  
  1562.  
  1563.  
  1564.