# source:sage/combinat/combinat.py@6995:8c4757db041e

Revision 6995:8c4757db041e, 72.1 KB checked in by William Stein <wstein@…>, 6 years ago (diff)

Working on fixing up the doctesting in the Sage library (removing sage.:; fixing ploting doctesting, so .show() works)

Line
1r"""
2Combinatorial Functions.
3
4AUTHORS:
5        -- David Joyner (2006-07), initial implementation.
6        -- William Stein (2006-07), editing of docs and code; many optimizations,
7                      refinements, and bug fixes in corner cases
8        -- DJ (2006-09): bug fix for combinations, added permutations_iterator,
9                      combinations_iterator from Python Cookbook, edited docs.
10
11This module implements some combinatorial functions, as listed
12below. For a more detailed description, see the relevant docstrings.
13
14Sequences:
15\begin{itemize}
16\item Bell numbers, \code{bell_number}
17
18\item Bernoulli numbers, \code{bernoulli_number} (though PARI's bernoulli is
19  better)
20
21\item Catalan numbers, \code{catalan_number} (not to be confused with the
22  Catalan constant)
23
24\item Eulerian/Euler numbers, \code{euler_number} (Maxima)
25
26\item Fibonacci numbers, \code{fibonacci} (PARI) and \code{fibonacci_number} (GAP)
27  The PARI version is better.
28
29\item Lucas numbers, \code{lucas_number1}, \code{lucas_number2}.
30
31\item Stirling numbers, \code{stirling_number1}, \code{stirling_number2}.
32\end{itemize}
33
34Set-theoretic constructions:
35\begin{itemize}
36\item Combinations of a multiset, \code{combinations}, \code{combinations_iterator},
37and \code{number_of_combinations}. These are unordered selections without
38repetition of k objects from a multiset S.
39
40\item Arrangements of a multiset, \code{arrangements} and \code{number_of_arrangements}
41  These are ordered selections without repetition of k objects from a
42  multiset S.
43
44\item Derangements of a multiset, \code{derangements} and \code{number_of_derangements}.
45
46\item Tuples of a multiset, \code{tuples} and \code{number_of_tuples}.
47  An ordered tuple of length k of set S is a ordered selection with
48  repetitions of S and is represented by a sorted list of length k
49  containing elements from S.
50
51\item Unordered tuples of a set, \code{unordered_tuple} and \code{number_of_unordered_tuples}.
52  An unordered tuple of length k of set S is a unordered selection with
53  repetitions of S and is represented by a sorted list of length k
54  containing elements from S.
55
56\item Permutations of a multiset, \code{permutations}, \code{permutations_iterator},
57\code{number_of_permutations}. A permutation is a list that contains exactly the same elements but
58possibly in different order.
59\end{itemize}
60
61Partitions:
62\begin{itemize}
63\item Partitions of a set, \code{partitions_set}, \code{number_of_partitions_set}.
64  An unordered partition of set S is a set of pairwise disjoint
65  nonempty sets with union S and is represented by a sorted list of
66  such sets.
67
68\item Partitions of an integer, \code{partitions_list}, \code{number_of_partitions_list}.
69  An unordered partition of n is an unordered sum
70  $n = p_1+p_2 +\ldots+ p_k$ of positive integers and is represented by
71  the list $p = [p_1,p_2,\ldots,p_k]$, in nonincreasing order, i.e.,
72  $p1\geq p_2 ...\geq p_k$.
73
74\item Ordered partitions of an integer, \code{ordered_partitions},
75  \code{number_of_ordered_partitions}.
76  An ordered partition of n is an ordered sum $n = p_1+p_2 +\ldots+ p_k$
77  of positive integers and is represented by
78  the list $p = [p_1,p_2,\ldots,p_k]$, in nonincreasing order, i.e.,
79  $p1\geq p_2 ...\geq p_k$.
80
81\item Restricted partitions of an integer, \code{partitions_restricted},
82  \code{number_of_partitions_restricted}.
83  An unordered restricted partition of n is an unordered sum
84  $n = p_1+p_2 +\ldots+ p_k$ of positive integers $p_i$ belonging to a
85  given set $S$, and is represented by the list $p = [p_1,p_2,\ldots,p_k]$,
86  in nonincreasing order, i.e., $p1\geq p_2 ...\geq p_k$.
87
88\item \code{partitions_greatest}
89   implements a special type of restricted partition.
90
91\item \code{partitions_greatest_eq} is another type of restricted partition.
92
93\item Tuples of partitions, \code{partition_tuples},
94  \code{number_of_partition_tuples}.
95  A $k$-tuple of partitions is represented by a list of all $k$-tuples
96  of partitions which together form a partition of $n$.
97
98\item Powers of a partition, \code{partition_power(pi, k)}.
99  The power of a partition corresponds to the $k$-th power of a
100  permutation with cycle structure $\pi$.
101
102\item Sign of a partition, \code{partition_sign( pi ) }
103  This means the sign of a permutation with cycle structure given by the
104  partition pi.
105
106\item Associated partition, \code{partition_associated( pi )}
107  The associated'' (also called conjugate'' in the literature)
108  partition of the partition pi which is obtained by transposing the
109  corresponding Ferrers diagram.
110
111\item Ferrers diagram, \code{ferrers_diagram}.
112  Analogous to the Young diagram of an irredicible representation
113  of $S_n$.
114  \end{itemize}
115
116Related functions:
117
118\begin{itemize}
119\item Bernoulli polynomials, \code{bernoulli_polynomial}
120\end{itemize}
121
122Implemented in other modules (listed for completeness):
123
124The \module{sage.rings.arith} module contains
125the following combinatorial functions:
126\begin{itemize}
127 \item binomial
128     the binomial coefficient (wrapped from PARI)
129 \item factorial (wrapped from PARI)
130 \item partition (from the Python Cookbook)
131     Generator of the list of all the partitions of the integer $n$.
132 \item \code{number_of_partitions} (wrapped from PARI)
133     the *number* of partitions:
134 \item \code{falling_factorial}
135     Definition: for integer $a \ge 0$ we have $x(x-1) \cdots (x-a+1)$.
136     In all other cases we use the GAMMA-function:
137     $\frac {\Gamma(x+1)} {\Gamma(x-a+1)}$.
138 \item \code{rising_factorial}
139     Definition: for integer $a \ge 0$ we have $x(x+1) \cdots (x+a-1)$.
140     In all other cases we use the GAMMA-function:
141     $\frac {\Gamma(x+a)} {\Gamma(x)}$.
142 \item gaussian_binomial
143     the gaussian binomial
144     $$145 \binom{n}{k}_q = \frac{(1-q^m)(1-q^{m-1})\cdots (1-q^{m-r+1})} 146 {(1-q)(1-q^2)\cdots (1-q^r)}. 147$$
148\end{itemize}
149
150The \module{sage.groups.perm_gps.permgroup_elements} contains the following
151combinatorial functions:
152\begin{itemize}
153\item matrix method of PermutationGroupElement yielding the permutation
154     matrix of the group element.
155\end{itemize}
156
157\begin{verbatim}
158TODO:
159   GUAVA commands:
160    * MOLS returns a list of n Mutually Orthogonal Latin Squares (MOLS).
161    * HadamardMat returns a Hadamard matrix of order n.
162    * VandermondeMat
163    * GrayMat returns a list of all different vectors of length n over
164      the field F, using Gray ordering.
165   Not in GAP:
166    * Rencontres numbers
167      http://en.wikipedia.org/wiki/Rencontres_number
168\end{verbatim}
169
170REFERENCES:
171      \url{http://en.wikipedia.org/wiki/Twelvefold_way} (general reference)
172
173"""
174
175#*****************************************************************************
176#       Copyright (C) 2006 David Joyner <wdjoyner@gmail.com>,
177#                     2007 Mike Hansen <mhansen@gmail.com>,
178#                     2006 William Stein <wstein@gmail.com>
179#
180#  Distributed under the terms of the GNU General Public License (GPL)
181#
182#    This code is distributed in the hope that it will be useful,
183#    but WITHOUT ANY WARRANTY; without even the implied warranty of
184#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
185#    General Public License for more details.
186#
187#  The full text of the GPL is available at:
188#
189#                  http://www.gnu.org/licenses/
190#*****************************************************************************
191
192import os
193
194from sage.interfaces.all import gap, maxima
195from sage.rings.all import QQ, RR, ZZ
196from sage.rings.arith import binomial
197from sage.misc.sage_eval import sage_eval
198from sage.libs.all import pari
199from sage.rings.arith import factorial
200from random import randint
201from sage.misc.misc import prod
202from sage.structure.sage_object import SageObject
203import __builtin__
204from sage.algebras.algebra import Algebra
205from sage.algebras.algebra_element import AlgebraElement
206import sage.structure.parent_base
207import partitions as partitions_ext
208
209######### combinatorial sequences
210
211def bell_number(n):
212    r"""
213    Returns the n-th Bell number (the number of ways to partition a
214    set of n elements into pairwise disjoint nonempty subsets).
215
216    If $n \leq 0$, returns $1$.
217
218    Wraps GAP's Bell.
219
220    EXAMPLES:
221        sage: bell_number(10)
222        115975
223        sage: bell_number(2)
224        2
225        sage: bell_number(-10)
226        1
227        sage: bell_number(1)
228        1
229        sage: bell_number(1/3)
230        Traceback (most recent call last):
231        ...
232        TypeError: no coercion of this rational to integer
233    """
234    ans=gap.eval("Bell(%s)"%ZZ(n))
235    return ZZ(eval(ans))
236
237## def bernoulli_number(n):
238##     r"""
239##     Returns the n-th Bernoulli number $B_n$; $B_n/n!$ is the
240##     coefficient of $x^n$ in the power series of $x/(e^x-1)$.
241##     Wraps GAP's Bernoulli.
242##     EXAMPLES:
243##         sage: bernoulli_number(50)
244##         495057205241079648212477525/66
245##     """
246##     ans=gap.eval("Bernoulli(%s)"%n)
247##     return QQ(ans)   ## use QQ, not eval, here
248
249def catalan_number(n):
250    r"""
251    Returns the n-th Catalan number
252
253      Catalan numbers: The $n$-th Catalan number is given directly in terms of
254      binomial coefficients by
255      $256 C_n = \frac{1}{n+1}{2n\choose n} = \frac{(2n)!}{(n+1)!\,n!} 257 \qquad\mbox{ for }n\ge 0. 258$
259
260      Consider the set $S = \{ 1, ..., n \}$. A {\it noncrossing partition} of $S$
261      is a partition in which no two blocks "cross" each other, i.e., if a
262      and b belong to one block and x and y to another, they are not arranged
263      in the order axby. $C_n$ is the number of noncrossing partitions of the set $S$.
264      There are many other interpretations (see REFERENCES).
265
266    When $n=-1$, this function raises a ZeroDivisionError; for other $n<0$ it
267    returns $0$.
268
269    EXAMPLES:
270        sage: [catalan_number(i) for i in range(7)]
271        [1, 1, 2, 5, 14, 42, 132]
272        sage: maxima.eval("-(1/2)*taylor (sqrt (1-4*x^2), x, 0, 15)")
273        '-1/2+x^2+x^4+2*x^6+5*x^8+14*x^10+42*x^12+132*x^14'
274        sage: [catalan_number(i) for i in range(-7,7) if i != -1]
275        [0, 0, 0, 0, 0, 0, 1, 1, 2, 5, 14, 42, 132]
276        sage: catalan_number(-1)
277        Traceback (most recent call last):
278        ...
279        ZeroDivisionError: Rational division by zero
280
281    REFERENCES:
282    \begin{itemize}
283      \item \url{http://en.wikipedia.org/wiki/Catalan_number}
284      \item \url{http://www-history.mcs.st-andrews.ac.uk/~history/Miscellaneous/CatalanNumbers/catalan.html}
285    \end{itemize}
286
287    """
288    n = ZZ(n)
289    return binomial(2*n,n)/(n+1)
290
291def euler_number(n):
292    """
293    Returns the n-th Euler number
294
295    IMPLEMENTATION: Wraps Maxima's euler.
296
297    EXAMPLES:
298        sage: [euler_number(i) for i in range(10)]
299        [1, 0, -1, 0, 5, 0, -61, 0, 1385, 0]
300        sage: maxima.eval("taylor (2/(exp(x)+exp(-x)), x, 0, 10)")
301        '1-x^2/2+5*x^4/24-61*x^6/720+277*x^8/8064-50521*x^10/3628800'
302        sage: [euler_number(i)/factorial(i) for i in range(11)]
303        [1, 0, -1/2, 0, 5/24, 0, -61/720, 0, 277/8064, 0, -50521/3628800]
304        sage: euler_number(-1)
305        Traceback (most recent call last):
306        ...
307        ValueError: n (=-1) must be a nonnegative integer
308
309    REFERENCES:
310        http://en.wikipedia.org/wiki/Euler_number
311    """
312    n = ZZ(n)
313    if n < 0:
314        raise ValueError, "n (=%s) must be a nonnegative integer"%n
315    return ZZ(maxima.eval("euler(%s)"%n))
316
317def fibonacci(n, algorithm="pari"):
318    """
319    Returns then n-th Fibonacci number. The Fibonacci sequence $F_n$
320    is defined by the initial conditions $F_1=F_2=1$ and the
321    recurrence relation $F_{n+2} = F_{n+1} + F_n$. For negative $n$ we
322    define $F_n = (-1)^{n+1}F_{-n}$, which is consistent with the
323    recurrence relation.
324
325    For negative $n$, define $F_{n} = -F_{|n|}$.
326
327    INPUT:
328        algorithm -- string:
329                     "pari" -- (default) -- use the PARI C library's fibo function.
330                     "gap" -- use GAP's Fibonacci function
331
332    NOTE: PARI is tens to hundreds of times faster than GAP here;
333          moreover, PARI works for every large input whereas GAP
334          doesn't.
335
336    EXAMPLES:
337        sage: fibonacci(10)
338        55
339        sage: fibonacci(10, algorithm='gap')
340        55
341
342        sage: fibonacci(-100)
343        -354224848179261915075
344        sage: fibonacci(100)
345        354224848179261915075
346
347        sage: fibonacci(0)
348        0
349        sage: fibonacci(1/2)
350        Traceback (most recent call last):
351        ...
352        TypeError: no coercion of this rational to integer
353    """
354    n = ZZ(n)
355    if algorithm == 'pari':
356        return ZZ(pari(n).fibonacci())
357    elif algorithm == 'gap':
358        return ZZ(gap.eval("Fibonacci(%s)"%n))
359    else:
360        raise ValueError, "no algorithm %s"%algorithm
361
362def lucas_number1(n,P,Q):
363    """
364    Returns then n-th Lucas number of the first kind'' (this is not
365    standard terminology). The Lucas sequence $L^{(1)}_n$ is defined
366    by the initial conditions $L^{(1)}_1=0$, $L^{(1)}_2=1$ and the recurrence
367    relation $L^{(1)}_{n+2} = P*L^{(1)}_{n+1} - Q*L^{(1)}_n$.
368
369    Wraps GAP's Lucas(...)[1].
370
371    P=1, Q=-1 fives the Fibonacci sequence.
372
373    INPUT:
374        n -- integer
375        P, Q -- integer or rational numbers
376
377    OUTPUT:
378        integer or rational number
379
380    EXAMPLES:
381        sage: lucas_number1(5,1,-1)
382        5
383        sage: lucas_number1(6,1,-1)
384        8
385        sage: lucas_number1(7,1,-1)
386        13
387        sage: lucas_number1(7,1,-2)
388        43
389
390        sage: lucas_number1(5,2,3/5)
391        229/25
392        sage: lucas_number1(5,2,1.5)
393        Traceback (most recent call last):
394        ...
395        TypeError: no implicit coercion of element to the rational numbers
396
397    There was a conjecture that the sequence $L_n$ defined by
398    $L_{n+2} = L_{n+1} + L_n$, $L_1=1$, $L_2=3$, has the property
399    that $n$ prime implies that $L_n$ is prime.
400
401        sage: lucas = lambda n:(5/2)*lucas_number1(n,1,-1)+(1/2)*lucas_number2(n,1,-1)
402        sage: [[lucas(n),is_prime(lucas(n)),n+1,is_prime(n+1)] for n in range(15)]
403        [[1, False, 1, False],
404         [3, True, 2, True],
405         [4, False, 3, True],
406         [7, True, 4, False],
407         [11, True, 5, True],
408         [18, False, 6, False],
409         [29, True, 7, True],
410         [47, True, 8, False],
411         [76, False, 9, False],
412         [123, False, 10, False],
413         [199, True, 11, True],
414         [322, False, 12, False],
415         [521, True, 13, True],
416         [843, False, 14, False],
417         [1364, False, 15, False]]
418
419    Can you use SAGE to find a counterexample to the conjecture?
420    """
421    ans=gap.eval("Lucas(%s,%s,%s)[1]"%(QQ._coerce_(P),QQ._coerce_(Q),ZZ(n)))
422    return sage_eval(ans)
423
424def lucas_number2(n,P,Q):
425    r"""
426    Returns then n-th Lucas number of the second kind'' (this is not
427    standard terminology). The Lucas sequence $L^{(2)}_n$ is defined
428    by the initial conditions $L^{(2)}_1=2$, $L^{(2)}_2=P$ and the recurrence
429    relation $L^{(2)}_{n+2} = P*L^{(2)}_{n+1} - Q*L^{(2)}_n$.
430
431    Wraps GAP's Lucas(...)[2].
432
433    INPUT:
434        n -- integer
435        P, Q -- integer or rational numbers
436
437    OUTPUT:
438        integer or rational number
439
440    EXAMPLES:
441        sage: [lucas_number2(i,1,-1) for i in range(10)]
442        [2, 1, 3, 4, 7, 11, 18, 29, 47, 76]
443        sage: [fibonacci(i-1)+fibonacci(i+1) for i in range(10)]
444        [2, 1, 3, 4, 7, 11, 18, 29, 47, 76]
445
446        sage: n = lucas_number2(5,2,3); n
447        2
448        sage: type(n)
449        <type 'sage.rings.integer.Integer'>
450        sage: n = lucas_number2(5,2,-3/9); n
451        418/9
452        sage: type(n)
453        <type 'sage.rings.rational.Rational'>
454
455    The case P=1, Q=-1 is the Lucas sequence in Brualdi's
456    {\bf Introductory Combinatorics}, 4th ed., Prentice-Hall, 2004:
457        sage: [lucas_number2(n,1,-1) for n in range(10)]
458        [2, 1, 3, 4, 7, 11, 18, 29, 47, 76]
459    """
460    ans=gap.eval("Lucas(%s,%s,%s)[2]"%(QQ._coerce_(P),QQ._coerce_(Q),ZZ(n)))
461    return sage_eval(ans)
462
463def stirling_number1(n,k):
464    """
465    Returns the n-th Stilling number $S_1(n,k)$ of the first kind (the number of
466    permutations of n points with k cycles).
467    Wraps GAP's Stirling1.
468
469    EXAMPLES:
470        sage: stirling_number1(3,2)
471        3
472        sage: stirling_number1(5,2)
473        50
474        sage: 9*stirling_number1(9,5)+stirling_number1(9,4)
475        269325
476        sage: stirling_number1(10,5)
477        269325
478
479    Indeed, $S_1(n,k) = S_1(n-1,k-1) + (n-1)S_1(n-1,k)$.
480    """
481    return ZZ(gap.eval("Stirling1(%s,%s)"%(ZZ(n),ZZ(k))))
482
483def stirling_number2(n,k):
484    """
485    Returns the n-th Stirling number $S_2(n,k)$ of the second kind (the
486    number of ways to partition a set of n elements into k pairwise
487    disjoint nonempty subsets). (The n-th Bell number is the sum of
488    the $S_2(n,k)$'s, $k=0,...,n$.)
489    Wraps GAP's Stirling2.
490
491    EXAMPLES:
492    Stirling numbers satisfy $S_2(n,k) = S_2(n-1,k-1) + kS_2(n-1,k)$:
493
494        sage: 5*stirling_number2(9,5) + stirling_number2(9,4)
495        42525
496        sage: stirling_number2(10,5)
497        42525
498
499        sage: n = stirling_number2(20,11); n
500        1900842429486
501        sage: type(n)
502        <type 'sage.rings.integer.Integer'>
503    """
504    return ZZ(gap.eval("Stirling2(%s,%s)"%(ZZ(n),ZZ(k))))
505
506def mod_stirling(q,n,k):
507    """
508    """
509    if k>n or k<0 or n<0:
510        raise ValueError, "n (= %s) and k (= %s) must greater than or equal to 0, and n must be greater than or equal to k"
511
512    if k == 0:
513        return 1
514    elif k == 1:
515        return (n**2+(2*q+1)*n)/2
516    elif k == n:
517        return prod( [ q+i for i in range(1, n+1) ] )
518    else:
519        return mod_stirling(q,n-1,k)+(q+n)*mod_stirling(q, n-1, k-1)
520
521
522
523
524class CombinatorialObject(SageObject):
525    def __init__(self, l):
526        self.list = l
527
528    def __str__(self):
529        return str(self.list)
530
531    def __repr__(self):
532        return self.list.__repr__()
533
534    def __eq__(self, other):
535        if isinstance(other, CombinatorialObject):
536            return self.list.__eq__(other.list)
537        else:
538            return self.list.__eq__(other)
539
540    def __lt__(self, other):
541        if isinstance(other, CombinatorialObject):
542            return self.list.__lt__(other.list)
543        else:
544            return self.list.__lt__(other)
545
546    def __le__(self, other):
547        if isinstance(other, CombinatorialObject):
548            return self.list.__le__(other.list)
549        else:
550            return self.list.__le__(other)
551
552    def __gt__(self, other):
553        if isinstance(other, CombinatorialObject):
554            return self.list.__gt__(other.list)
555        else:
556            return self.list.__gt__(other)
557
558    def __ge__(self, other):
559        if isinstance(other, CombinatorialObject):
560            return self.list.__ge__(other.list)
561        else:
562            return self.list.__ge__(other)
563
564    def __ne__(self, other):
565        if isinstance(other, CombinatorialObject):
566            return self.list.__ne__(other.list)
567        else:
568            return self.list.__ne__(other)
569
570    def __add__(self, other):
571        return self.list + other
572
573    def __hash__(self):
574        return str(self.list).__hash__()
575
576    #def __cmp__(self, other):
577    #    return self.list.__cmp__(other)
578
579    def __len__(self):
580        return self.list.__len__()
581
582    def __getitem__(self, key):
583        return self.list.__getitem__(key)
584
585    def __iter__(self):
586        return self.list.__iter__()
587
588    def __contains__(self, item):
589        return self.list.__contains__(item)
590
591
592    def index(self, key):
593        return self.list.index(key)
594
595
596class CombinatorialClass(SageObject):
597    def __len__(self):
598        """
599        Returns the number of elements in the combinatorial class.
600
601        EXAMPLES:
602            sage: len(Partitions(5))
603            7
604        """
605        return self.count()
606
607    def __getitem__(self, i):
608        """
609        Returns the combinatorial object of rank i.
610
611        EXAMPLES:
612            sage: p5 = Partitions(5)
613            sage: p5[0]
614            [5]
615            sage: p5[6]
616            [1, 1, 1, 1, 1]
617            sage: p5[7]
618            Traceback (most recent call last):
619            ...
620            ValueError: the value must be between 0 and 6 inclusive
621        """
622        return self.unrank(i)
623
624    def __str__(self):
625        """
626        Returns a string representation of self.
627
628        EXAMPLES:
629            sage: str(Partitions(5))
630            'Partitions of the integer 5'
631        """
632        return self.__repr__()
633
634    def __repr__(self):
635        """
636        EXAMPLES:
637            sage: repr(Partitions(5))
638            'Partitions of the integer 5'
639        """
640        if hasattr(self, '_name') and self._name:
641            return self._name
642        else:
643            return "Combinatorial Class -- REDEFINE ME!"
644
645    def __contains__(self, x):
646        """
647        Tests whether or not the combinatorial class contains the
648        object x.  This raises a NotImplementedError as a default
649        since _all_ subclasses of CombinatorialClass should
650        override this.
651
652        Note that we could replace this with a default implementation
653        that just iterates through the elements of the combinatorial
654        class and checks for equality.  However, since we use __contains__
655        for type checking, this operation should be cheap and should be
656        implemented manually for each combinatorial class.
657        """
658        raise NotImplementedError
659
660    def __iter__(self):
661        """
662        Allows the combinatorial class to be treated as an iterator.
663
664        EXAMPLES:
665            sage: p5 = Partitions(5)
666            sage: [i for i in p5]
667            [[5], [4, 1], [3, 2], [3, 1, 1], [2, 2, 1], [2, 1, 1, 1], [1, 1, 1, 1, 1]]
668        """
669        return self.iterator()
670
671    def __cmp__(self, x):
672        """
673        Compares two different combinatorial classes.  For now, the comparison
674        is done just on their repr's.
675
676        EXAMPLES:
677            sage: p5 = Partitions(5)
678            sage: p6 = Partitions(6)
679            sage: repr(p5) == repr(p6)
680            False
681            sage: p5 == p6
682            False
683        """
684        return cmp(repr(self), repr(x))
685
686    def __count_from_iterator(self):
687        """
688        Default implmentation of count which just goes through the iterator
689        of the combinatorial class to count the number of objects.
690        """
691        c = 0
692        for x in self.iterator():
693            c += 1
694        return c
695    count = __count_from_iterator
696
697    def __call__(self, x):
698        """
699        Returns x as an element of the combinatorial class's object class.
700
701        EXAMPLES:
702            sage: p5 = Partitions(5)
703            sage: a = [2,2,1]
704            sage: type(a)
705            <type 'list'>
706            sage: a = p5(a)
707            sage: type(a)
708            <class 'sage.combinat.partition.Partition_class'>
709            sage: p5([2,1])
710            Traceback (most recent call last):
711            ...
712            ValueError: [2, 1] not in Partitions of the integer 5
713        """
714        if x in self:
715            return self.object_class(x)
716        else:
717            raise ValueError, "%s not in %s"%(x, self)
718
719    def __list_from_iterator(self):
720        """
721        The default implementation of list which builds the list from
722        the iterator.
723        """
724        return [x for x in self.iterator()]
725
726    #Set list to the default implementation
727    list  = __list_from_iterator
728
729    #Set the default object class to be CombinatorialObject
730    object_class = CombinatorialObject
731
732    def __iterator_from_next(self):
733        """
734        An iterator to use when .first() and .next() are provided.
735        """
736        f = self.first()
737        yield f
738        while True:
739            try:
740                f = self.next(f)
741            except:
742                break
743
744            if f == None:
745                break
746            else:
747                yield f
748
749    def __iterator_from_previous(self):
750        """
751        An iterator to use when .last() and .previous() are provided.
752        """
753        l = self.last()
754        yield l
755        while True:
756            try:
757                l = self.previous(l)
758            except:
759                break
760
761            if l == None:
762                break
763            else:
764                yield l
765
766    def __iterator_from_unrank(self):
767        """
768        An iterator to use when .unrank() is provided.
769        """
770        r = 0
771        u = self.unrank(r)
772        yield f
773        while True:
774            r += 1
775            try:
776                u = self.unrank(l)
777            except:
778                break
779
780            if u == None:
781                break
782            else:
783                yield u
784
785    def __iterator_from_list(self):
786        """
787        An iterator to use when .list() is provided()
788        """
789        for x in self.list():
790            yield x
791
792    def iterator(self):
793        """
794        Default implementation of iterator.
795        """
796        #Check to see if .first() and .next() are overridden in the subclass
797        if ( self.first != self.__first_from_iterator and
798             self.next  != self.__next_from_iterator ):
799            return self.__iterator_from_next()
800        #Check to see if .last() and .previous() are overridden in the subclass
801        elif ( self.last != self.__last_from_iterator and
802               self.previous != self.__previous_from_iterator):
803            return self.__iterator_from_previous()
804        #Check to see if .unrank() is overridden in the subclass
805        elif self.unrank != self.__unrank_from_iterator:
806            return self.__iterator_from_unrank()
807        #Finally, check to see if .list() is overridden in the subclass
808        elif self.list != self.__list_from_iterator:
809            return self.__iterator_from_list()
810        else:
811            raise NotImplementedError, "iterator called but not implemented"
812
813
814    def __list_from_unrank_and_count(self):
815        return [ self.unrank(i) for i in range(self.count()) ]
816
817    def __unrank_from_iterator(self, r):
818        """
819        Default implementation of unrank which goes through the iterator.
820        """
821        counter = 0
822        for u in self.iterator():
823            if counter == r:
824                return u
825            counter += 1
826        raise ValueError, "the value must be between %s and %s inclusive"%(0,counter-1)
827
828    #Set the default implementation of unrank
829    unrank = __unrank_from_iterator
830
831
832    def __random_from_unrank(self):
833        """
834        Default implementation of random which uses unrank.
835        """
836        c = self.count()
837        r = randint(0, c-1)
838        if hasattr(self, 'object_class'):
839            return self.object_class(self.unrank(r))
840        else:
841            return self.unrank(r)
842
843    #Set the default implementation of random
844    random = __random_from_unrank
845
846
847    def __rank_from_iterator(self, obj):
848        r = 0
849        for i in self.iterator():
850            if i == obj:
851                return r
852            r += 1
853
854    rank =  __rank_from_iterator
855
856    def __first_from_iterator(self):
857        for i in self.iterator():
858            return i
859
860    first = __first_from_iterator
861
862    def __last_from_iterator(self):
863        for i in self.iterator():
864            pass
865        return i
866
867    last = __last_from_iterator
868
869    def __next_from_iterator(self, obj):
870        if hasattr(obj, 'next'):
871            res = obj.next()
872            if res:
873                return res
874            else:
875                return None
876        found = False
877        for i in self.iterator():
878            if found:
879                return i
880            if i == obj:
881                found = True
882        return None
883
884    next = __next_from_iterator
885
886    def __previous_from_iterator(self, obj):
887        if hasattr(obj, 'previous'):
888            res = obj.previous()
889            if res:
890                return res
891            else:
892                return None
893        prev = None
894        for i in self.iterator():
895            if i == obj:
896                break
897            prev = i
898        return prev
899
900    previous = __previous_from_iterator
901
902    def filter(self, f, name=None):
903        """
904        Returns the combinatorial subclass of f which consists of
905        the elements x of self such that f(x) is True.
906
907        EXAMPLES:
908            sage: P = Permutations(3).filter(lambda x: x.avoids([1,2]))
909            sage: P.list()
910            [[3, 2, 1]]
911        """
912        return FilteredCombinatorialClass(self, f, name=name)
913
914    def union(self, right_cc, name=None):
915        """
916        Returns the combinatorial class representing the union
917        of self and right_cc.
918
919        EXAMPLES:
920            sage: P = Permutations(2).union(Permutations(1))
921            sage: P.list()
922            [[1, 2], [2, 1], [1]]
923        """
924        return UnionCombinatorialClass(self, right_cc, name=name)
925
926class FilteredCombinatorialClass(CombinatorialClass):
927    def __init__(self, combinatorial_class, f, name=None):
928        """
929        TESTS:
930            sage: P = Permutations(3).filter(lambda x: x.avoids([1,2]))
931        """
932        self.f = f
933        self.combinatorial_class = combinatorial_class
934        self._name = name
935
936    def __repr__(self):
937        if self._name:
938            return self._name
939        else:
940            return "Filtered sublass of " + repr(self.combinatorial_class)
941
942    def __contains__(self, x):
943        return x in self.combinatorial_class and self.f(x)
944
945    def count(self):
946        c = 0
947        for x in self.iterator():
948            c += 1
949        return c
950
951    def list(self):
952        res = []
953        for x in self.combinatorial_class.iterator():
954            if self.f(x):
955                res.append(x)
956        return res
957
958    def iterator(self):
959        for x in self.combinatorial_class.iterator():
960            if self.f(x):
961                yield x
962
963    def __unrank_from_iterator(self, r):
964        """
965        Default implementation of unrank which goes through the iterator.
966        """
967        counter = 0
968        for u in self.iterator():
969            if counter == r:
970                return u
971            counter += 1
972        raise ValueError, "the value must be between %s and %s inclusive"%(0,counter-1)
973
974    #Set the default implementation of unrank
975    unrank = __unrank_from_iterator
976
977
978    def __random_from_unrank(self):
979        """
980        Default implementation of random which uses unrank.
981        """
982        c = self.count()
983        r = randint(0, c-1)
984        if hasattr(self, 'object_class'):
985            return self.object_class(self.unrank(r))
986        else:
987            return self.unrank(r)
988
989    #Set the default implementation of random
990    random = __random_from_unrank
991
992
993    def __rank_from_iterator(self, obj):
994        r = 0
995        for i in self.iterator():
996            if i == obj:
997                return r
998            r += 1
999
1000    rank =  __rank_from_iterator
1001
1002    def __first_from_iterator(self):
1003        for i in self.iterator():
1004            return i
1005
1006    first = __first_from_iterator
1007
1008    def __last_from_iterator(self):
1009        for i in self.iterator():
1010            pass
1011        return i
1012
1013    last = __last_from_iterator
1014
1015    def __next_from_iterator(self, obj):
1016        if hasattr(obj, 'next'):
1017            res = obj.next()
1018            if res:
1019                return res
1020            else:
1021                return None
1022        found = False
1023        for i in self.iterator():
1024            if found:
1025                return i
1026            if i == obj:
1027                found = True
1028        return None
1029
1030    next = __next_from_iterator
1031
1032    def __previous_from_iterator(self, obj):
1033        if hasattr(obj, 'previous'):
1034            res = obj.previous()
1035            if res:
1036                return res
1037            else:
1038                return None
1039        prev = None
1040        for i in self.iterator():
1041            if i == obj:
1042                break
1043            prev = i
1044        return prev
1045
1046    previous = __previous_from_iterator
1047
1048
1049class UnionCombinatorialClass(CombinatorialClass):
1050    def __init__(self, left_cc, right_cc, name=None):
1051        """
1052        TESTS:
1053            sage: P = Permutations(3).union(Permutations(2))
1054            sage: P == loads(dumps(P))
1055            True
1056        """
1057        self.left_cc = left_cc
1058        self.right_cc = right_cc
1059        self._name = name
1060
1061    def __repr__(self):
1062        """
1063        TESTS:
1064            sage: print repr(Permutations(3).union(Permutations(2)))
1065            Union combinatorial class of
1066                Standard permutations of 3
1067            and
1068                Standard permutations of 2
1069
1070        """
1071        if self._name:
1072            return self._name
1073        else:
1074            return "Union combinatorial class of \n    %s\nand\n    %s"%(self.left_cc, self.right_cc)
1075
1076    def __contains__(self, x):
1077        return x in self.left_cc or x in self.right_cc
1078
1079    def count(self):
1080        return self.left_cc.count() + self.right_cc.count()
1081
1082    def list(self):
1083        res = []
1084        for x in self.left_cc.iterator():
1085            res.append(x)
1086        for x in self.right_cc.iterator():
1087            res.append(x)
1088        return res
1089
1090    def iterator(self):
1091        for x in self.left_cc.iterator():
1092            yield x
1093        for x in self.right_cc.iterator():
1094            yield x
1095
1096    def __unrank_from_iterator(self, r):
1097        """
1098        Default implementation of unrank which goes through the iterator.
1099        """
1100        counter = 0
1101        for u in self.iterator():
1102            if counter == r:
1103                return u
1104            counter += 1
1105        raise ValueError, "the value must be between %s and %s inclusive"%(0,counter-1)
1106
1107    #Set the default implementation of unrank
1108    unrank = __unrank_from_iterator
1109
1110    def __random_from_unrank(self):
1111        """
1112        Default implementation of random which uses unrank.
1113        """
1114        c = self.count()
1115        r = randint(0, c-1)
1116        if hasattr(self, 'object_class'):
1117            return self.object_class(self.unrank(r))
1118        else:
1119            return self.unrank(r)
1120
1121    #Set the default implementation of random
1122    random = __random_from_unrank
1123
1124    def __rank_from_iterator(self, obj):
1125        r = 0
1126        for i in self.iterator():
1127            if i == obj:
1128                return r
1129            r += 1
1130
1131    rank =  __rank_from_iterator
1132
1133    def __first_from_iterator(self):
1134        for i in self.iterator():
1135            return i
1136
1137    first = __first_from_iterator
1138
1139    def __last_from_iterator(self):
1140        for i in self.iterator():
1141            pass
1142        return i
1143
1144    last = __last_from_iterator
1145
1146    def __next_from_iterator(self, obj):
1147        if hasattr(obj, 'next'):
1148            res = obj.next()
1149            if res:
1150                return res
1151            else:
1152                return None
1153        found = False
1154        for i in self.iterator():
1155            if found:
1156                return i
1157            if i == obj:
1158                found = True
1159        return None
1160
1161    next = __next_from_iterator
1162
1163    def __previous_from_iterator(self, obj):
1164        if hasattr(obj, 'previous'):
1165            res = obj.previous()
1166            if res:
1167                return res
1168            else:
1169                return None
1170        prev = None
1171        for i in self.iterator():
1172            if i == obj:
1173                break
1174            prev = i
1175        return prev
1176
1177    previous = __previous_from_iterator
1178
1179
1180def hurwitz_zeta(s,x,N):
1181    """
1182    Returns the value of the $\zeta(s,x)$ to $N$ decimals, where s and x are real.
1183
1184    The Hurwitz zeta function is one of the many zeta functions. It defined as
1185    $1186 \zeta(s,x) = \sum_{k=0}^\infty (k+x)^{-s}. 1187$
1188    When $x = 1$, this coincides with Riemann's zeta function. The Dirichlet L-functions
1189    may be expressed as a linear combination of Hurwitz zeta functions.
1190
1191    EXAMPLES:
1192        sage: hurwitz_zeta(3,1/2,6)
1193        8.41439000000000
1194        sage: hurwitz_zeta(1.1,1/2,6)
1195        12.1041000000000
1196        sage: hurwitz_zeta(1.1,1/2,50)
1197        12.103813495683744469025853545548130581952676591199
1198
1199    REFERENCES:
1200        http://en.wikipedia.org/wiki/Hurwitz_zeta_function
1201
1202    """
1203    maxima.eval('load ("bffac")')
1204    s = maxima.eval("bfhzeta (%s,%s,%s)"%(s,x,N))
1205
1206    #Handle the case where there is a 'b' in the string
1207    #'1.2000b0' means 1.2000 and
1208    #'1.2000b1' means 12.000
1209    i = s.rfind('b')
1210    if i == -1:
1211        return sage_eval(s)
1212    else:
1213        if s[i+1:] == '0':
1214            return sage_eval(s[:i])
1215        else:
1216            return sage_eval(s[:i])*10**sage_eval(s[i+1:])
1217
1218    return s  ## returns an odd string
1219
1220
1221#####################################################
1222#### combinatorial sets/lists
1223
1224def combinations(mset,k):
1225    r"""
1226    A {\it combination} of a multiset (a list of objects which may
1227    contain the same object several times) mset is an unordered
1228    selection without repetitions and is represented by a sorted
1229    sublist of mset.  Returns the set of all combinations of the
1230    multiset mset with k elements.
1231
1232    WARNING: Wraps GAP's Combinations.  Hence mset must be a list of
1233    objects that have string representations that can be interpreted by
1234    the GAP intepreter.  If mset consists of at all complicated SAGE
1235    objects, this function does *not* do what you expect.  A proper
1236    function should be written! (TODO!)
1237
1238    EXAMPLES:
1239        sage: mset = [1,1,2,3,4,4,5]
1240        sage: combinations(mset,2)
1241        [[1, 1],
1242         [1, 2],
1243         [1, 3],
1244         [1, 4],
1245         [1, 5],
1246         [2, 3],
1247         [2, 4],
1248         [2, 5],
1249         [3, 4],
1250         [3, 5],
1251         [4, 4],
1252         [4, 5]]
1253         sage: mset = ["d","a","v","i","d"]
1254         sage: combinations(mset,3)
1255         ['add', 'adi', 'adv', 'aiv', 'ddi', 'ddv', 'div']
1256
1257    NOTE: For large lists, this raises a string error.
1258    """
1259    ans=gap.eval("Combinations(%s,%s)"%(mset,ZZ(k))).replace("\n","")
1260    return eval(ans)
1261
1262def combinations_iterator(mset,n=None):
1263    """
1264    Posted by Raymond Hettinger, 2006/03/23, to the Python Cookbook:
1265    http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/474124
1266
1267    Much faster than combinations.
1268
1269    EXAMPLES:
1270        sage: X = combinations_iterator([1,2,3,4,5],3)
1271        sage: [x for x in X]
1272        [[1, 2, 3],
1273         [1, 2, 4],
1274         [1, 2, 5],
1275         [1, 3, 4],
1276         [1, 3, 5],
1277         [1, 4, 5],
1278         [2, 3, 4],
1279         [2, 3, 5],
1280         [2, 4, 5],
1281         [3, 4, 5]]
1282    """
1283    items = mset
1284    if n is None:
1285        n = len(items)
1286    for i in range(len(items)):
1287        v = items[i:i+1]
1288        if n == 1:
1289            yield v
1290        else:
1291            rest = items[i+1:]
1292            for c in combinations_iterator(rest, n-1):
1293                yield v + c
1294
1295
1296def number_of_combinations(mset,k):
1297    """
1298    Returns the size of combinations(mset,k).
1299    IMPLEMENTATION: Wraps GAP's NrCombinations.
1300
1301
1302    NOTE: mset must be a list of integers or strings (i.e., this is very restricted).
1303
1304    EXAMPLES:
1305        sage: mset = [1,1,2,3,4,4,5]
1306        sage: number_of_combinations(mset,2)
1307        12
1308    """
1309    return ZZ(gap.eval("NrCombinations(%s,%s)"%(mset,ZZ(k))))
1310
1311def arrangements(mset,k):
1312    r"""
1313    An arrangement of mset is an ordered selection without repetitions
1314    and is represented by a list that contains only elements from
1315    mset, but maybe in a different order.
1316
1317    \code{arrangements} returns the set of arrangements of the
1318    multiset mset that contain k elements.
1319
1320    IMPLEMENTATION: Wraps GAP's Arrangements.
1321
1322    WARNING: Wraps GAP -- hence mset must be a list of objects that
1323    have string representations that can be interpreted by the GAP
1324    intepreter.  If mset consists of at all complicated SAGE objects,
1325    this function does *not* do what you expect.  A proper function
1326    should be written! (TODO!)
1327
1328    EXAMPLES:
1329        sage: mset = [1,1,2,3,4,4,5]
1330        sage: arrangements(mset,2)
1331        [[1, 1],
1332         [1, 2],
1333         [1, 3],
1334         [1, 4],
1335         [1, 5],
1336         [2, 1],
1337         [2, 3],
1338         [2, 4],
1339         [2, 5],
1340         [3, 1],
1341         [3, 2],
1342         [3, 4],
1343         [3, 5],
1344         [4, 1],
1345         [4, 2],
1346         [4, 3],
1347         [4, 4],
1348         [4, 5],
1349         [5, 1],
1350         [5, 2],
1351         [5, 3],
1352         [5, 4]]
1353         sage: arrangements( ["c","a","t"], 2 )
1354         ['ac', 'at', 'ca', 'ct', 'ta', 'tc']
1355         sage: arrangements( ["c","a","t"], 3 )
1356         ['act', 'atc', 'cat', 'cta', 'tac', 'tca']
1357    """
1358    ans=gap.eval("Arrangements(%s,%s)"%(mset,k))
1359    return eval(ans)
1360
1361def number_of_arrangements(mset,k):
1362    """
1363    Returns the size of arrangements(mset,k).
1364    Wraps GAP's NrArrangements.
1365
1366    EXAMPLES:
1367        sage: mset = [1,1,2,3,4,4,5]
1368        sage: number_of_arrangements(mset,2)
1369        22
1370    """
1371    return ZZ(gap.eval("NrArrangements(%s,%s)"%(mset,ZZ(k))))
1372
1373def derangements(mset):
1374    """
1375    A derangement is a fixed point free permutation of list and is
1376    represented by a list that contains exactly the same elements as
1377    mset, but possibly in different order.  Derangements returns the
1378    set of all derangements of a multiset.
1379
1380    Wraps GAP's Derangements.
1381
1382    WARNING: Wraps GAP -- hence mset must be a list of objects that
1383    have string representations that can be interpreted by the GAP
1384    intepreter.  If mset consists of at all complicated SAGE objects,
1385    this function does *not* do what you expect.  A proper function
1386    should be written! (TODO!)
1387
1388    EXAMPLES:
1389        sage: mset = [1,2,3,4]
1390        sage: derangements(mset)
1391        [[2, 1, 4, 3],
1392         [2, 3, 4, 1],
1393         [2, 4, 1, 3],
1394         [3, 1, 4, 2],
1395         [3, 4, 1, 2],
1396         [3, 4, 2, 1],
1397         [4, 1, 2, 3],
1398         [4, 3, 1, 2],
1399         [4, 3, 2, 1]]
1400         sage: derangements(["c","a","t"])
1401         ['atc', 'tca']
1402
1403    """
1404    ans=gap.eval("Derangements(%s)"%mset)
1405    return eval(ans)
1406
1407def number_of_derangements(mset):
1408    """
1409    Returns the size of derangements(mset).
1410    Wraps GAP's NrDerangements.
1411
1412    EXAMPLES:
1413        sage: mset = [1,2,3,4]
1414        sage: number_of_derangements(mset)
1415        9
1416    """
1417    ans=gap.eval("NrDerangements(%s)"%mset)
1418    return ZZ(ans)
1419
1420def tuples(S,k):
1421    """
1422    An ordered tuple of length k of set is an ordered selection with
1423    repetition and is represented by a list of length k containing
1424    elements of set.
1425    tuples returns the set of all ordered tuples of length k of the set.
1426
1427    EXAMPLES:
1428        sage: S = [1,2]
1429        sage: tuples(S,3)
1430        [[1, 1, 1], [2, 1, 1], [1, 2, 1], [2, 2, 1], [1, 1, 2], [2, 1, 2], [1, 2, 2], [2, 2, 2]]
1431        sage: mset = ["s","t","e","i","n"]
1432        sage: tuples(mset,2)
1433        [['s', 's'], ['t', 's'], ['e', 's'], ['i', 's'], ['n', 's'], ['s', 't'], ['t', 't'],
1434         ['e', 't'], ['i', 't'], ['n', 't'], ['s', 'e'], ['t', 'e'], ['e', 'e'], ['i', 'e'],
1435         ['n', 'e'], ['s', 'i'], ['t', 'i'], ['e', 'i'], ['i', 'i'], ['n', 'i'], ['s', 'n'],
1436         ['t', 'n'], ['e', 'n'], ['i', 'n'], ['n', 'n']]
1437
1438    The Set(...) comparisons are necessary because finite fields are not
1439    enumerated in a standard order.
1440        sage: K.<a> = GF(4, 'a')
1441        sage: mset = [x for x in K if x!=0]
1442        sage: tuples(mset,2)
1443        [[a, a], [a + 1, a], [1, a], [a, a + 1], [a + 1, a + 1], [1, a + 1], [a, 1], [a + 1, 1], [1, 1]]
1444
1445    AUTHOR: Jon Hanke (2006-08?)
1446    """
1447    import copy
1448    if k<=0:
1449        return [[]]
1450    if k==1:
1451        return [[x] for x in S]
1452    ans = []
1453    for s in S:
1454        for x in tuples(S,k-1):
1455            y = copy.copy(x)
1456            y.append(s)
1457            ans.append(y)
1458    return ans
1459    ## code wrapping GAP's Tuples:
1460    #ans=gap.eval("Tuples(%s,%s)"%(S,k))
1461    #return eval(ans)
1462
1463
1464def number_of_tuples(S,k):
1465    """
1466    Returns the size of tuples(S,k).
1467    Wraps GAP's NrTuples.
1468
1469    EXAMPLES:
1470        sage: S = [1,2,3,4,5]
1471        sage: number_of_tuples(S,2)
1472        25
1473        sage: S = [1,1,2,3,4,5]
1474        sage: number_of_tuples(S,2)
1475        25
1476
1477    """
1478    ans=gap.eval("NrTuples(%s,%s)"%(S,ZZ(k)))
1479    return ZZ(ans)
1480
1481def unordered_tuples(S,k):
1482    """
1483    An unordered tuple of length k of set is a unordered selection
1484    with repetitions of set and is represented by a sorted list of
1485    length k containing elements from set.
1486
1487    unordered_tuples returns the set of all unordered tuples of length k
1488    of the set.
1489    Wraps GAP's UnorderedTuples.
1490
1491    WARNING: Wraps GAP -- hence mset must be a list of objects that
1492    have string representations that can be interpreted by the GAP
1493    intepreter.  If mset consists of at all complicated SAGE objects,
1494    this function does *not* do what you expect.  A proper function
1495    should be written! (TODO!)
1496
1497    EXAMPLES:
1498        sage: S = [1,2]
1499        sage: unordered_tuples(S,3)
1500        [[1, 1, 1], [1, 1, 2], [1, 2, 2], [2, 2, 2]]
1501        sage: unordered_tuples(["a","b","c"],2)
1502        ['aa', 'ab', 'ac', 'bb', 'bc', 'cc']
1503
1504    """
1505    ans=gap.eval("UnorderedTuples(%s,%s)"%(S,ZZ(k)))
1506    return eval(ans)
1507
1508def number_of_unordered_tuples(S,k):
1509    """
1510    Returns the size of unordered_tuples(S,k).
1511    Wraps GAP's NrUnorderedTuples.
1512
1513    EXAMPLES:
1514        sage: S = [1,2,3,4,5]
1515        sage: number_of_unordered_tuples(S,2)
1516        15
1517    """
1518    ans=gap.eval("NrUnorderedTuples(%s,%s)"%(S,ZZ(k)))
1519    return ZZ(ans)
1520
1521def permutations(mset):
1522    """
1523    A {\it permutation} is represented by a list that contains exactly the same
1524    elements as mset, but possibly in different order. If mset is a
1525    proper set there are $|mset| !$ such permutations. Otherwise if the
1526    first elements appears $k_1$ times, the second element appears $k_2$ times
1527    and so on, the number of permutations is $|mset|! / (k_1! k_2! \ldots)$,
1528    which is sometimes called a {\it multinomial coefficient}.
1529
1530    permutations returns the set of all permutations of a multiset.
1531    Wraps GAP's PermutationsList.
1532
1533    WARNING: Wraps GAP -- hence mset must be a list of objects that
1534    have string representations that can be interpreted by the GAP
1535    intepreter.  If mset consists of at all complicated SAGE objects,
1536    this function does *not* do what you expect.  A proper function
1537    should be written! (TODO!)
1538
1539    EXAMPLES:
1540        sage: mset = [1,1,2,2,2]
1541        sage: permutations(mset)
1542        [[1, 1, 2, 2, 2],
1543         [1, 2, 1, 2, 2],
1544         [1, 2, 2, 1, 2],
1545         [1, 2, 2, 2, 1],
1546         [2, 1, 1, 2, 2],
1547         [2, 1, 2, 1, 2],
1548         [2, 1, 2, 2, 1],
1549         [2, 2, 1, 1, 2],
1550         [2, 2, 1, 2, 1],
1551         [2, 2, 2, 1, 1]]
1552
1553    """
1554    ans=gap.eval("PermutationsList(%s)"%mset)
1555    return eval(ans)
1556
1557def permutations_iterator(mset,n=None):
1558    """
1559    Posted by Raymond Hettinger, 2006/03/23, to the Python Cookbook:
1560    http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/474124
1561
1562    Note-- This function considers repeated elements as different entries,
1563    so for example:
1564        sage: from sage.combinat.combinat import permutations, permutations_iterator
1565        sage: mset = [1,2,2]
1566        sage: permutations(mset)
1567        [[1, 2, 2], [2, 1, 2], [2, 2, 1]]
1568        sage: for p in permutations_iterator(mset): print p
1569        [1, 2, 2]
1570        [1, 2, 2]
1571        [2, 1, 2]
1572        [2, 2, 1]
1573        [2, 1, 2]
1574        [2, 2, 1]
1575
1576
1577    EXAMPLES:
1578        sage: X = permutations_iterator(range(3),2)
1579        sage: [x for x in X]
1580        [[0, 1], [0, 2], [1, 0], [1, 2], [2, 0], [2, 1]]
1581    """
1582    items = mset
1583    if n is None:
1584        n = len(items)
1585    for i in range(len(items)):
1586        v = items[i:i+1]
1587        if n == 1:
1588            yield v
1589        else:
1590            rest = items[:i] + items[i+1:]
1591            for p in permutations_iterator(rest, n-1):
1592                yield v + p
1593
1594def number_of_permutations(mset):
1595    """
1596    Returns the size of permutations(mset).
1597
1598    AUTHOR: Robert L. Miller
1599
1600    EXAMPLES:
1601        sage: mset = [1,1,2,2,2]
1602        sage: number_of_permutations(mset)
1603        10
1604
1605    """
1606    from sage.rings.arith import factorial, prod
1607    m = len(mset)
1608    n = []
1609    seen = []
1610    for element in mset:
1611        try:
1612            n[seen.index(element)] += 1
1613        except:
1614            n.append(1)
1615            seen.append(element)
1616    return factorial(m)/prod([factorial(k) for k in n])
1617
1618def cyclic_permutations(mset):
1619    """
1620    Returns a list of all cyclic permutations of mset. Treats mset as a list,
1621    not a set, i.e. entries with the same value are distinct.
1622
1623    AUTHOR: Emily Kirkman
1624
1625    EXAMPLES:
1626        sage: from sage.combinat.combinat import cyclic_permutations, cyclic_permutations_iterator
1627        sage: cyclic_permutations(range(4))
1628        [[0, 1, 2, 3], [0, 1, 3, 2], [0, 2, 1, 3], [0, 2, 3, 1], [0, 3, 1, 2], [0, 3, 2, 1]]
1629        sage: for cycle in cyclic_permutations(['a', 'b', 'c']):
1630        ...       print cycle
1631        ['a', 'b', 'c']
1632        ['a', 'c', 'b']
1633
1634    Note that lists with repeats are not handled intuitively:
1635        sage: cyclic_permutations([1,1,1])
1636        [[1, 1, 1], [1, 1, 1]]
1637
1638    """
1639    return list(cyclic_permutations_iterator(mset))
1640
1641def cyclic_permutations_iterator(mset):
1642    """
1643    Iterates over all cyclic permutations of mset in cycle notation. Treats
1644    mset as a list, not a set, i.e. entries with the same value are distinct.
1645
1646    AUTHOR: Emily Kirkman
1647
1648    EXAMPLES:
1649        sage: from sage.combinat.combinat import cyclic_permutations, cyclic_permutations_iterator
1650        sage: cyclic_permutations(range(4))
1651        [[0, 1, 2, 3], [0, 1, 3, 2], [0, 2, 1, 3], [0, 2, 3, 1], [0, 3, 1, 2], [0, 3, 2, 1]]
1652        sage: for cycle in cyclic_permutations(['a', 'b', 'c']):
1653        ...       print cycle
1654        ['a', 'b', 'c']
1655        ['a', 'c', 'b']
1656
1657    Note that lists with repeats are not handled intuitively:
1658        sage: cyclic_permutations([1,1,1])
1659        [[1, 1, 1], [1, 1, 1]]
1660
1661    """
1662    if len(mset) > 2:
1663        for perm in permutations_iterator(mset[1:]):
1664            yield [mset[0]] + perm
1665    else:
1666        yield mset
1667
1668#### partitions
1669
1670def partitions_set(S,k=None, use_file=True):
1671    r"""
1672    An {\it unordered partition} of a set $S$ is a set of pairwise disjoint
1673    nonempty subsets with union $S$ and is represented by a sorted
1674    list of such subsets.
1675
1676    partitions_set returns the set of all unordered partitions of the
1677    list $S$ of increasing positive integers into k pairwise disjoint
1678    nonempty sets. If k is omitted then all partitions are returned.
1679
1680    The Bell number $B_n$, named in honor of Eric Temple Bell, is
1681    the number of different partitions of a set with n elements.
1682
1683    WARNING: Wraps GAP -- hence S must be a list of objects that have
1684    string representations that can be interpreted by the GAP
1685    intepreter.  If mset consists of at all complicated SAGE objects,
1686    this function does *not* do what you expect.  A proper function
1687    should be written! (TODO!)
1688
1689    WARNING: This function is inefficient.  The runtime is dominated
1690    by parsing the output from GAP.
1691
1692    Wraps GAP's PartitionsSet.
1693
1694    EXAMPLES:
1695        sage: S = [1,2,3,4]
1696        sage: partitions_set(S,2)
1697        [[[1], [2, 3, 4]],
1698         [[1, 2], [3, 4]],
1699         [[1, 2, 3], [4]],
1700         [[1, 2, 4], [3]],
1701         [[1, 3], [2, 4]],
1702         [[1, 3, 4], [2]],
1703         [[1, 4], [2, 3]]]
1704
1705    REFERENCES:
1706       http://en.wikipedia.org/wiki/Partition_of_a_set
1707    """
1708    if k is None:
1709        ans=gap("PartitionsSet(%s)"%S).str(use_file=use_file)
1710    else:
1711        ans=gap("PartitionsSet(%s,%s)"%(S,k)).str(use_file=use_file)
1712    return eval(ans)
1713
1714def number_of_partitions_set(S,k):
1715    r"""
1716    Returns the size of \code{partitions_set(S,k)}.  Wraps GAP's
1717    NrPartitionsSet.
1718
1719    The Stirling number of the second kind is the number of partitions
1720    of a set of size n into k blocks.
1721
1722    EXAMPLES:
1723        sage: mset = [1,2,3,4]
1724        sage: number_of_partitions_set(mset,2)
1725        7
1726        sage: stirling_number2(4,2)
1727        7
1728
1729    REFERENCES
1730        http://en.wikipedia.org/wiki/Partition_of_a_set
1731
1732    """
1733    if k is None:
1734        ans=gap.eval("NrPartitionsSet(%s)"%S)
1735    else:
1736        ans=gap.eval("NrPartitionsSet(%s,%s)"%(S,ZZ(k)))
1737    return ZZ(ans)
1738
1739def partitions_list(n,k=None):
1740    r"""
1741    An {\it unordered partition of $n$} is an unordered sum
1742    $n = p_1+p_2 +\ldots+ p_k$ of positive integers and is represented by
1743    the list $p = [p_1,p_2,\ldots,p_k]$, in nonincreasing order, i.e.,
1744    $p1\geq p_2 ...\geq p_k$.
1745
1746    INPUT:
1747        n -- a positive integer
1748
1749    \code{partitions_list(n,k)} returns the list of all (unordered)
1750    partitions of the positive integer n into sums with k summands. If
1751    k is omitted then all partitions are returned.
1752
1753    Do not call partitions_list with an n much larger than 40, in
1754    which case there are 37338 partitions, since the list will simply
1755    become too large.
1756
1757    Wraps GAP's Partitions.
1758
1759    The function \code{partitions} (a wrapper for the corresponding
1760    PARI function) returns not a list but rather a generator for a
1761    list. It is also a function of only one argument.
1762
1763    EXAMPLES:
1764        sage: partitions_list(10,2)
1765        [[5, 5], [6, 4], [7, 3], [8, 2], [9, 1]]
1766        sage: partitions_list(5)
1767        [[1, 1, 1, 1, 1], [2, 1, 1, 1], [2, 2, 1], [3, 1, 1], [3, 2], [4, 1], [5]]
1768
1769    However, partitions(5) returns <generator object at ...>''.
1770    """
1771    n = ZZ(n)
1772    if n <= 0:
1773        raise ValueError, "n (=%s) must be a positive integer"%n
1774    if k is None:
1775        ans=gap.eval("Partitions(%s)"%(n))
1776    else:
1777        ans=gap.eval("Partitions(%s,%s)"%(n,k))
1778    return eval(ans.replace('\n',''))
1779
1780def number_of_partitions(n,k=None, algorithm='default'):
1781    r"""
1782    Returns the size of partitions_list(n,k).
1783
1784    INPUT:
1785        n -- an integer
1786        k -- (default: None); if specified, instead returns the
1787             cardinality of the set of all (unordered) partitions of
1788             the positive integer n into sums with k summands.
1789        algorithm -- (default: 'default')
1790            'default' -- If k is not None, then use Gap (very slow).
1791                         If k is None, use Jon Bober's highly
1792                         optimized implementation (this is the fastest
1793                         code in the world for this problem).
1794            'bober' -- use Jonathon Bober's implementation
1795            'gap' -- use GAP (VERY *slow*)
1796            'pari' -- use PARI.  Speed seems the same as GAP until $n$ is
1797                      in the thousands, in which case PARI is faster. *But*
1798                      PARI has a bug, e.g., on 64-bit Linux PARI-2.3.2
1799                      outputs numbpart(147007)%1000 as 536, but it
1800                      should be 533!.  So do not use this option.
1801
1802    IMPLEMENTATION: Wraps GAP's NrPartitions or PARI's numbpart function.
1803
1804    Use the function \code{partitions(n)} to return a generator over
1805    all partitions of $n$.
1806
1807    It is possible to associate with every partition of the integer n
1808    a conjugacy class of permutations in the symmetric group on n
1809    points and vice versa.  Therefore p(n) = NrPartitions(n) is the
1810    number of conjugacy classes of the symmetric group on n points.
1811
1812    EXAMPLES:
1813        sage: v = list(partitions(5)); v
1814        [(1, 1, 1, 1, 1), (1, 1, 1, 2), (1, 2, 2), (1, 1, 3), (2, 3), (1, 4), (5,)]
1815        sage: len(v)
1816        7
1817        sage: number_of_partitions(5, algorithm='gap')
1818        7
1819        sage: number_of_partitions(5, algorithm='pari')
1820        7
1821        sage: number_of_partitions(5, algorithm='bober')
1822        7
1823
1824    The input must be a nonnegative integer or a ValueError is raised.
1825        sage: number_of_partitions(-5)
1826        Traceback (most recent call last):
1827        ...
1828        ValueError: n (=-5) must be a nonnegative integer
1829
1830        sage: number_of_partitions(10,2)
1831        5
1832        sage: number_of_partitions(10)
1833        42
1834        sage: number_of_partitions(3)
1835        3
1836        sage: number_of_partitions(10)
1837        42
1838        sage: number_of_partitions(3, algorithm='pari')
1839        3
1840        sage: number_of_partitions(10, algorithm='pari')
1841        42
1842        sage: number_of_partitions(40)
1843        37338
1844        sage: number_of_partitions(100)
1845        190569292
1846        sage: number_of_partitions(100000)
1847        27493510569775696512677516320986352688173429315980054758203125984302147328114964173055050741660736621590157844774296248940493063070200461792764493033510116079342457190155718943509725312466108452006369558934464248716828789832182345009262853831404597021307130674510624419227311238999702284408609370935531629697851569569892196108480158600569421098519
1848
1849    A generating function for p(n) is given by the reciprocal of
1850    Euler's function:
1851
1852    $1853 \sum_{n=0}^\infty p(n)x^n = \prod_{k=1}^\infty \left(\frac {1}{1-x^k} \right). 1854$
1855
1856    We use SAGE to verify that the first several coefficients do
1857    instead agree:
1858
1859        sage: q = PowerSeriesRing(QQ, 'q', default_prec=9).gen()
1860        sage: prod([(1-q^k)^(-1) for k in range(1,9)])  ## partial product of
1861        1 + q + 2*q^2 + 3*q^3 + 5*q^4 + 7*q^5 + 11*q^6 + 15*q^7 + 22*q^8 + O(q^9)
1862        sage: [number_of_partitions(k) for k in range(2,10)]
1863        [2, 3, 5, 7, 11, 15, 22, 30]
1864
1865    REFERENCES:
1866        http://en.wikipedia.org/wiki/Partition_%28number_theory%29
1867
1868    TESTS:
1869        sage: n = 500 + randint(0,500)
1870        sage: number_of_partitions( n - (n % 385) + 369) % 385 == 0
1871        True
1872        sage: n = 1500 + randint(0,1500)
1873        sage: number_of_partitions( n - (n % 385) + 369) % 385 == 0
1874        True
1875        sage: n = 1000000 + randint(0,1000000)
1876        sage: number_of_partitions( n - (n % 385) + 369) % 385 == 0
1877        True
1878        sage: n = 1000000 + randint(0,1000000)
1879        sage: number_of_partitions( n - (n % 385) + 369) % 385 == 0
1880        True
1881        sage: n = 1000000 + randint(0,1000000)
1882        sage: number_of_partitions( n - (n % 385) + 369) % 385 == 0
1883        True
1884        sage: n = 1000000 + randint(0,1000000)
1885        sage: number_of_partitions( n - (n % 385) + 369) % 385 == 0
1886        True
1887        sage: n = 1000000 + randint(0,1000000)
1888        sage: number_of_partitions( n - (n % 385) + 369) % 385 == 0
1889        True
1890        sage: n = 1000000 + randint(0,1000000)
1891        sage: number_of_partitions( n - (n % 385) + 369) % 385 == 0
1892        True
1893        sage: n = 100000000 + randint(0,100000000)
1894        sage: number_of_partitions( n - (n % 385) + 369) % 385 == 0
1895        True
1896        sage: n = 1000000000 + randint(0,1000000000)
1897        sage: number_of_partitions( n - (n % 385) + 369) % 385 == 0      # takes a long time
1898        True
1899
1900    Another consistency test for n up to 500:
1901        sage: len([n for n in [1..500] if number_of_partitions(n) != number_of_partitions(n,algorithm='pari')])
1902        0
1903    """
1904    n = ZZ(n)
1905    if n < 0:
1906        raise ValueError, "n (=%s) must be a nonnegative integer"%n
1907    elif n == 0:
1908        return ZZ(1)
1909
1910    if algorithm == 'default':
1911        if k is None:
1912            algorithm = 'bober'
1913        else:
1914            algorithm = 'gap'
1915
1916    if algorithm == 'gap':
1917        if k is None:
1918            ans=gap.eval("NrPartitions(%s)"%(ZZ(n)))
1919        else:
1920            ans=gap.eval("NrPartitions(%s,%s)"%(ZZ(n),ZZ(k)))
1921        return ZZ(ans)
1922
1923    if k is not None:
1924        raise ValueError, "only the GAP algorithm works if k is specified."
1925
1926    if algorithm == 'bober':
1927        return partitions_ext.number_of_partitions(n)
1928
1929    elif algorithm == 'pari':
1930        return ZZ(pari(ZZ(n)).numbpart())
1931
1932    raise ValueError, "unknown algorithm '%s'"%algorithm
1933
1934def partitions(n):
1935    r"""
1936    Generator of all the partitions of the integer $n$.
1937
1938    INPUT:
1939        n -- int
1940
1941    To compute the number of partitions of $n$ use
1942    \code{number_of_partitions(n)}.
1943
1944    EXAMPLES:
1945        sage: partitions(3)          # random location
1946        <generator object at 0xab3b3eac>
1947        sage: list(partitions(3))
1948        [(1, 1, 1), (1, 2), (3,)]
1949
1950
1951    AUTHOR: Adapted from David Eppstein, Jan Van lent, George Yoshida;
1952    Python Cookbook 2, Recipe 19.16.
1953    """
1954    n == ZZ(n)
1955    # base case of the recursion: zero is the sum of the empty tuple
1956    if n == 0:
1957        yield ( )
1958        return
1959    # modify the partitions of n-1 to form the partitions of n
1960    for p in partitions(n-1):
1961        yield (1,) + p
1962        if p and (len(p) < 2 or p[1] > p[0]):
1963            yield (p[0] + 1,) + p[1:]
1964
1965def cyclic_permutations_of_partition(partition):
1966    """
1967    Returns all combinations of cyclic permutations of each cell of the
1968    partition.
1969
1970    AUTHOR: Robert L. Miller
1971
1972    EXAMPLES:
1973        sage: from sage.combinat.combinat import cyclic_permutations_of_partition
1974        sage: cyclic_permutations_of_partition([[1,2,3,4],[5,6,7]])
1975        [[[1, 2, 3, 4], [5, 6, 7]],
1976         [[1, 2, 4, 3], [5, 6, 7]],
1977         [[1, 3, 2, 4], [5, 6, 7]],
1978         [[1, 3, 4, 2], [5, 6, 7]],
1979         [[1, 4, 2, 3], [5, 6, 7]],
1980         [[1, 4, 3, 2], [5, 6, 7]],
1981         [[1, 2, 3, 4], [5, 7, 6]],
1982         [[1, 2, 4, 3], [5, 7, 6]],
1983         [[1, 3, 2, 4], [5, 7, 6]],
1984         [[1, 3, 4, 2], [5, 7, 6]],
1985         [[1, 4, 2, 3], [5, 7, 6]],
1986         [[1, 4, 3, 2], [5, 7, 6]]]
1987
1988    Note that repeated elements are not considered equal:
1989        sage: cyclic_permutations_of_partition([[1,2,3],[4,4,4]])
1990        [[[1, 2, 3], [4, 4, 4]],
1991         [[1, 3, 2], [4, 4, 4]],
1992         [[1, 2, 3], [4, 4, 4]],
1993         [[1, 3, 2], [4, 4, 4]]]
1994
1995    """
1996    return list(cyclic_permutations_of_partition_iterator(partition))
1997
1998def cyclic_permutations_of_partition_iterator(partition):
1999    """
2000    Iterates over all combinations of cyclic permutations of each cell of the
2001    partition.
2002
2003    AUTHOR: Robert L. Miller
2004
2005    EXAMPLES:
2006        sage: from sage.combinat.combinat import cyclic_permutations_of_partition
2007        sage: cyclic_permutations_of_partition([[1,2,3,4],[5,6,7]])
2008        [[[1, 2, 3, 4], [5, 6, 7]],
2009         [[1, 2, 4, 3], [5, 6, 7]],
2010         [[1, 3, 2, 4], [5, 6, 7]],
2011         [[1, 3, 4, 2], [5, 6, 7]],
2012         [[1, 4, 2, 3], [5, 6, 7]],
2013         [[1, 4, 3, 2], [5, 6, 7]],
2014         [[1, 2, 3, 4], [5, 7, 6]],
2015         [[1, 2, 4, 3], [5, 7, 6]],
2016         [[1, 3, 2, 4], [5, 7, 6]],
2017         [[1, 3, 4, 2], [5, 7, 6]],
2018         [[1, 4, 2, 3], [5, 7, 6]],
2019         [[1, 4, 3, 2], [5, 7, 6]]]
2020
2021    Note that repeated elements are not considered equal:
2022        sage: cyclic_permutations_of_partition([[1,2,3],[4,4,4]])
2023        [[[1, 2, 3], [4, 4, 4]],
2024         [[1, 3, 2], [4, 4, 4]],
2025         [[1, 2, 3], [4, 4, 4]],
2026         [[1, 3, 2], [4, 4, 4]]]
2027
2028    """
2029    if len(partition) == 1:
2030        for i in cyclic_permutations_iterator(partition[0]):
2031            yield [i]
2032    else:
2033        for right in cyclic_permutations_of_partition_iterator(partition[1:]):
2034            for perm in cyclic_permutations_iterator(partition[0]):
2035                yield [perm] + right
2036
2037def ferrers_diagram(pi):
2038    """
2039    Return the Ferrers diagram of pi.
2040
2041    INPUT:
2042        pi -- a partition, given as a list of integers.
2043
2044    EXAMPLES:
2045        sage: print ferrers_diagram([5,5,2,1])
2046        *****
2047        *****
2048        **
2049        *
2050        sage: pi = partitions_list(10)[30] ## [6,1,1,1,1]
2051        sage: print ferrers_diagram(pi)
2052        ******
2053        *
2054        *
2055        *
2056        *
2057        sage: pi = partitions_list(10)[33] ## [6, 3, 1]
2058        sage: print ferrers_diagram(pi)
2059        ******
2060        ***
2061        *
2062    """
2063    return '\n'.join(['*'*p for p in pi])
2064
2065
2066def ordered_partitions(n,k=None):
2067    r"""
2068    An {\it ordered partition of $n$} is an ordered sum
2069    $$2070 n = p_1+p_2 + \cdots + p_k 2071$$
2072    of positive integers and is represented by the list $p = [p_1,p_2,\cdots ,p_k]$.
2073    If $k$ is omitted then all ordered partitions are returned.
2074
2075    \code{ordered_partitions(n,k)} returns the set of all (ordered)
2076    partitions of the positive integer n into sums with k summands.
2077
2078    Do not call \code{ordered_partitions} with an n much larger than
2079    15, since the list will simply become too large.
2080
2081    Wraps GAP's OrderedPartitions.
2082
2083    The number of ordered partitions $T_n$ of $\{ 1, 2, ..., n \}$ has the
2084    generating function is
2085    $2086 \sum_n {T_n \over n!} x^n = {1 \over 2-e^x}. 2087$
2088
2089    EXAMPLES:
2090        sage: ordered_partitions(10,2)
2091        [[1, 9], [2, 8], [3, 7], [4, 6], [5, 5], [6, 4], [7, 3], [8, 2], [9, 1]]
2092
2093        sage: ordered_partitions(4)
2094        [[1, 1, 1, 1], [1, 1, 2], [1, 2, 1], [1, 3], [2, 1, 1], [2, 2], [3, 1], [4]]
2095
2096    REFERENCES:
2097        http://en.wikipedia.org/wiki/Ordered_partition_of_a_set
2098
2099    """
2100    if k is None:
2101        ans=gap.eval("OrderedPartitions(%s)"%(ZZ(n)))
2102    else:
2103        ans=gap.eval("OrderedPartitions(%s,%s)"%(ZZ(n),ZZ(k)))
2104    return eval(ans.replace('\n',''))
2105
2106def number_of_ordered_partitions(n,k=None):
2107    """
2108    Returns the size of ordered_partitions(n,k).
2109    Wraps GAP's NrOrderedPartitions.
2110
2111    It is possible to associate with every partition of the integer n a conjugacy
2112    class of permutations in the symmetric group on n points and vice versa.
2113    Therefore p(n) = NrPartitions(n) is the number of conjugacy classes of the
2114    symmetric group on n points.
2115
2116
2117    EXAMPLES:
2118        sage: number_of_ordered_partitions(10,2)
2119        9
2120        sage: number_of_ordered_partitions(15)
2121        16384
2122    """
2123    if k is None:
2124        ans=gap.eval("NrOrderedPartitions(%s)"%(n))
2125    else:
2126        ans=gap.eval("NrOrderedPartitions(%s,%s)"%(n,k))
2127    return ZZ(ans)
2128
2129def partitions_greatest(n,k):
2130    """
2131    Returns the set of all (unordered) restricted'' partitions of the integer n having
2132    parts less than or equal to the integer k.
2133
2134    Wraps GAP's PartitionsGreatestLE.
2135
2136    EXAMPLES:
2137        sage: partitions_greatest(10,2)
2138        [[1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
2139         [2, 1, 1, 1, 1, 1, 1, 1, 1],
2140         [2, 2, 1, 1, 1, 1, 1, 1],
2141         [2, 2, 2, 1, 1, 1, 1],
2142         [2, 2, 2, 2, 1, 1],
2143         [2, 2, 2, 2, 2]]
2144    """
2145    return eval(gap.eval("PartitionsGreatestLE(%s,%s)"%(ZZ(n),ZZ(k))))
2146
2147def partitions_greatest_eq(n,k):
2148    """
2149    Returns the set of all (unordered) restricted'' partitions of the
2150    integer n having at least one part equal to the integer k.
2151
2152    Wraps GAP's PartitionsGreatestEQ.
2153
2154    EXAMPLES:
2155        sage: partitions_greatest_eq(10,2)
2156        [[2, 1, 1, 1, 1, 1, 1, 1, 1],
2157         [2, 2, 1, 1, 1, 1, 1, 1],
2158         [2, 2, 2, 1, 1, 1, 1],
2159         [2, 2, 2, 2, 1, 1],
2160         [2, 2, 2, 2, 2]]
2161
2162    """
2163    ans = gap.eval("PartitionsGreatestEQ(%s,%s)"%(n,k))
2164    return eval(ans)
2165
2166def partitions_restricted(n,S,k=None):
2167    r"""
2168    A {\it restricted partition} is, like an ordinary partition, an
2169    unordered sum $n = p_1+p_2+\ldots+p_k$ of positive integers and is
2170    represented by the list $p = [p_1,p_2,\ldots,p_k]$, in nonincreasing
2171    order. The difference is that here the $p_i$ must be elements
2172    from the set $S$, while for ordinary partitions they may be
2173    elements from $[1..n]$.
2174
2175    Returns the set of all restricted partitions of the positive integer
2176    n into sums with k summands with the summands of the partition coming
2177    from the set $S$. If k is not given all restricted partitions for all
2178    k are returned.
2179
2180    Wraps GAP's RestrictedPartitions.
2181
2182    EXAMPLES:
2183        sage: partitions_restricted(8,[1,3,5,7])
2184        [[1, 1, 1, 1, 1, 1, 1, 1],
2185         [3, 1, 1, 1, 1, 1],
2186         [3, 3, 1, 1],
2187         [5, 1, 1, 1],
2188         [5, 3],
2189         [7, 1]]
2190        sage: partitions_restricted(8,[1,3,5,7],2)
2191        [[5, 3], [7, 1]]
2192    """
2193    if k is None:
2194        ans=gap.eval("RestrictedPartitions(%s,%s)"%(n,S))
2195    else:
2196        ans=gap.eval("RestrictedPartitions(%s,%s,%s)"%(n,S,k))
2197    return eval(ans)
2198
2199def number_of_partitions_restricted(n,S,k=None):
2200    """
2201    Returns the size of partitions_restricted(n,S,k).
2202    Wraps GAP's NrRestrictedPartitions.
2203
2204    EXAMPLES:
2205        sage: number_of_partitions_restricted(8,[1,3,5,7])
2206        6
2207        sage: number_of_partitions_restricted(8,[1,3,5,7],2)
2208        2
2209
2210    """
2211    if k is None:
2212        ans=gap.eval("NrRestrictedPartitions(%s,%s)"%(ZZ(n),S))
2213    else:
2214        ans=gap.eval("NrRestrictedPartitions(%s,%s,%s)"%(ZZ(n),S,ZZ(k)))
2215    return ZZ(ans)
2216
2217def partitions_tuples(n,k):
2218    """
2219    partition_tuples( n, k ) returns the list of all k-tuples of partitions
2220    which together form a partition of n.
2221
2222    k-tuples of partitions describe the classes and the characters of
2223    wreath products of groups with k conjugacy classes with the symmetric
2224    group $S_n$.
2225
2226    Wraps GAP's PartitionTuples.
2227
2228    EXAMPLES:
2229        sage: partitions_tuples(3,2)
2230        [[[1, 1, 1], []],
2231         [[1, 1], [1]],
2232         [[1], [1, 1]],
2233         [[], [1, 1, 1]],
2234         [[2, 1], []],
2235         [[1], [2]],
2236         [[2], [1]],
2237         [[], [2, 1]],
2238         [[3], []],
2239         [[], [3]]]
2240    """
2241    ans=gap.eval("PartitionTuples(%s,%s)"%(ZZ(n),ZZ(k)))
2242    return eval(ans)
2243
2244def number_of_partitions_tuples(n,k):
2245    r"""
2246    number_of_partition_tuples( n, k ) returns the number of partition_tuples(n,k).
2247
2248    Wraps GAP's NrPartitionTuples.
2249
2250    EXAMPLES:
2251        sage: number_of_partitions_tuples(3,2)
2252        10
2253        sage: number_of_partitions_tuples(8,2)
2254        185
2255
2256    Now we compare that with the result of the following GAP
2257    computation:
2258 \begin{verbatim}
2259        gap> S8:=Group((1,2,3,4,5,6,7,8),(1,2));
2260        Group([ (1,2,3,4,5,6,7,8), (1,2) ])
2261        gap> C2:=Group((1,2));
2262        Group([ (1,2) ])
2263        gap> W:=WreathProduct(C2,S8);
2264        <permutation group of size 10321920 with 10 generators>
2265        gap> Size(W);
2266        10321920     ## = 2^8*Factorial(8), which is good:-)
2267        gap> Size(ConjugacyClasses(W));
2268        185
2269\end{verbatim}
2270    """
2271    ans=gap.eval("NrPartitionTuples(%s,%s)"%(ZZ(n),ZZ(k)))
2272    return ZZ(ans)
2273
2274def partition_power(pi,k):
2275    """
2276    partition_power( pi, k ) returns the partition corresponding to the
2277    $k$-th power of a permutation with cycle structure pi
2278    (thus describes the powermap of symmetric groups).
2279
2280    Wraps GAP's PowerPartition.
2281
2282    EXAMPLES:
2283        sage: partition_power([5,3],1)
2284        [5, 3]
2285        sage: partition_power([5,3],2)
2286        [5, 3]
2287        sage: partition_power([5,3],3)
2288        [5, 1, 1, 1]
2289        sage: partition_power([5,3],4)
2290        [5, 3]
2291
2292     Now let us compare this to the power map on $S_8$:
2293
2294        sage: G = SymmetricGroup(8)
2295        sage: g = G([(1,2,3,4,5),(6,7,8)])
2296        sage: g
2297        (1,2,3,4,5)(6,7,8)
2298        sage: g^2
2299        (1,3,5,2,4)(6,8,7)
2300        sage: g^3
2301        (1,4,2,5,3)
2302        sage: g^4
2303        (1,5,4,3,2)(6,7,8)
2304
2305    """
2306    ans=gap.eval("PowerPartition(%s,%s)"%(pi,ZZ(k)))
2307    return eval(ans)
2308
2309def partition_sign(pi):
2310    r"""
2311    partition_sign( pi ) returns the sign of a permutation with cycle structure
2312    given by the partition pi.
2313
2314    This function corresponds to a homomorphism from the symmetric group
2315    $S_n$ into the cyclic group of order 2, whose kernel is exactly the
2316    alternating group $A_n$. Partitions of sign $1$ are called {\it even partitions}
2317    while partitions of sign $-1$ are called {\it odd}.
2318
2319    Wraps GAP's SignPartition.
2320
2321    EXAMPLES:
2322        sage: partition_sign([5,3])
2323        1
2324        sage: partition_sign([5,2])
2325        -1
2326
2327    {\it Zolotarev's lemma} states that the Legendre symbol
2328    $\left(\frac{a}{p}\right)$ for an integer $a \pmod p$ ($p$ a prime number),
2329    can be computed as sign(p_a), where sign denotes the sign of a permutation
2330    and p_a the permutation of the residue classes $\pmod p$ induced by
2331    modular multiplication by $a$, provided $p$ does not divide $a$.
2332
2333    We verify this in some examples.
2334
2335        sage: F = GF(11)
2336        sage: a = F.multiplicative_generator();a
2337        2
2338        sage: plist = [int(a*F(x)) for x in range(1,11)]; plist
2339        [2, 4, 6, 8, 10, 1, 3, 5, 7, 9]
2340
2341    This corresponds ot the permutation (1, 2, 4, 8, 5, 10, 9, 7, 3, 6)
2342    (acting the set $\{1,2,...,10\}$) and to the partition [10].
2343
2344        sage: p = PermutationGroupElement('(1, 2, 4, 8, 5, 10, 9, 7, 3, 6)')
2345        sage: p.sign()
2346        -1
2347        sage: partition_sign([10])
2348        -1
2349        sage: kronecker_symbol(11,2)
2350        -1
2351
2352    Now replace $2$ by $3$:
2353
2354        sage: plist = [int(F(3*x)) for x in range(1,11)]; plist
2355        [3, 6, 9, 1, 4, 7, 10, 2, 5, 8]
2356        sage: range(1,11)
2357        [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
2358        sage: p = PermutationGroupElement('(3,4,8,7,9)')
2359        sage: p.sign()
2360        1
2361        sage: kronecker_symbol(3,11)
2362        1
2363        sage: partition_sign([5,1,1,1,1,1])
2364        1
2365
2366    In both cases, Zolotarev holds.
2367
2368    REFERENCES:
2369        http://en.wikipedia.org/wiki/Zolotarev's_lemma
2370    """
2371    ans=gap.eval("SignPartition(%s)"%(pi))
2372    return sage_eval(ans)
2373
2374def partition_associated(pi):
2375    """
2376    partition_associated( pi ) returns the associated'' (also called
2377    conjugate'' in the literature) partition of the partition pi which is
2378    obtained by transposing the corresponding Ferrers diagram.
2379
2380    Wraps GAP's AssociatedPartition.
2381
2382    EXAMPLES:
2383        sage: partition_associated([2,2])
2384        [2, 2]
2385        sage: partition_associated([6,3,1])
2386        [3, 2, 2, 1, 1, 1]
2387        sage: print ferrers_diagram([6,3,1])
2388        ******
2389        ***
2390        *
2391        sage: print ferrers_diagram([3,2,2,1,1,1])
2392        ***
2393        **
2394        **
2395        *
2396        *
2397        *
2398    """
2399    ans=gap.eval("AssociatedPartition(%s)"%(pi))
2400    return eval(ans)
Note: See TracBrowser for help on using the repository browser.