Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
sagemath
GitHub Repository: sagemath/sagelib
Path: blob/master/sage/combinat/combinat.py
4045 views
1
r"""
2
Combinatorial Functions
3
4
AUTHORS:
5
6
- David Joyner (2006-07): initial implementation.
7
8
- William Stein (2006-07): editing of docs and code; many
9
optimizations, refinements, and bug fixes in corner cases
10
11
- David Joyner (2006-09): bug fix for combinations, added
12
permutations_iterator, combinations_iterator from Python Cookbook,
13
edited docs.
14
15
- David Joyner (2007-11): changed permutations, added hadamard_matrix
16
17
- Florent Hivert (2009-02): combinatorial class cleanup
18
19
- Fredrik Johansson (2010-07): fast implementation of ``stirling_number2``
20
21
This module implements some combinatorial functions, as listed
22
below. For a more detailed description, see the relevant
23
docstrings.
24
25
Sequences:
26
27
28
- Bell numbers, ``bell_number``
29
30
- Bernoulli numbers, ``bernoulli_number`` (though
31
PARI's bernoulli is better)
32
33
- Catalan numbers, ``catalan_number`` (not to be
34
confused with the Catalan constant)
35
36
- Eulerian/Euler numbers, ``euler_number`` (Maxima)
37
38
- Fibonacci numbers, ``fibonacci`` (PARI) and
39
``fibonacci_number`` (GAP) The PARI version is
40
better.
41
42
- Lucas numbers, ``lucas_number1``,
43
``lucas_number2``.
44
45
- Stirling numbers, ``stirling_number1``,
46
``stirling_number2``.
47
48
49
Set-theoretic constructions:
50
51
52
- Combinations of a multiset, ``combinations``,
53
``combinations_iterator``, and
54
``number_of_combinations``. These are unordered
55
selections without repetition of k objects from a multiset S.
56
57
- Arrangements of a multiset, ``arrangements`` and
58
``number_of_arrangements`` These are ordered
59
selections without repetition of k objects from a multiset S.
60
61
- Derangements of a multiset, ``derangements`` and
62
``number_of_derangements``.
63
64
- Tuples of a multiset, ``tuples`` and
65
``number_of_tuples``. An ordered tuple of length k of
66
set S is a ordered selection with repetitions of S and is
67
represented by a sorted list of length k containing elements from
68
S.
69
70
- Unordered tuples of a set, ``unordered_tuple`` and
71
``number_of_unordered_tuples``. An unordered tuple
72
of length k of set S is an unordered selection with repetitions of S
73
and is represented by a sorted list of length k containing elements
74
from S.
75
76
- Permutations of a multiset, ``permutations``,
77
``permutations_iterator``,
78
``number_of_permutations``. A permutation is a list
79
that contains exactly the same elements but possibly in different
80
order.
81
82
83
Partitions:
84
85
86
- Partitions of a set, ``partitions_set``,
87
``number_of_partitions_set``. An unordered partition
88
of set S is a set of pairwise disjoint nonempty sets with union S
89
and is represented by a sorted list of such sets.
90
91
- Partitions of an integer, ``partitions_list``,
92
``number_of_partitions_list``. An unordered
93
partition of n is an unordered sum
94
`n = p_1+p_2 +\ldots+ p_k` of positive integers and is
95
represented by the list `p = [p_1,p_2,\ldots,p_k]`, in
96
nonincreasing order, i.e., `p1\geq p_2 ...\geq p_k`.
97
98
- Ordered partitions of an integer,
99
``ordered_partitions``,
100
``number_of_ordered_partitions``. An ordered
101
partition of n is an ordered sum
102
`n = p_1+p_2 +\ldots+ p_k` of positive integers and is
103
represented by the list `p = [p_1,p_2,\ldots,p_k]`, in
104
nonincreasing order, i.e., `p1\geq p_2 ...\geq p_k`.
105
106
- Restricted partitions of an integer,
107
``partitions_restricted``,
108
``number_of_partitions_restricted``. An unordered
109
restricted partition of n is an unordered sum
110
`n = p_1+p_2 +\ldots+ p_k` of positive integers
111
`p_i` belonging to a given set `S`, and is
112
represented by the list `p = [p_1,p_2,\ldots,p_k]`, in
113
nonincreasing order, i.e., `p1\geq p_2 ...\geq p_k`.
114
115
- ``partitions_greatest`` implements a special type
116
of restricted partition.
117
118
- ``partitions_greatest_eq`` is another type of
119
restricted partition.
120
121
- Tuples of partitions, ``partition_tuples``,
122
``number_of_partition_tuples``. A `k`-tuple
123
of partitions is represented by a list of all `k`-tuples of
124
partitions which together form a partition of `n`.
125
126
- Powers of a partition, ``partition_power(pi, k)``.
127
The power of a partition corresponds to the `k`-th power of
128
a permutation with cycle structure `\pi`.
129
130
- Sign of a partition, ``partition_sign( pi )`` This
131
means the sign of a permutation with cycle structure given by the
132
partition pi.
133
134
- Associated partition, ``partition_associated( pi
135
)`` The "associated" (also called "conjugate" in the
136
literature) partition of the partition pi which is obtained by
137
transposing the corresponding Ferrers diagram.
138
139
- Ferrers diagram, ``ferrers_diagram``. Analogous to
140
the Young diagram of an irreducible representation of
141
`S_n`.
142
143
144
Related functions:
145
146
147
- Bernoulli polynomials, ``bernoulli_polynomial``
148
149
150
Implemented in other modules (listed for completeness):
151
152
The ``sage.rings.arith`` module contains the following
153
combinatorial functions:
154
155
156
- binomial the binomial coefficient (wrapped from PARI)
157
158
- factorial (wrapped from PARI)
159
160
- partition (from the Python Cookbook) Generator of the list of
161
all the partitions of the integer `n`.
162
163
- ``number_of_partitions`` (wrapped from PARI) the
164
*number* of partitions:
165
166
- ``falling_factorial`` Definition: for integer
167
`a \ge 0` we have `x(x-1) \cdots (x-a+1)`. In all
168
other cases we use the GAMMA-function:
169
`\frac {\Gamma(x+1)} {\Gamma(x-a+1)}`.
170
171
- ``rising_factorial`` Definition: for integer
172
`a \ge 0` we have `x(x+1) \cdots (x+a-1)`. In all
173
other cases we use the GAMMA-function:
174
`\frac {\Gamma(x+a)} {\Gamma(x)}`.
175
176
- gaussian_binomial the gaussian binomial
177
178
.. math::
179
180
\binom{n}{k}_q = \frac{(1-q^m)(1-q^{m-1})\cdots (1-q^{m-r+1})} {(1-q)(1-q^2)\cdots (1-q^r)}.
181
182
183
184
185
The ``sage.groups.perm_gps.permgroup_elements``
186
contains the following combinatorial functions:
187
188
189
- matrix method of PermutationGroupElement yielding the
190
permutation matrix of the group element.
191
192
.. TODO::
193
194
GUAVA commands:
195
* MOLS returns a list of n Mutually Orthogonal Latin Squares (MOLS).
196
* VandermondeMat
197
* GrayMat returns a list of all different vectors of length n over
198
the field F, using Gray ordering.
199
Not in GAP:
200
* Rencontres numbers
201
http://en.wikipedia.org/wiki/Rencontres_number
202
203
REFERENCES:
204
205
- http://en.wikipedia.org/wiki/Twelvefold_way (general reference)
206
"""
207
208
#*****************************************************************************
209
# Copyright (C) 2006 David Joyner <[email protected]>,
210
# 2007 Mike Hansen <[email protected]>,
211
# 2006 William Stein <[email protected]>
212
#
213
# Distributed under the terms of the GNU General Public License (GPL)
214
#
215
# This code is distributed in the hope that it will be useful,
216
# but WITHOUT ANY WARRANTY; without even the implied warranty of
217
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
218
# General Public License for more details.
219
#
220
# The full text of the GPL is available at:
221
#
222
# http://www.gnu.org/licenses/
223
#*****************************************************************************
224
from sage.interfaces.all import gap, maxima
225
from sage.rings.all import QQ, ZZ, Integer
226
from sage.rings.arith import bernoulli, binomial
227
from sage.rings.polynomial.polynomial_element import Polynomial
228
from sage.misc.sage_eval import sage_eval
229
from sage.libs.all import pari
230
from sage.misc.prandom import randint
231
from sage.misc.misc import prod
232
from sage.structure.sage_object import SageObject
233
from sage.structure.parent import Parent
234
from sage.misc.lazy_attribute import lazy_attribute
235
from sage.misc.misc import deprecation
236
from combinat_cython import _stirling_number2
237
######### combinatorial sequences
238
239
def bell_number(n):
240
r"""
241
Returns the n-th Bell number (the number of ways to partition a set
242
of n elements into pairwise disjoint nonempty subsets).
243
244
If `n \leq 0`, returns `1`.
245
246
Wraps GAP's Bell.
247
248
EXAMPLES::
249
250
sage: bell_number(10)
251
115975
252
sage: bell_number(2)
253
2
254
sage: bell_number(-10)
255
1
256
sage: bell_number(1)
257
1
258
sage: bell_number(1/3)
259
Traceback (most recent call last):
260
...
261
TypeError: no conversion of this rational to integer
262
"""
263
ans=gap.eval("Bell(%s)"%ZZ(n))
264
return ZZ(ans)
265
266
## def bernoulli_number(n):
267
## r"""
268
## Returns the n-th Bernoulli number $B_n$; $B_n/n!$ is the
269
## coefficient of $x^n$ in the power series of $x/(e^x-1)$.
270
## Wraps GAP's Bernoulli.
271
## EXAMPLES:
272
## sage: bernoulli_number(50)
273
## 495057205241079648212477525/66
274
## """
275
## ans=gap.eval("Bernoulli(%s)"%n)
276
## return QQ(ans) ## use QQ, not eval, here
277
278
def catalan_number(n):
279
r"""
280
Returns the n-th Catalan number
281
282
Catalan numbers: The `n`-th Catalan number is given
283
directly in terms of binomial coefficients by
284
285
.. math::
286
287
C_n = \frac{1}{n+1}{2n\choose n} = \frac{(2n)!}{(n+1)!\,n!} \qquad\mbox{ for }\quad n\ge 0.
288
289
290
291
Consider the set `S = \{ 1, ..., n \}`. A noncrossing
292
partition of `S` is a partition in which no two blocks
293
"cross" each other, i.e., if a and b belong to one block and x and
294
y to another, they are not arranged in the order axby.
295
`C_n` is the number of noncrossing partitions of the set
296
`S`. There are many other interpretations (see
297
REFERENCES).
298
299
When `n=-1`, this function raises a ZeroDivisionError; for
300
other `n<0` it returns `0`.
301
302
INPUT:
303
304
- ``n`` - integer
305
306
OUTPUT: integer
307
308
309
310
EXAMPLES::
311
312
sage: [catalan_number(i) for i in range(7)]
313
[1, 1, 2, 5, 14, 42, 132]
314
sage: taylor((-1/2)*sqrt(1 - 4*x^2), x, 0, 15)
315
132*x^14 + 42*x^12 + 14*x^10 + 5*x^8 + 2*x^6 + x^4 + x^2 - 1/2
316
sage: [catalan_number(i) for i in range(-7,7) if i != -1]
317
[0, 0, 0, 0, 0, 0, 1, 1, 2, 5, 14, 42, 132]
318
sage: catalan_number(-1)
319
Traceback (most recent call last):
320
...
321
ZeroDivisionError
322
sage: [catalan_number(n).mod(2) for n in range(16)]
323
[1, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1]
324
325
REFERENCES:
326
327
- http://en.wikipedia.org/wiki/Catalan_number
328
329
- http://www-history.mcs.st-andrews.ac.uk/~history/Miscellaneous/CatalanNumbers/catalan.html
330
"""
331
from sage.rings.arith import binomial
332
n = ZZ(n)
333
return binomial(2*n,n).divide_knowing_divisible_by(n+1)
334
335
def euler_number(n):
336
"""
337
Returns the n-th Euler number
338
339
IMPLEMENTATION: Wraps Maxima's euler.
340
341
EXAMPLES::
342
343
sage: [euler_number(i) for i in range(10)]
344
[1, 0, -1, 0, 5, 0, -61, 0, 1385, 0]
345
sage: maxima.eval("taylor (2/(exp(x)+exp(-x)), x, 0, 10)")
346
'1-x^2/2+5*x^4/24-61*x^6/720+277*x^8/8064-50521*x^10/3628800'
347
sage: [euler_number(i)/factorial(i) for i in range(11)]
348
[1, 0, -1/2, 0, 5/24, 0, -61/720, 0, 277/8064, 0, -50521/3628800]
349
sage: euler_number(-1)
350
Traceback (most recent call last):
351
...
352
ValueError: n (=-1) must be a nonnegative integer
353
354
REFERENCES:
355
356
- http://en.wikipedia.org/wiki/Euler_number
357
"""
358
n = ZZ(n)
359
if n < 0:
360
raise ValueError, "n (=%s) must be a nonnegative integer"%n
361
return ZZ(maxima.eval("euler(%s)"%n))
362
363
def fibonacci(n, algorithm="pari"):
364
"""
365
Returns the n-th Fibonacci number. The Fibonacci sequence
366
`F_n` is defined by the initial conditions
367
`F_1=F_2=1` and the recurrence relation
368
`F_{n+2} = F_{n+1} + F_n`. For negative `n` we
369
define `F_n = (-1)^{n+1}F_{-n}`, which is consistent with
370
the recurrence relation.
371
372
INPUT:
373
374
375
- ``algorithm`` - string:
376
377
- ``"pari"`` - (default) - use the PARI C library's
378
fibo function.
379
380
- ``"gap"`` - use GAP's Fibonacci function
381
382
383
.. note::
384
385
PARI is tens to hundreds of times faster than GAP here;
386
moreover, PARI works for every large input whereas GAP doesn't.
387
388
EXAMPLES::
389
390
sage: fibonacci(10)
391
55
392
sage: fibonacci(10, algorithm='gap')
393
55
394
395
::
396
397
sage: fibonacci(-100)
398
-354224848179261915075
399
sage: fibonacci(100)
400
354224848179261915075
401
402
::
403
404
sage: fibonacci(0)
405
0
406
sage: fibonacci(1/2)
407
Traceback (most recent call last):
408
...
409
TypeError: no conversion of this rational to integer
410
"""
411
n = ZZ(n)
412
if algorithm == 'pari':
413
return ZZ(pari(n).fibonacci())
414
elif algorithm == 'gap':
415
return ZZ(gap.eval("Fibonacci(%s)"%n))
416
else:
417
raise ValueError, "no algorithm %s"%algorithm
418
419
def lucas_number1(n,P,Q):
420
"""
421
Returns the n-th Lucas number "of the first kind" (this is not
422
standard terminology). The Lucas sequence `L^{(1)}_n` is
423
defined by the initial conditions `L^{(1)}_1=0`,
424
`L^{(1)}_2=1` and the recurrence relation
425
`L^{(1)}_{n+2} = P*L^{(1)}_{n+1} - Q*L^{(1)}_n`.
426
427
Wraps GAP's Lucas(...)[1].
428
429
P=1, Q=-1 gives the Fibonacci sequence.
430
431
INPUT:
432
433
434
- ``n`` - integer
435
436
- ``P, Q`` - integer or rational numbers
437
438
439
OUTPUT: integer or rational number
440
441
EXAMPLES::
442
443
sage: lucas_number1(5,1,-1)
444
5
445
sage: lucas_number1(6,1,-1)
446
8
447
sage: lucas_number1(7,1,-1)
448
13
449
sage: lucas_number1(7,1,-2)
450
43
451
452
::
453
454
sage: lucas_number1(5,2,3/5)
455
229/25
456
sage: lucas_number1(5,2,1.5)
457
Traceback (most recent call last):
458
...
459
TypeError: no canonical coercion from Real Field with 53 bits of precision to Rational Field
460
461
There was a conjecture that the sequence `L_n` defined by
462
`L_{n+2} = L_{n+1} + L_n`, `L_1=1`,
463
`L_2=3`, has the property that `n` prime implies
464
that `L_n` is prime.
465
466
::
467
468
sage: lucas = lambda n : Integer((5/2)*lucas_number1(n,1,-1)+(1/2)*lucas_number2(n,1,-1))
469
sage: [[lucas(n),is_prime(lucas(n)),n+1,is_prime(n+1)] for n in range(15)]
470
[[1, False, 1, False],
471
[3, True, 2, True],
472
[4, False, 3, True],
473
[7, True, 4, False],
474
[11, True, 5, True],
475
[18, False, 6, False],
476
[29, True, 7, True],
477
[47, True, 8, False],
478
[76, False, 9, False],
479
[123, False, 10, False],
480
[199, True, 11, True],
481
[322, False, 12, False],
482
[521, True, 13, True],
483
[843, False, 14, False],
484
[1364, False, 15, False]]
485
486
Can you use Sage to find a counterexample to the conjecture?
487
"""
488
ans=gap.eval("Lucas(%s,%s,%s)[1]"%(QQ._coerce_(P),QQ._coerce_(Q),ZZ(n)))
489
return sage_eval(ans)
490
491
def lucas_number2(n,P,Q):
492
r"""
493
Returns the n-th Lucas number "of the second kind" (this is not
494
standard terminology). The Lucas sequence `L^{(2)}_n` is
495
defined by the initial conditions `L^{(2)}_1=2`,
496
`L^{(2)}_2=P` and the recurrence relation
497
`L^{(2)}_{n+2} = P*L^{(2)}_{n+1} - Q*L^{(2)}_n`.
498
499
Wraps GAP's Lucas(...)[2].
500
501
INPUT:
502
503
504
- ``n`` - integer
505
506
- ``P, Q`` - integer or rational numbers
507
508
509
OUTPUT: integer or rational number
510
511
EXAMPLES::
512
513
sage: [lucas_number2(i,1,-1) for i in range(10)]
514
[2, 1, 3, 4, 7, 11, 18, 29, 47, 76]
515
sage: [fibonacci(i-1)+fibonacci(i+1) for i in range(10)]
516
[2, 1, 3, 4, 7, 11, 18, 29, 47, 76]
517
518
::
519
520
sage: n = lucas_number2(5,2,3); n
521
2
522
sage: type(n)
523
<type 'sage.rings.integer.Integer'>
524
sage: n = lucas_number2(5,2,-3/9); n
525
418/9
526
sage: type(n)
527
<type 'sage.rings.rational.Rational'>
528
529
The case P=1, Q=-1 is the Lucas sequence in Brualdi's Introductory
530
Combinatorics, 4th ed., Prentice-Hall, 2004::
531
532
sage: [lucas_number2(n,1,-1) for n in range(10)]
533
[2, 1, 3, 4, 7, 11, 18, 29, 47, 76]
534
"""
535
ans=gap.eval("Lucas(%s,%s,%s)[2]"%(QQ._coerce_(P),QQ._coerce_(Q),ZZ(n)))
536
return sage_eval(ans)
537
538
def stirling_number1(n,k):
539
"""
540
Returns the n-th Stilling number `S_1(n,k)` of the first
541
kind (the number of permutations of n points with k cycles). Wraps
542
GAP's Stirling1.
543
544
EXAMPLES::
545
546
sage: stirling_number1(3,2)
547
3
548
sage: stirling_number1(5,2)
549
50
550
sage: 9*stirling_number1(9,5)+stirling_number1(9,4)
551
269325
552
sage: stirling_number1(10,5)
553
269325
554
555
Indeed, `S_1(n,k) = S_1(n-1,k-1) + (n-1)S_1(n-1,k)`.
556
"""
557
return ZZ(gap.eval("Stirling1(%s,%s)"%(ZZ(n),ZZ(k))))
558
559
def stirling_number2(n, k, algorithm=None):
560
"""
561
Returns the n-th Stirling number `S_2(n,k)` of the second
562
kind (the number of ways to partition a set of n elements into k
563
pairwise disjoint nonempty subsets). (The n-th Bell number is the
564
sum of the `S_2(n,k)`'s, `k=0,...,n`.)
565
566
INPUT:
567
568
* ``n`` - nonnegative machine-size integer
569
* ``k`` - nonnegative machine-size integer
570
* ``algorithm``:
571
572
* None (default) - use native implementation
573
* ``"maxima"`` - use Maxima's stirling2 function
574
* ``"gap"`` - use GAP's Stirling2 function
575
576
EXAMPLES:
577
578
Print a table of the first several Stirling numbers of the second kind::
579
580
sage: for n in range(10):
581
... for k in range(10):
582
... print str(stirling_number2(n,k)).rjust(k and 6),
583
... print
584
...
585
1 0 0 0 0 0 0 0 0 0
586
0 1 0 0 0 0 0 0 0 0
587
0 1 1 0 0 0 0 0 0 0
588
0 1 3 1 0 0 0 0 0 0
589
0 1 7 6 1 0 0 0 0 0
590
0 1 15 25 10 1 0 0 0 0
591
0 1 31 90 65 15 1 0 0 0
592
0 1 63 301 350 140 21 1 0 0
593
0 1 127 966 1701 1050 266 28 1 0
594
0 1 255 3025 7770 6951 2646 462 36 1
595
596
Stirling numbers satisfy `S_2(n,k) = S_2(n-1,k-1) + kS_2(n-1,k)`::
597
598
sage: 5*stirling_number2(9,5) + stirling_number2(9,4)
599
42525
600
sage: stirling_number2(10,5)
601
42525
602
603
TESTS::
604
605
sage: stirling_number2(500,501)
606
0
607
sage: stirling_number2(500,500)
608
1
609
sage: stirling_number2(500,499)
610
124750
611
sage: stirling_number2(500,498)
612
7739801875
613
sage: stirling_number2(500,497)
614
318420320812125
615
sage: stirling_number2(500,0)
616
0
617
sage: stirling_number2(500,1)
618
1
619
sage: stirling_number2(500,2)
620
1636695303948070935006594848413799576108321023021532394741645684048066898202337277441635046162952078575443342063780035504608628272942696526664263794687
621
sage: stirling_number2(500,3)
622
6060048632644989473730877846590553186337230837666937173391005972096766698597315914033083073801260849147094943827552228825899880265145822824770663507076289563105426204030498939974727520682393424986701281896187487826395121635163301632473646
623
sage: stirling_number2(500,30)
624
13707767141249454929449108424328432845001327479099713037876832759323918134840537229737624018908470350134593241314462032607787062188356702932169472820344473069479621239187226765307960899083230982112046605340713218483809366970996051181537181362810003701997334445181840924364501502386001705718466534614548056445414149016614254231944272872440803657763210998284198037504154374028831561296154209804833852506425742041757849726214683321363035774104866182331315066421119788248419742922490386531970053376982090046434022248364782970506521655684518998083846899028416459701847828711541840099891244700173707021989771147674432503879702222276268661726508226951587152781439224383339847027542755222936463527771486827849728880
625
sage: stirling_number2(500,31)
626
5832088795102666690960147007601603328246123996896731854823915012140005028360632199516298102446004084519955789799364757997824296415814582277055514048635928623579397278336292312275467402957402880590492241647229295113001728653772550743446401631832152281610081188041624848850056657889275564834450136561842528589000245319433225808712628826136700651842562516991245851618481622296716433577650218003181535097954294609857923077238362717189185577756446945178490324413383417876364657995818830270448350765700419876347023578011403646501685001538551891100379932684279287699677429566813471166558163301352211170677774072447414719380996777162087158124939742564291760392354506347716119002497998082844612434332155632097581510486912
627
sage: n = stirling_number2(20,11)
628
sage: n
629
1900842429486
630
sage: type(n)
631
<type 'sage.rings.integer.Integer'>
632
sage: n = stirling_number2(20,11,algorithm='gap')
633
sage: n
634
1900842429486
635
sage: type(n)
636
<type 'sage.rings.integer.Integer'>
637
sage: n = stirling_number2(20,11,algorithm='maxima')
638
sage: n
639
1900842429486
640
sage: type(n)
641
<type 'sage.rings.integer.Integer'>
642
643
Sage's implementation splitting the computation of the Stirling
644
numbers of the second kind in two cases according to `n`, let us
645
check the result it gives agree with both maxima and gap.
646
647
For `n<200`::
648
649
sage: for n in Subsets(range(100,200), 5).random_element():
650
... for k in Subsets(range(n), 5).random_element():
651
... s_sage = stirling_number2(n,k)
652
... s_maxima = stirling_number2(n,k, algorithm = "maxima")
653
... s_gap = stirling_number2(n,k, algorithm = "gap")
654
... if not (s_sage == s_maxima and s_sage == s_gap):
655
... print "Error with n<200"
656
657
For `n\geq 200`::
658
659
sage: for n in Subsets(range(200,300), 5).random_element():
660
... for k in Subsets(range(n), 5).random_element():
661
... s_sage = stirling_number2(n,k)
662
... s_maxima = stirling_number2(n,k, algorithm = "maxima")
663
... s_gap = stirling_number2(n,k, algorithm = "gap")
664
... if not (s_sage == s_maxima and s_sage == s_gap):
665
... print "Error with n<200"
666
667
668
TESTS:
669
670
Checking an exception is raised whenever a wrong value is given
671
for ``algorithm``::
672
673
sage: s_sage = stirling_number2(50,3, algorithm = "CloudReading")
674
Traceback (most recent call last):
675
...
676
ValueError: unknown algorithm: CloudReading
677
"""
678
if algorithm is None:
679
return _stirling_number2(n, k)
680
elif algorithm == 'gap':
681
return ZZ(gap.eval("Stirling2(%s,%s)"%(ZZ(n),ZZ(k))))
682
elif algorithm == 'maxima':
683
return ZZ(maxima.eval("stirling2(%s,%s)"%(ZZ(n),ZZ(k))))
684
else:
685
raise ValueError("unknown algorithm: %s" % algorithm)
686
687
class CombinatorialObject(SageObject):
688
def __init__(self, l):
689
"""
690
CombinatorialObject provides a thin wrapper around a list. The main
691
differences are that __setitem__ is disabled so that
692
CombinatorialObjects are shallowly immutable, and the intention is
693
that they are semantically immutable.
694
695
Because of this, CombinatorialObjects provide a __hash__
696
function which computes the hash of the string representation of a
697
list and the hash of its parent's class. Thus, each
698
CombinatorialObject should have a unique string representation.
699
700
INPUT:
701
702
- ``l`` - a list or any object that can be convert to a list by
703
``list``
704
705
EXAMPLES::
706
707
sage: c = CombinatorialObject([1,2,3])
708
sage: c == loads(dumps(c))
709
True
710
sage: c._list
711
[1, 2, 3]
712
sage: c._hash is None
713
True
714
"""
715
if isinstance(l, list):
716
self._list = l
717
else:
718
self._list = list(l)
719
self._hash = None
720
721
def __str__(self):
722
"""
723
EXAMPLES::
724
725
sage: c = CombinatorialObject([1,2,3])
726
sage: str(c)
727
'[1, 2, 3]'
728
"""
729
return str(self._list)
730
731
def _repr_(self):
732
"""
733
EXAMPLES::
734
735
sage: c = CombinatorialObject([1,2,3])
736
sage: c.__repr__()
737
'[1, 2, 3]'
738
"""
739
return self._list.__repr__()
740
741
def __eq__(self, other):
742
"""
743
EXAMPLES::
744
745
sage: c = CombinatorialObject([1,2,3])
746
sage: d = CombinatorialObject([2,3,4])
747
sage: c == [1,2,3]
748
True
749
sage: c == [2,3,4]
750
False
751
sage: c == d
752
False
753
"""
754
755
if isinstance(other, CombinatorialObject):
756
return self._list.__eq__(other._list)
757
else:
758
return self._list.__eq__(other)
759
760
def __lt__(self, other):
761
"""
762
EXAMPLES::
763
764
sage: c = CombinatorialObject([1,2,3])
765
sage: d = CombinatorialObject([2,3,4])
766
sage: c < d
767
True
768
sage: c < [2,3,4]
769
True
770
"""
771
772
if isinstance(other, CombinatorialObject):
773
return self._list.__lt__(other._list)
774
else:
775
return self._list.__lt__(other)
776
777
def __le__(self, other):
778
"""
779
EXAMPLES::
780
781
sage: c = CombinatorialObject([1,2,3])
782
sage: d = CombinatorialObject([2,3,4])
783
sage: c <= c
784
True
785
sage: c <= d
786
True
787
sage: c <= [1,2,3]
788
True
789
"""
790
if isinstance(other, CombinatorialObject):
791
return self._list.__le__(other._list)
792
else:
793
return self._list.__le__(other)
794
795
def __gt__(self, other):
796
"""
797
EXAMPLES::
798
799
sage: c = CombinatorialObject([1,2,3])
800
sage: d = CombinatorialObject([2,3,4])
801
sage: c > c
802
False
803
sage: c > d
804
False
805
sage: c > [1,2,3]
806
False
807
"""
808
if isinstance(other, CombinatorialObject):
809
return self._list.__gt__(other._list)
810
else:
811
return self._list.__gt__(other)
812
813
def __ge__(self, other):
814
"""
815
EXAMPLES::
816
817
sage: c = CombinatorialObject([1,2,3])
818
sage: d = CombinatorialObject([2,3,4])
819
sage: c >= c
820
True
821
sage: c >= d
822
False
823
sage: c >= [1,2,3]
824
True
825
"""
826
if isinstance(other, CombinatorialObject):
827
return self._list.__ge__(other._list)
828
else:
829
return self._list.__ge__(other)
830
831
def __ne__(self, other):
832
"""
833
EXAMPLES::
834
835
sage: c = CombinatorialObject([1,2,3])
836
sage: d = CombinatorialObject([2,3,4])
837
sage: c != c
838
False
839
sage: c != d
840
True
841
sage: c != [1,2,3]
842
False
843
"""
844
if isinstance(other, CombinatorialObject):
845
return self._list.__ne__(other._list)
846
else:
847
return self._list.__ne__(other)
848
849
def __add__(self, other):
850
"""
851
EXAMPLES::
852
853
sage: c = CombinatorialObject([1,2,3])
854
sage: c + [4]
855
[1, 2, 3, 4]
856
sage: type(_)
857
<type 'list'>
858
"""
859
return self._list + other
860
861
def __hash__(self):
862
"""
863
Computes the hash of self by computing the hash of the string
864
representation of self._list. The hash is cached and stored in
865
self._hash.
866
867
EXAMPLES::
868
869
sage: c = CombinatorialObject([1,2,3])
870
sage: c._hash is None
871
True
872
sage: hash(c) #random
873
1335416675971793195
874
sage: c._hash #random
875
1335416675971793195
876
"""
877
if self._hash is None:
878
self._hash = str(self._list).__hash__()
879
return self._hash
880
881
def __len__(self):
882
"""
883
EXAMPLES::
884
885
sage: c = CombinatorialObject([1,2,3])
886
sage: len(c)
887
3
888
sage: c.__len__()
889
3
890
"""
891
return self._list.__len__()
892
893
def __getitem__(self, key):
894
"""
895
EXAMPLES::
896
897
sage: c = CombinatorialObject([1,2,3])
898
sage: c[0]
899
1
900
sage: c[1:]
901
[2, 3]
902
sage: type(_)
903
<type 'list'>
904
"""
905
return self._list.__getitem__(key)
906
907
def __iter__(self):
908
"""
909
EXAMPLES::
910
911
sage: c = CombinatorialObject([1,2,3])
912
sage: list(iter(c))
913
[1, 2, 3]
914
"""
915
return self._list.__iter__()
916
917
def __contains__(self, item):
918
"""
919
EXAMPLES::
920
921
sage: c = CombinatorialObject([1,2,3])
922
sage: 1 in c
923
True
924
sage: 5 in c
925
False
926
"""
927
return self._list.__contains__(item)
928
929
930
def index(self, key):
931
"""
932
EXAMPLES::
933
934
sage: c = CombinatorialObject([1,2,3])
935
sage: c.index(1)
936
0
937
sage: c.index(3)
938
2
939
"""
940
return self._list.index(key)
941
942
943
from sage.misc.classcall_metaclass import ClasscallMetaclass
944
from sage.categories.enumerated_sets import EnumeratedSets
945
class CombinatorialClass(Parent):
946
"""
947
This class is deprecated, and will disappear as soon as all derived
948
classes in Sage's library will have been fixed. Please derive
949
directly from Parent and use the category :class:`EnumeratedSets`,
950
:class:`FiniteEnumeratedSets`, or :class:`InfiniteEnumeratedSets`, as
951
appropriate.
952
953
For examples, see::
954
955
sage: FiniteEnumeratedSets().example()
956
An example of a finite enumerated set: {1,2,3}
957
sage: InfiniteEnumeratedSets().example()
958
An example of an infinite enumerated set: the non negative integers
959
"""
960
__metaclass__ = ClasscallMetaclass
961
962
def __init__(self, category = None, *keys, **opts):
963
"""
964
TESTS::
965
966
sage: C = sage.combinat.combinat.CombinatorialClass()
967
sage: C.category()
968
Category of enumerated sets
969
sage: C.__class__
970
<class 'sage.combinat.combinat.CombinatorialClass_with_category'>
971
sage: isinstance(C, Parent)
972
True
973
sage: C = sage.combinat.combinat.CombinatorialClass(category = FiniteEnumeratedSets())
974
sage: C.category()
975
Category of finite enumerated sets
976
"""
977
Parent.__init__(self, category = EnumeratedSets().or_subcategory(category))
978
979
980
def __len__(self):
981
"""
982
__len__ has been removed ! to get the number of element in a
983
combinatorial class, use .cardinality instead.
984
985
986
TEST::
987
988
sage: len(Partitions(5))
989
Traceback (most recent call last):
990
...
991
AttributeError: __len__ has been removed; use .cardinality() instead
992
"""
993
raise AttributeError, "__len__ has been removed; use .cardinality() instead"
994
995
def count(self):
996
"""
997
Deprecated ! Please use ``.cardinality`` instead.
998
999
1000
TEST::
1001
1002
sage: class C(CombinatorialClass):
1003
... def __iter__(self):
1004
... return iter([1,2,3])
1005
...
1006
sage: C().count() #indirect doctest
1007
doctest:1: DeprecationWarning: The usage of count for combinatorial classes is deprecated. Please use cardinality
1008
3
1009
"""
1010
deprecation("The usage of count for combinatorial classes is deprecated. Please use cardinality")
1011
return self.cardinality()
1012
1013
def is_finite(self):
1014
"""
1015
Returns whether self is finite or not.
1016
1017
EXAMPLES::
1018
1019
sage: Partitions(5).is_finite()
1020
True
1021
sage: Permutations().is_finite()
1022
False
1023
"""
1024
from sage.rings.all import infinity
1025
return self.cardinality() != infinity
1026
1027
def __getitem__(self, i):
1028
"""
1029
Returns the combinatorial object of rank i.
1030
1031
EXAMPLES::
1032
1033
sage: p5 = Partitions(5)
1034
sage: p5[0]
1035
[5]
1036
sage: p5[6]
1037
[1, 1, 1, 1, 1]
1038
sage: p5[7]
1039
Traceback (most recent call last):
1040
...
1041
ValueError: the value must be between 0 and 6 inclusive
1042
"""
1043
return self.unrank(i)
1044
1045
def __str__(self):
1046
"""
1047
Returns a string representation of self.
1048
1049
EXAMPLES::
1050
1051
sage: str(Partitions(5))
1052
'Partitions of the integer 5'
1053
"""
1054
return self.__repr__()
1055
1056
def _repr_(self):
1057
"""
1058
EXAMPLES::
1059
1060
sage: repr(Partitions(5)) # indirect doctest
1061
'Partitions of the integer 5'
1062
"""
1063
if hasattr(self, '_name') and self._name:
1064
return self._name
1065
else:
1066
return "Combinatorial Class -- REDEFINE ME!"
1067
1068
def __contains__(self, x):
1069
"""
1070
Tests whether or not the combinatorial class contains the object x.
1071
This raises a NotImplementedError as a default since _all_
1072
subclasses of CombinatorialClass should override this.
1073
1074
Note that we could replace this with a default implementation that
1075
just iterates through the elements of the combinatorial class and
1076
checks for equality. However, since we use __contains__ for
1077
type checking, this operation should be cheap and should be
1078
implemented manually for each combinatorial class.
1079
1080
EXAMPLES::
1081
1082
sage: C = CombinatorialClass()
1083
sage: x in C
1084
Traceback (most recent call last):
1085
...
1086
NotImplementedError
1087
"""
1088
raise NotImplementedError
1089
1090
def __cmp__(self, x):
1091
"""
1092
Compares two different combinatorial classes. For now, the
1093
comparison is done just on their repr's.
1094
1095
EXAMPLES::
1096
1097
sage: p5 = Partitions(5)
1098
sage: p6 = Partitions(6)
1099
sage: repr(p5) == repr(p6)
1100
False
1101
sage: p5 == p6
1102
False
1103
"""
1104
return cmp(repr(self), repr(x))
1105
1106
def __cardinality_from_iterator(self):
1107
"""
1108
Default implementation of cardinality which just goes through the iterator
1109
of the combinatorial class to count the number of objects.
1110
1111
EXAMPLES::
1112
1113
sage: class C(CombinatorialClass):
1114
... def __iter__(self):
1115
... return iter([1,2,3])
1116
...
1117
sage: C().cardinality() #indirect doctest
1118
3
1119
"""
1120
c = Integer(0)
1121
one = Integer(1)
1122
for _ in self:
1123
c += one
1124
return c
1125
cardinality = __cardinality_from_iterator
1126
1127
# __call__, element_class, and _element_constructor_ are poor
1128
# man's versions of those from Parent. This is for transition,
1129
# until all combinatorial classes are proper parents (in Parent)
1130
# and use coercion, etcc
1131
1132
def __call__(self, x):
1133
"""
1134
Returns x as an element of the combinatorial class's object class.
1135
1136
EXAMPLES::
1137
1138
sage: p5 = Partitions(5)
1139
sage: a = [2,2,1]
1140
sage: type(a)
1141
<type 'list'>
1142
sage: a = p5(a)
1143
sage: type(a)
1144
<class 'sage.combinat.partition.Partition_class'>
1145
sage: p5([2,1])
1146
Traceback (most recent call last):
1147
...
1148
ValueError: [2, 1] not in Partitions of the integer 5
1149
"""
1150
if x in self:
1151
return self._element_constructor_(x)
1152
else:
1153
raise ValueError, "%s not in %s"%(x, self)
1154
1155
Element = CombinatorialObject # mostly for backward compatibility
1156
1157
@lazy_attribute
1158
def element_class(self):
1159
"""
1160
This function is a temporary helper so that a CombinatorialClass
1161
behaves as a parent for creating elements. This will disappear when
1162
combinatorial classes will be turned into actual parents (in the
1163
category EnumeratedSets).
1164
1165
TESTS::
1166
1167
sage: P5 = Partitions(5)
1168
sage: P5.element_class
1169
<class 'sage.combinat.partition.Partition_class'>
1170
"""
1171
# assert not isinstance(self, Parent) # Raises an alert if we override the proper definition from Parent
1172
if hasattr(self, "object_class"):
1173
from sage.misc.misc import deprecation
1174
deprecation("Using object_class for specifying the class of the elements of a combinatorial class is deprecated. Please use Element instead")
1175
return self.Element
1176
1177
def _element_constructor_(self, x):
1178
"""
1179
This function is a temporary helper so that a CombinatorialClass
1180
behaves as a parent for creating elements. This will disappear when
1181
combinatorial classes will be turned into actual parents (in the
1182
category EnumeratedSets).
1183
1184
TESTS::
1185
1186
sage: P5 = Partitions(5)
1187
sage: p = P5([3,2]) # indirect doctest
1188
sage: type(p)
1189
<class 'sage.combinat.partition.Partition_class'>
1190
"""
1191
# assert not isinstance(self, Parent) # Raises an alert if we override the proper definition from Parent
1192
return self.element_class(x)
1193
1194
def __list_from_iterator(self):
1195
"""
1196
The default implementation of list which builds the list from the
1197
iterator.
1198
1199
EXAMPLES::
1200
1201
sage: class C(CombinatorialClass):
1202
... def __iter__(self):
1203
... return iter([1,2,3])
1204
...
1205
sage: C().list() #indirect doctest
1206
[1, 2, 3]
1207
"""
1208
return [x for x in self]
1209
1210
#Set list to the default implementation
1211
list = __list_from_iterator
1212
1213
#Set the default object class to be CombinatorialObject
1214
Element = CombinatorialObject
1215
1216
def __iterator_from_next(self):
1217
"""
1218
An iterator to use when .first() and .next() are provided.
1219
1220
EXAMPLES::
1221
1222
sage: C = CombinatorialClass()
1223
sage: C.first = lambda: 0
1224
sage: C.next = lambda c: c+1
1225
sage: it = iter(C) # indirect doctest
1226
sage: [it.next() for _ in range(4)]
1227
[0, 1, 2, 3]
1228
"""
1229
f = self.first()
1230
yield f
1231
while True:
1232
try:
1233
f = self.next(f)
1234
except (TypeError, ValueError ):
1235
break
1236
1237
if f is None or f is False :
1238
break
1239
else:
1240
yield f
1241
1242
def __iterator_from_previous(self):
1243
"""
1244
An iterator to use when .last() and .previous() are provided. Note
1245
that this requires the combinatorial class to be finite. It is not
1246
recommended to implement combinatorial classes using last and
1247
previous.
1248
1249
EXAMPLES::
1250
1251
sage: C = CombinatorialClass()
1252
sage: C.last = lambda: 4
1253
sage: def prev(c):
1254
... if c <= 1:
1255
... return None
1256
... else:
1257
... return c-1
1258
...
1259
sage: C.previous = prev
1260
sage: it = iter(C) # indirect doctest
1261
sage: [it.next() for _ in range(4)]
1262
[1, 2, 3, 4]
1263
"""
1264
l = self.last()
1265
li = [l]
1266
while True:
1267
try:
1268
l = self.previous(l)
1269
except (TypeError, ValueError):
1270
break
1271
1272
if l == None:
1273
break
1274
else:
1275
li.append(l)
1276
return reversed(li)
1277
1278
def __iterator_from_unrank(self):
1279
"""
1280
An iterator to use when .unrank() is provided.
1281
1282
EXAMPLES::
1283
1284
sage: C = CombinatorialClass()
1285
sage: l = [1,2,3]
1286
sage: C.unrank = lambda c: l[c]
1287
sage: list(C) # indirect doctest
1288
[1, 2, 3]
1289
"""
1290
r = 0
1291
u = self.unrank(r)
1292
yield u
1293
while True:
1294
r += 1
1295
try:
1296
u = self.unrank(r)
1297
except (TypeError, ValueError, IndexError):
1298
break
1299
1300
if u == None:
1301
break
1302
else:
1303
yield u
1304
1305
def __iterator_from_list(self):
1306
"""
1307
An iterator to use when .list() is provided()
1308
1309
EXAMPLES::
1310
1311
sage: C = CombinatorialClass()
1312
sage: C.list = lambda: [1, 2, 3]
1313
sage: list(C) # indirect doctest
1314
[1, 2, 3]
1315
"""
1316
for x in self.list():
1317
yield x
1318
1319
def iterator(self):
1320
"""
1321
Iterator is deprecated for combinatorial classes.
1322
1323
EXAMPLES::
1324
1325
sage: p5 = Partitions(3)
1326
sage: it = p5.iterator()
1327
doctest:1: DeprecationWarning: The usage of iterator for combinatorial classes is deprecated. Please use the class itself
1328
sage: [i for i in it]
1329
[[3], [2, 1], [1, 1, 1]]
1330
sage: [i for i in p5]
1331
[[3], [2, 1], [1, 1, 1]]
1332
"""
1333
from sage.misc.misc import deprecation
1334
deprecation("The usage of iterator for combinatorial classes is deprecated. Please use the class itself")
1335
return self.__iter__()
1336
1337
def __iter__(self):
1338
"""
1339
Allows the combinatorial class to be treated as an iterator. Default
1340
implementation.
1341
1342
EXAMPLES::
1343
1344
sage: p5 = Partitions(5)
1345
sage: [i for i in p5]
1346
[[5], [4, 1], [3, 2], [3, 1, 1], [2, 2, 1], [2, 1, 1, 1], [1, 1, 1, 1, 1]]
1347
sage: C = CombinatorialClass()
1348
sage: iter(C)
1349
Traceback (most recent call last):
1350
...
1351
NotImplementedError: iterator called but not implemented
1352
"""
1353
#Check to see if .first() and .next() are overridden in the subclass
1354
if ( self.first != self.__first_from_iterator and
1355
self.next != self.__next_from_iterator ):
1356
return self.__iterator_from_next()
1357
#Check to see if .last() and .previous() are overridden in the subclass
1358
elif ( self.last != self.__last_from_iterator and
1359
self.previous != self.__previous_from_iterator):
1360
return self.__iterator_from_previous()
1361
#Check to see if .unrank() is overridden in the subclass
1362
elif self.unrank != self.__unrank_from_iterator:
1363
return self.__iterator_from_unrank()
1364
#Finally, check to see if .list() is overridden in the subclass
1365
elif self.list != self.__list_from_iterator:
1366
return self.__iterator_from_list()
1367
else:
1368
raise NotImplementedError, "iterator called but not implemented"
1369
1370
def __unrank_from_iterator(self, r):
1371
"""
1372
Default implementation of unrank which goes through the iterator.
1373
1374
EXAMPLES::
1375
1376
sage: C = CombinatorialClass()
1377
sage: C.list = lambda: [1,2,3]
1378
sage: C.unrank(1) # indirect doctest
1379
2
1380
"""
1381
counter = 0
1382
for u in self:
1383
if counter == r:
1384
return u
1385
counter += 1
1386
raise ValueError, "the value must be between %s and %s inclusive"%(0,counter-1)
1387
1388
#Set the default implementation of unrank
1389
unrank = __unrank_from_iterator
1390
1391
1392
def __random_element_from_unrank(self):
1393
"""
1394
Default implementation of random which uses unrank.
1395
1396
EXAMPLES::
1397
1398
sage: C = CombinatorialClass()
1399
sage: C.list = lambda: [1,2,3]
1400
sage: C.random_element() # indirect doctest
1401
1
1402
"""
1403
c = self.cardinality()
1404
r = randint(0, c-1)
1405
return self.unrank(r)
1406
1407
1408
#Set the default implementation of random
1409
random_element = __random_element_from_unrank
1410
1411
def random(self):
1412
"""
1413
Deprecated. Use self.random_element() instead.
1414
1415
EXAMPLES::
1416
1417
sage: C = CombinatorialClass()
1418
sage: C.random()
1419
Traceback (most recent call last):
1420
...
1421
NotImplementedError: Deprecated: use random_element() instead
1422
"""
1423
raise NotImplementedError, "Deprecated: use random_element() instead"
1424
1425
def __rank_from_iterator(self, obj):
1426
"""
1427
Default implementation of rank which uses iterator.
1428
1429
EXAMPLES::
1430
1431
sage: C = CombinatorialClass()
1432
sage: C.list = lambda: [1,2,3]
1433
sage: C.rank(3) # indirect doctest
1434
2
1435
"""
1436
r = 0
1437
for i in self:
1438
if i == obj:
1439
return r
1440
r += 1
1441
raise ValueError
1442
1443
rank = __rank_from_iterator
1444
1445
def __first_from_iterator(self):
1446
"""
1447
Default implementation for first which uses iterator.
1448
1449
EXAMPLES::
1450
1451
sage: C = CombinatorialClass()
1452
sage: C.list = lambda: [1,2,3]
1453
sage: C.first() # indirect doctest
1454
1
1455
"""
1456
for i in self:
1457
return i
1458
1459
first = __first_from_iterator
1460
1461
def __last_from_iterator(self):
1462
"""
1463
Default implementation for first which uses iterator.
1464
1465
EXAMPLES::
1466
1467
sage: C = CombinatorialClass()
1468
sage: C.list = lambda: [1,2,3]
1469
sage: C.last() # indirect doctest
1470
3
1471
"""
1472
for i in self:
1473
pass
1474
return i
1475
1476
last = __last_from_iterator
1477
1478
def __next_from_iterator(self, obj):
1479
"""
1480
Default implementation for next which uses iterator.
1481
1482
EXAMPLES::
1483
1484
sage: C = CombinatorialClass()
1485
sage: C.list = lambda: [1,2,3]
1486
sage: C.next(2) # indirect doctest
1487
3
1488
"""
1489
found = False
1490
for i in self:
1491
if found:
1492
return i
1493
if i == obj:
1494
found = True
1495
return None
1496
1497
next = __next_from_iterator
1498
1499
def __previous_from_iterator(self, obj):
1500
"""
1501
Default implementation for next which uses iterator.
1502
1503
EXAMPLES::
1504
1505
sage: C = CombinatorialClass()
1506
sage: C.list = lambda: [1,2,3]
1507
sage: C.previous(2) # indirect doctest
1508
1
1509
"""
1510
prev = None
1511
for i in self:
1512
if i == obj:
1513
break
1514
prev = i
1515
return prev
1516
1517
previous = __previous_from_iterator
1518
1519
def filter(self, f, name=None):
1520
"""
1521
Returns the combinatorial subclass of f which consists of the
1522
elements x of self such that f(x) is True.
1523
1524
EXAMPLES::
1525
1526
sage: P = Permutations(3).filter(lambda x: x.avoids([1,2]))
1527
sage: P.list()
1528
[[3, 2, 1]]
1529
"""
1530
return FilteredCombinatorialClass(self, f, name=name)
1531
1532
def union(self, right_cc, name=None):
1533
"""
1534
Returns the combinatorial class representing the union of self and
1535
right_cc.
1536
1537
EXAMPLES::
1538
1539
sage: P = Permutations(2).union(Permutations(1))
1540
sage: P.list()
1541
[[1, 2], [2, 1], [1]]
1542
"""
1543
if not isinstance(right_cc, CombinatorialClass):
1544
raise TypeError, "right_cc must be a CombinatorialClass"
1545
return UnionCombinatorialClass(self, right_cc, name=name)
1546
1547
def map(self, f, name=None):
1548
r"""
1549
Returns the image `\{f(x) | x \in \text{self}\}` of this combinatorial
1550
class by `f`, as a combinatorial class.
1551
1552
`f` is supposed to be injective.
1553
1554
EXAMPLES::
1555
1556
sage: R = Permutations(3).map(attrcall('reduced_word')); R
1557
Image of Standard permutations of 3 by *.reduced_word()
1558
sage: R.cardinality()
1559
6
1560
sage: R.list()
1561
[[], [2], [1], [1, 2], [2, 1], [2, 1, 2]]
1562
sage: [ r for r in R]
1563
[[], [2], [1], [1, 2], [2, 1], [2, 1, 2]]
1564
1565
If the function is not injective, then there may be repeated elements:
1566
sage: P = Partitions(4)
1567
sage: P.list()
1568
[[4], [3, 1], [2, 2], [2, 1, 1], [1, 1, 1, 1]]
1569
sage: P.map(len).list()
1570
[1, 2, 2, 3, 4]
1571
1572
TESTS::
1573
1574
sage: R = Permutations(3).map(attrcall('reduced_word'))
1575
sage: R == loads(dumps(R))
1576
True
1577
"""
1578
return MapCombinatorialClass(self, f, name)
1579
1580
class FilteredCombinatorialClass(CombinatorialClass):
1581
def __init__(self, combinatorial_class, f, name=None):
1582
"""
1583
A filtered combinatorial class F is a subset of another
1584
combinatorial class C specified by a function f that takes in an
1585
element c of C and returns True if and only if c is in F.
1586
1587
TESTS::
1588
1589
sage: Permutations(3).filter(lambda x: x.avoids([1,2]))
1590
Filtered sublass of Standard permutations of 3
1591
"""
1592
self.f = f
1593
self.combinatorial_class = combinatorial_class
1594
self._name = name
1595
1596
def __repr__(self):
1597
"""
1598
EXAMPLES::
1599
1600
sage: P = Permutations(3).filter(lambda x: x.avoids([1,2]))
1601
sage: P.__repr__()
1602
'Filtered sublass of Standard permutations of 3'
1603
sage: P._name = 'Permutations avoiding [1, 2]'
1604
sage: P.__repr__()
1605
'Permutations avoiding [1, 2]'
1606
"""
1607
if self._name:
1608
return self._name
1609
else:
1610
return "Filtered sublass of " + repr(self.combinatorial_class)
1611
1612
def __contains__(self, x):
1613
"""
1614
EXAMPLES::
1615
1616
sage: P = Permutations(3).filter(lambda x: x.avoids([1,2]))
1617
sage: 'cat' in P
1618
False
1619
sage: [4,3,2,1] in P
1620
False
1621
sage: Permutation([1,2,3]) in P
1622
False
1623
sage: Permutation([3,2,1]) in P
1624
True
1625
"""
1626
return x in self.combinatorial_class and self.f(x)
1627
1628
def cardinality(self):
1629
"""
1630
EXAMPLES::
1631
1632
sage: P = Permutations(3).filter(lambda x: x.avoids([1,2]))
1633
sage: P.cardinality()
1634
1
1635
"""
1636
c = 0
1637
for _ in self:
1638
c += 1
1639
return c
1640
1641
def __iter__(self):
1642
"""
1643
EXAMPLES::
1644
1645
sage: P = Permutations(3).filter(lambda x: x.avoids([1,2]))
1646
sage: list(P)
1647
[[3, 2, 1]]
1648
"""
1649
for x in self.combinatorial_class:
1650
if self.f(x):
1651
yield x
1652
1653
class UnionCombinatorialClass(CombinatorialClass):
1654
def __init__(self, left_cc, right_cc, name=None):
1655
"""
1656
A UnionCombinatorialClass is a union of two other combinatorial
1657
classes.
1658
1659
TESTS::
1660
1661
sage: P = Permutations(3).union(Permutations(2))
1662
sage: P == loads(dumps(P))
1663
True
1664
"""
1665
self.left_cc = left_cc
1666
self.right_cc = right_cc
1667
self._name = name
1668
1669
def __repr__(self):
1670
"""
1671
TESTS::
1672
1673
sage: print repr(Permutations(3).union(Permutations(2)))
1674
Union combinatorial class of
1675
Standard permutations of 3
1676
and
1677
Standard permutations of 2
1678
"""
1679
if self._name:
1680
return self._name
1681
else:
1682
return "Union combinatorial class of \n %s\nand\n %s"%(self.left_cc, self.right_cc)
1683
1684
def __contains__(self, x):
1685
"""
1686
EXAMPLES::
1687
1688
sage: P = Permutations(3).union(Permutations(2))
1689
sage: [1,2] in P
1690
True
1691
sage: [3,2,1] in P
1692
True
1693
sage: [1,2,3,4] in P
1694
False
1695
"""
1696
return x in self.left_cc or x in self.right_cc
1697
1698
def cardinality(self):
1699
"""
1700
EXAMPLES::
1701
1702
sage: P = Permutations(3).union(Permutations(2))
1703
sage: P.cardinality()
1704
8
1705
"""
1706
return self.left_cc.cardinality() + self.right_cc.cardinality()
1707
1708
def list(self):
1709
"""
1710
EXAMPLES::
1711
1712
sage: P = Permutations(3).union(Permutations(2))
1713
sage: P.list()
1714
[[1, 2, 3],
1715
[1, 3, 2],
1716
[2, 1, 3],
1717
[2, 3, 1],
1718
[3, 1, 2],
1719
[3, 2, 1],
1720
[1, 2],
1721
[2, 1]]
1722
"""
1723
return self.left_cc.list() + self.right_cc.list()
1724
1725
1726
def __iter__(self):
1727
"""
1728
EXAMPLES::
1729
1730
sage: P = Permutations(3).union(Permutations(2))
1731
sage: list(P)
1732
[[1, 2, 3],
1733
[1, 3, 2],
1734
[2, 1, 3],
1735
[2, 3, 1],
1736
[3, 1, 2],
1737
[3, 2, 1],
1738
[1, 2],
1739
[2, 1]]
1740
"""
1741
for x in self.left_cc:
1742
yield x
1743
for x in self.right_cc:
1744
yield x
1745
1746
def first(self):
1747
"""
1748
EXAMPLES::
1749
1750
sage: P = Permutations(3).union(Permutations(2))
1751
sage: P.first()
1752
[1, 2, 3]
1753
"""
1754
return self.left_cc.first()
1755
1756
def last(self):
1757
"""
1758
EXAMPLES::
1759
1760
sage: P = Permutations(3).union(Permutations(2))
1761
sage: P.last()
1762
[2, 1]
1763
"""
1764
return self.right_cc.last()
1765
1766
def rank(self, x):
1767
"""
1768
EXAMPLES::
1769
1770
sage: P = Permutations(3).union(Permutations(2))
1771
sage: P.rank(Permutation([2,1]))
1772
7
1773
sage: P.rank(Permutation([1,2,3]))
1774
0
1775
"""
1776
try:
1777
return self.left_cc.rank(x)
1778
except (TypeError, ValueError):
1779
return self.left_cc.cardinality() + self.right_cc.rank(x)
1780
1781
def unrank(self, x):
1782
"""
1783
EXAMPLES::
1784
1785
sage: P = Permutations(3).union(Permutations(2))
1786
sage: P.unrank(7)
1787
[2, 1]
1788
sage: P.unrank(0)
1789
[1, 2, 3]
1790
"""
1791
try:
1792
return self.left_cc.unrank(x)
1793
except (TypeError, ValueError):
1794
return self.right_cc.unrank(x - self.left_cc.cardinality())
1795
1796
1797
##############################################################################
1798
class MapCombinatorialClass(CombinatorialClass):
1799
r"""
1800
A MapCombinatorialClass models the image of a combinatorial
1801
class through a function which is assumed to be injective
1802
1803
See CombinatorialClass.map for examples
1804
"""
1805
def __init__(self, cc, f, name=None):
1806
"""
1807
TESTS::
1808
1809
sage: Partitions(3).map(attrcall('conjugate'))
1810
Image of Partitions of the integer 3 by *.conjugate()
1811
"""
1812
self.cc = cc
1813
self.f = f
1814
self._name = name
1815
1816
def __repr__(self):
1817
"""
1818
TESTS::
1819
1820
sage: Partitions(3).map(attrcall('conjugate'))
1821
Image of Partitions of the integer 3 by *.conjugate()
1822
1823
"""
1824
if self._name:
1825
return self._name
1826
else:
1827
return "Image of %s by %s"%(self.cc, self.f)
1828
1829
def cardinality(self):
1830
"""
1831
Returns the cardinality of this combinatorial class
1832
1833
EXAMPLES::
1834
1835
sage: R = Permutations(10).map(attrcall('reduced_word'))
1836
sage: R.cardinality()
1837
3628800
1838
1839
"""
1840
return self.cc.cardinality()
1841
1842
def __iter__(self):
1843
"""
1844
Returns an iterator over the elements of this combinatorial class
1845
1846
EXAMPLES::
1847
1848
sage: R = Permutations(10).map(attrcall('reduced_word'))
1849
sage: R.cardinality()
1850
3628800
1851
"""
1852
for x in self.cc:
1853
yield self.f(x)
1854
1855
def an_element(self):
1856
"""
1857
Returns an element of this combinatorial class
1858
1859
EXAMPLES::
1860
1861
sage: R = SymmetricGroup(10).map(attrcall('reduced_word'))
1862
sage: R.an_element()
1863
[9, 8, 7, 6, 5, 4, 3, 2, 1]
1864
"""
1865
return self.f(self.cc.an_element())
1866
1867
##############################################################################
1868
from sage.rings.all import infinity
1869
class InfiniteAbstractCombinatorialClass(CombinatorialClass):
1870
r"""
1871
This is an internal class that should not be used directly. A class which
1872
inherits from InfiniteAbstractCombinatorialClass inherits the standard
1873
methods list and count.
1874
1875
If self._infinite_cclass_slice exists then self.__iter__ returns an
1876
iterator for self, otherwise raise NotImplementedError. The method
1877
self._infinite_cclass_slice is supposed to accept any integer as an
1878
argument and return something which is iterable.
1879
"""
1880
def cardinality(self):
1881
"""
1882
Counts the elements of the combinatorial class.
1883
1884
EXAMPLES:
1885
sage: R = InfiniteAbstractCombinatorialClass()
1886
sage: R.cardinality()
1887
+Infinity
1888
"""
1889
return infinity
1890
1891
def list(self):
1892
"""
1893
Returns an error since self is an infinite combinatorial class.
1894
1895
EXAMPLES:
1896
sage: R = InfiniteAbstractCombinatorialClass()
1897
sage: R.list()
1898
Traceback (most recent call last):
1899
...
1900
NotImplementedError: infinite list
1901
"""
1902
raise NotImplementedError, "infinite list"
1903
1904
def __iter__(self):
1905
"""
1906
Returns an iterator for the infinite combinatorial class self if
1907
possible or raise a NotImplementedError.
1908
1909
EXAMPLES:
1910
sage: R = InfiniteAbstractCombinatorialClass()
1911
sage: iter(R).next()
1912
Traceback (most recent call last):
1913
...
1914
NotImplementedError
1915
1916
sage: c = iter(Compositions()) # indirect doctest
1917
sage: c.next(), c.next(), c.next(), c.next(), c.next(), c.next()
1918
([], [1], [1, 1], [2], [1, 1, 1], [1, 2])
1919
sage: c.next(), c.next(), c.next(), c.next(), c.next(), c.next()
1920
([2, 1], [3], [1, 1, 1, 1], [1, 1, 2], [1, 2, 1], [1, 3])
1921
"""
1922
try:
1923
finite = self._infinite_cclass_slice
1924
except AttributeError:
1925
raise NotImplementedError
1926
i = 0
1927
while True:
1928
for c in finite(i):
1929
yield c
1930
i+=1
1931
1932
1933
1934
def hurwitz_zeta(s,x,N):
1935
"""
1936
Returns the value of the `\zeta(s,x)` to `N`
1937
decimals, where s and x are real.
1938
1939
The Hurwitz zeta function is one of the many zeta functions. It
1940
defined as
1941
1942
.. math::
1943
1944
\zeta(s,x) = \sum_{k=0}^\infty (k+x)^{-s}.
1945
1946
1947
When `x = 1`, this coincides with Riemann's zeta function.
1948
The Dirichlet L-functions may be expressed as a linear combination
1949
of Hurwitz zeta functions.
1950
1951
Note that if you use floating point inputs, then the results may be
1952
slightly off.
1953
1954
EXAMPLES::
1955
1956
sage: hurwitz_zeta(3,1/2,6)
1957
8.41439000000000
1958
sage: hurwitz_zeta(11/10,1/2,6)
1959
12.1041000000000
1960
sage: hurwitz_zeta(11/10,1/2,50)
1961
12.10381349568375510570907741296668061903364861809
1962
1963
REFERENCES:
1964
1965
- http://en.wikipedia.org/wiki/Hurwitz_zeta_function
1966
"""
1967
maxima.eval('load ("bffac")')
1968
s = maxima.eval("bfhzeta (%s,%s,%s)"%(s,x,N))
1969
1970
#Handle the case where there is a 'b' in the string
1971
#'1.2000b0' means 1.2000 and
1972
#'1.2000b1' means 12.000
1973
i = s.rfind('b')
1974
if i == -1:
1975
return sage_eval(s)
1976
else:
1977
if s[i+1:] == '0':
1978
return sage_eval(s[:i])
1979
else:
1980
return sage_eval(s[:i])*10**sage_eval(s[i+1:])
1981
1982
return s ## returns an odd string
1983
1984
1985
#####################################################
1986
#### combinatorial sets/lists
1987
1988
def combinations(mset,k):
1989
r"""
1990
A combination of a multiset (a list of objects which may contain
1991
the same object several times) mset is an unordered selection
1992
without repetitions and is represented by a sorted sublist of mset.
1993
Returns the set of all combinations of the multiset mset with k
1994
elements.
1995
1996
.. warning::
1997
1998
Wraps GAP's Combinations. Hence mset must be a list of objects
1999
that have string representations that can be interpreted by the
2000
GAP interpreter. If mset consists of at all complicated Sage
2001
objects, this function does *not* do what you expect. A proper
2002
function should be written! (TODO!)
2003
2004
EXAMPLES::
2005
2006
sage: mset = [1,1,2,3,4,4,5]
2007
sage: combinations(mset,2)
2008
[[1, 1],
2009
[1, 2],
2010
[1, 3],
2011
[1, 4],
2012
[1, 5],
2013
[2, 3],
2014
[2, 4],
2015
[2, 5],
2016
[3, 4],
2017
[3, 5],
2018
[4, 4],
2019
[4, 5]]
2020
sage: mset = ["d","a","v","i","d"]
2021
sage: combinations(mset,3)
2022
['add', 'adi', 'adv', 'aiv', 'ddi', 'ddv', 'div']
2023
2024
.. note::
2025
2026
For large lists, this raises a string error.
2027
"""
2028
ans=gap.eval("Combinations(%s,%s)"%(mset,ZZ(k))).replace("\n","")
2029
return eval(ans)
2030
2031
def combinations_iterator(mset,n=None):
2032
"""
2033
Posted by Raymond Hettinger, 2006/03/23, to the Python Cookbook:
2034
http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/474124
2035
2036
Much faster than combinations.
2037
2038
EXAMPLES::
2039
2040
sage: X = combinations_iterator([1,2,3,4,5],3)
2041
sage: [x for x in X]
2042
[[1, 2, 3],
2043
[1, 2, 4],
2044
[1, 2, 5],
2045
[1, 3, 4],
2046
[1, 3, 5],
2047
[1, 4, 5],
2048
[2, 3, 4],
2049
[2, 3, 5],
2050
[2, 4, 5],
2051
[3, 4, 5]]
2052
"""
2053
items = mset
2054
if n is None:
2055
n = len(items)
2056
for i in range(len(items)):
2057
v = items[i:i+1]
2058
if n == 1:
2059
yield v
2060
else:
2061
rest = items[i+1:]
2062
for c in combinations_iterator(rest, n-1):
2063
yield v + c
2064
2065
2066
def number_of_combinations(mset,k):
2067
"""
2068
Returns the size of combinations(mset,k). IMPLEMENTATION: Wraps
2069
GAP's NrCombinations.
2070
2071
.. note::
2072
2073
``mset`` must be a list of integers or strings (i.e., this is
2074
very restricted).
2075
2076
EXAMPLES::
2077
2078
sage: mset = [1,1,2,3,4,4,5]
2079
sage: number_of_combinations(mset,2)
2080
12
2081
"""
2082
return ZZ(gap.eval("NrCombinations(%s,%s)"%(mset,ZZ(k))))
2083
2084
def arrangements(mset,k):
2085
r"""
2086
An arrangement of mset is an ordered selection without repetitions
2087
and is represented by a list that contains only elements from mset,
2088
but maybe in a different order.
2089
2090
``arrangements`` returns the set of arrangements of the
2091
multiset mset that contain k elements.
2092
2093
IMPLEMENTATION: Wraps GAP's Arrangements.
2094
2095
.. warning::
2096
2097
Wraps GAP - hence mset must be a list of objects that have
2098
string representations that can be interpreted by the GAP
2099
interpreter. If mset consists of at all complicated Sage
2100
objects, this function does *not* do what you expect. A proper
2101
function should be written! (TODO!)
2102
2103
EXAMPLES::
2104
2105
sage: mset = [1,1,2,3,4,4,5]
2106
sage: arrangements(mset,2)
2107
[[1, 1],
2108
[1, 2],
2109
[1, 3],
2110
[1, 4],
2111
[1, 5],
2112
[2, 1],
2113
[2, 3],
2114
[2, 4],
2115
[2, 5],
2116
[3, 1],
2117
[3, 2],
2118
[3, 4],
2119
[3, 5],
2120
[4, 1],
2121
[4, 2],
2122
[4, 3],
2123
[4, 4],
2124
[4, 5],
2125
[5, 1],
2126
[5, 2],
2127
[5, 3],
2128
[5, 4]]
2129
sage: arrangements( ["c","a","t"], 2 )
2130
['ac', 'at', 'ca', 'ct', 'ta', 'tc']
2131
sage: arrangements( ["c","a","t"], 3 )
2132
['act', 'atc', 'cat', 'cta', 'tac', 'tca']
2133
"""
2134
ans=gap.eval("Arrangements(%s,%s)"%(mset,k))
2135
return eval(ans)
2136
2137
def number_of_arrangements(mset,k):
2138
"""
2139
Returns the size of arrangements(mset,k). Wraps GAP's
2140
NrArrangements.
2141
2142
EXAMPLES::
2143
2144
sage: mset = [1,1,2,3,4,4,5]
2145
sage: number_of_arrangements(mset,2)
2146
22
2147
"""
2148
return ZZ(gap.eval("NrArrangements(%s,%s)"%(mset,ZZ(k))))
2149
2150
def derangements(mset):
2151
"""
2152
A derangement is a fixed point free permutation of list and is
2153
represented by a list that contains exactly the same elements as
2154
mset, but possibly in different order. Derangements returns the set
2155
of all derangements of a multiset.
2156
2157
Wraps GAP's Derangements.
2158
2159
.. warning::
2160
2161
Wraps GAP - hence mset must be a list of objects that have
2162
string representations that can be interpreted by the GAP
2163
interpreter. If mset consists of at all complicated Sage
2164
objects, this function does *not* do what you expect. A proper
2165
function should be written! (TODO!)
2166
2167
EXAMPLES::
2168
2169
sage: mset = [1,2,3,4]
2170
sage: derangements(mset)
2171
[[2, 1, 4, 3],
2172
[2, 3, 4, 1],
2173
[2, 4, 1, 3],
2174
[3, 1, 4, 2],
2175
[3, 4, 1, 2],
2176
[3, 4, 2, 1],
2177
[4, 1, 2, 3],
2178
[4, 3, 1, 2],
2179
[4, 3, 2, 1]]
2180
sage: derangements(["c","a","t"])
2181
['atc', 'tca']
2182
"""
2183
ans=gap.eval("Derangements(%s)"%mset)
2184
return eval(ans)
2185
2186
def number_of_derangements(mset):
2187
"""
2188
Returns the size of derangements(mset). Wraps GAP's
2189
NrDerangements.
2190
2191
EXAMPLES::
2192
2193
sage: mset = [1,2,3,4]
2194
sage: number_of_derangements(mset)
2195
9
2196
"""
2197
ans=gap.eval("NrDerangements(%s)"%mset)
2198
return ZZ(ans)
2199
2200
def tuples(S,k):
2201
"""
2202
An ordered tuple of length k of set is an ordered selection with
2203
repetition and is represented by a list of length k containing
2204
elements of set. tuples returns the set of all ordered tuples of
2205
length k of the set.
2206
2207
EXAMPLES::
2208
2209
sage: S = [1,2]
2210
sage: tuples(S,3)
2211
[[1, 1, 1], [2, 1, 1], [1, 2, 1], [2, 2, 1], [1, 1, 2], [2, 1, 2], [1, 2, 2], [2, 2, 2]]
2212
sage: mset = ["s","t","e","i","n"]
2213
sage: tuples(mset,2)
2214
[['s', 's'], ['t', 's'], ['e', 's'], ['i', 's'], ['n', 's'], ['s', 't'], ['t', 't'],
2215
['e', 't'], ['i', 't'], ['n', 't'], ['s', 'e'], ['t', 'e'], ['e', 'e'], ['i', 'e'],
2216
['n', 'e'], ['s', 'i'], ['t', 'i'], ['e', 'i'], ['i', 'i'], ['n', 'i'], ['s', 'n'],
2217
['t', 'n'], ['e', 'n'], ['i', 'n'], ['n', 'n']]
2218
2219
The Set(...) comparisons are necessary because finite fields are
2220
not enumerated in a standard order.
2221
2222
::
2223
2224
sage: K.<a> = GF(4, 'a')
2225
sage: mset = [x for x in K if x!=0]
2226
sage: tuples(mset,2)
2227
[[a, a], [a + 1, a], [1, a], [a, a + 1], [a + 1, a + 1], [1, a + 1], [a, 1], [a + 1, 1], [1, 1]]
2228
2229
AUTHORS:
2230
2231
- Jon Hanke (2006-08)
2232
"""
2233
import copy
2234
if k<=0:
2235
return [[]]
2236
if k==1:
2237
return [[x] for x in S]
2238
ans = []
2239
for s in S:
2240
for x in tuples(S,k-1):
2241
y = copy.copy(x)
2242
y.append(s)
2243
ans.append(y)
2244
return ans
2245
## code wrapping GAP's Tuples:
2246
#ans=gap.eval("Tuples(%s,%s)"%(S,k))
2247
#return eval(ans)
2248
2249
2250
def number_of_tuples(S,k):
2251
"""
2252
Returns the size of tuples(S,k). Wraps GAP's NrTuples.
2253
2254
EXAMPLES::
2255
2256
sage: S = [1,2,3,4,5]
2257
sage: number_of_tuples(S,2)
2258
25
2259
sage: S = [1,1,2,3,4,5]
2260
sage: number_of_tuples(S,2)
2261
25
2262
"""
2263
ans=gap.eval("NrTuples(%s,%s)"%(S,ZZ(k)))
2264
return ZZ(ans)
2265
2266
def unordered_tuples(S,k):
2267
"""
2268
An unordered tuple of length k of set is a unordered selection with
2269
repetitions of set and is represented by a sorted list of length k
2270
containing elements from set.
2271
2272
unordered_tuples returns the set of all unordered tuples of length
2273
k of the set. Wraps GAP's UnorderedTuples.
2274
2275
.. warning::
2276
2277
Wraps GAP - hence mset must be a list of objects that have
2278
string representations that can be interpreted by the GAP
2279
interpreter. If mset consists of at all complicated Sage
2280
objects, this function does *not* do what you expect. A proper
2281
function should be written! (TODO!)
2282
2283
EXAMPLES::
2284
2285
sage: S = [1,2]
2286
sage: unordered_tuples(S,3)
2287
[[1, 1, 1], [1, 1, 2], [1, 2, 2], [2, 2, 2]]
2288
sage: unordered_tuples(["a","b","c"],2)
2289
['aa', 'ab', 'ac', 'bb', 'bc', 'cc']
2290
"""
2291
ans=gap.eval("UnorderedTuples(%s,%s)"%(S,ZZ(k)))
2292
return eval(ans)
2293
2294
def number_of_unordered_tuples(S,k):
2295
"""
2296
Returns the size of unordered_tuples(S,k). Wraps GAP's
2297
NrUnorderedTuples.
2298
2299
EXAMPLES::
2300
2301
sage: S = [1,2,3,4,5]
2302
sage: number_of_unordered_tuples(S,2)
2303
15
2304
"""
2305
ans=gap.eval("NrUnorderedTuples(%s,%s)"%(S,ZZ(k)))
2306
return ZZ(ans)
2307
2308
def permutations(mset):
2309
"""
2310
A permutation is represented by a list that contains exactly the
2311
same elements as mset, but possibly in different order. If mset is
2312
a proper set there are `|mset| !` such permutations.
2313
Otherwise if the first elements appears `k_1` times, the
2314
second element appears `k_2` times and so on, the number
2315
of permutations is `|mset|! / (k_1! k_2! \ldots)`, which
2316
is sometimes called a multinomial coefficient.
2317
2318
permutations returns the set of all permutations of a multiset.
2319
Calls a function written by Mike Hansen, not GAP.
2320
2321
EXAMPLES::
2322
2323
sage: mset = [1,1,2,2,2]
2324
sage: permutations(mset)
2325
[[1, 1, 2, 2, 2],
2326
[1, 2, 1, 2, 2],
2327
[1, 2, 2, 1, 2],
2328
[1, 2, 2, 2, 1],
2329
[2, 1, 1, 2, 2],
2330
[2, 1, 2, 1, 2],
2331
[2, 1, 2, 2, 1],
2332
[2, 2, 1, 1, 2],
2333
[2, 2, 1, 2, 1],
2334
[2, 2, 2, 1, 1]]
2335
sage: MS = MatrixSpace(GF(2),2,2)
2336
sage: A = MS([1,0,1,1])
2337
sage: permutations(A.rows())
2338
[[(1, 0), (1, 1)], [(1, 1), (1, 0)]]
2339
"""
2340
from sage.combinat.permutation import Permutations
2341
ans = Permutations(mset)
2342
return ans.list()
2343
2344
def permutations_iterator(mset,n=None):
2345
"""
2346
Do not use this function. It will be deprecated in future version
2347
of Sage and eventually removed. Use Permutations instead; instead
2348
of
2349
2350
for p in permutations_iterator(range(1, m+1), n)
2351
2352
use
2353
2354
for p in Permutations(m, n).
2355
2356
Note that Permutations, unlike this function, treats repeated
2357
elements as identical.
2358
2359
If you insist on using this now:
2360
2361
Returns an iterator (http://docs.python.org/lib/typeiter.html)
2362
which can be used in place of permutations(mset) if all you need it
2363
for is a 'for' loop.
2364
2365
Posted by Raymond Hettinger, 2006/03/23, to the Python Cookbook:
2366
http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/474124
2367
2368
Note- This function considers repeated elements as different
2369
entries, so for example::
2370
2371
sage: from sage.combinat.combinat import permutations, permutations_iterator
2372
sage: mset = [1,2,2]
2373
sage: permutations(mset)
2374
[[1, 2, 2], [2, 1, 2], [2, 2, 1]]
2375
sage: for p in permutations_iterator(mset): print p
2376
[1, 2, 2]
2377
[1, 2, 2]
2378
[2, 1, 2]
2379
[2, 2, 1]
2380
[2, 1, 2]
2381
[2, 2, 1]
2382
2383
EXAMPLES::
2384
2385
sage: X = permutations_iterator(range(3),2)
2386
sage: [x for x in X]
2387
[[0, 1], [0, 2], [1, 0], [1, 2], [2, 0], [2, 1]]
2388
"""
2389
items = mset
2390
if n is None:
2391
n = len(items)
2392
for i in range(len(items)):
2393
v = items[i:i+1]
2394
if n == 1:
2395
yield v
2396
else:
2397
rest = items[:i] + items[i+1:]
2398
for p in permutations_iterator(rest, n-1):
2399
yield v + p
2400
2401
def number_of_permutations(mset):
2402
"""
2403
Do not use this function. It will be deprecated in future version
2404
of Sage and eventually removed. Use Permutations instead; instead
2405
of
2406
2407
number_of_permutations(mset)
2408
2409
use
2410
2411
Permutations(mset).cardinality().
2412
2413
If you insist on using this now:
2414
2415
Returns the size of permutations(mset).
2416
2417
AUTHORS:
2418
2419
- Robert L. Miller
2420
2421
EXAMPLES::
2422
2423
sage: mset = [1,1,2,2,2]
2424
sage: number_of_permutations(mset)
2425
10
2426
"""
2427
from sage.rings.arith import factorial
2428
m = len(mset)
2429
n = []
2430
seen = []
2431
for element in mset:
2432
try:
2433
n[seen.index(element)] += 1
2434
except (IndexError, ValueError):
2435
n.append(1)
2436
seen.append(element)
2437
return factorial(m)/prod([factorial(k) for k in n])
2438
2439
def cyclic_permutations(mset):
2440
"""
2441
Returns a list of all cyclic permutations of mset. Treats mset as a
2442
list, not a set, i.e. entries with the same value are distinct.
2443
2444
AUTHORS:
2445
2446
- Emily Kirkman
2447
2448
EXAMPLES::
2449
2450
sage: from sage.combinat.combinat import cyclic_permutations, cyclic_permutations_iterator
2451
sage: cyclic_permutations(range(4))
2452
[[0, 1, 2, 3], [0, 1, 3, 2], [0, 2, 1, 3], [0, 2, 3, 1], [0, 3, 1, 2], [0, 3, 2, 1]]
2453
sage: for cycle in cyclic_permutations(['a', 'b', 'c']):
2454
... print cycle
2455
['a', 'b', 'c']
2456
['a', 'c', 'b']
2457
2458
Note that lists with repeats are not handled intuitively::
2459
2460
sage: cyclic_permutations([1,1,1])
2461
[[1, 1, 1], [1, 1, 1]]
2462
"""
2463
return list(cyclic_permutations_iterator(mset))
2464
2465
def cyclic_permutations_iterator(mset):
2466
"""
2467
Iterates over all cyclic permutations of mset in cycle notation.
2468
Treats mset as a list, not a set, i.e. entries with the same value
2469
are distinct.
2470
2471
AUTHORS:
2472
2473
- Emily Kirkman
2474
2475
EXAMPLES::
2476
2477
sage: from sage.combinat.combinat import cyclic_permutations, cyclic_permutations_iterator
2478
sage: cyclic_permutations(range(4))
2479
[[0, 1, 2, 3], [0, 1, 3, 2], [0, 2, 1, 3], [0, 2, 3, 1], [0, 3, 1, 2], [0, 3, 2, 1]]
2480
sage: for cycle in cyclic_permutations(['a', 'b', 'c']):
2481
... print cycle
2482
['a', 'b', 'c']
2483
['a', 'c', 'b']
2484
2485
Note that lists with repeats are not handled intuitively::
2486
2487
sage: cyclic_permutations([1,1,1])
2488
[[1, 1, 1], [1, 1, 1]]
2489
"""
2490
if len(mset) > 2:
2491
for perm in permutations_iterator(mset[1:]):
2492
yield [mset[0]] + perm
2493
else:
2494
yield mset
2495
2496
def bell_polynomial(n, k):
2497
r"""
2498
This function returns the Bell Polynomial
2499
2500
.. math::
2501
2502
B_{n,k}(x_1, x_2, \ldots, x_{n-k+1}) = \sum_{\sum{j_i}=k, \sum{i j_i}
2503
=n} \frac{n!}{j_1!j_2!\ldots} \frac{x_1}{1!}^j_1 \frac{x_2}{2!}^j_2
2504
\ldots
2505
2506
INPUT:
2507
2508
- ``n`` - integer
2509
2510
- ``k`` - integer
2511
2512
OUTPUT:
2513
2514
- polynomial expression (SymbolicArithmetic)
2515
2516
EXAMPLES::
2517
2518
sage: bell_polynomial(6,2)
2519
10*x_3^2 + 15*x_2*x_4 + 6*x_1*x_5
2520
sage: bell_polynomial(6,3)
2521
15*x_2^3 + 60*x_1*x_2*x_3 + 15*x_1^2*x_4
2522
2523
REFERENCES:
2524
2525
- E.T. Bell, "Partition Polynomials"
2526
2527
AUTHORS:
2528
2529
- Blair Sutton (2009-01-26)
2530
"""
2531
from sage.combinat.partition import Partitions
2532
from sage.rings.arith import factorial
2533
vars = ZZ[tuple(['x_'+str(i) for i in range(1, n-k+2)])].gens()
2534
result = 0
2535
for p in Partitions(n, length=k):
2536
factorial_product = 1
2537
power_factorial_product = 1
2538
for part, count in p.to_exp_dict().iteritems():
2539
factorial_product *= factorial(count)
2540
power_factorial_product *= factorial(part)**count
2541
2542
coefficient = factorial(n) / (factorial_product * power_factorial_product)
2543
result += coefficient * prod([vars[i-1] for i in p])
2544
2545
return result
2546
2547
def fibonacci_sequence(start, stop=None, algorithm=None):
2548
r"""
2549
Returns an iterator over the Fibonacci sequence, for all fibonacci
2550
numbers `f_n` from ``n = start`` up to (but
2551
not including) ``n = stop``
2552
2553
INPUT:
2554
2555
2556
- ``start`` - starting value
2557
2558
- ``stop`` - stopping value
2559
2560
- ``algorithm`` - default (None) - passed on to
2561
fibonacci function (or not passed on if None, i.e., use the
2562
default).
2563
2564
2565
EXAMPLES::
2566
2567
sage: fibs = [i for i in fibonacci_sequence(10, 20)]
2568
sage: fibs
2569
[55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181]
2570
2571
::
2572
2573
sage: sum([i for i in fibonacci_sequence(100, 110)])
2574
69919376923075308730013
2575
2576
.. seealso::
2577
2578
:func:`fibonacci_xrange`
2579
2580
AUTHORS:
2581
2582
- Bobby Moretti
2583
"""
2584
if stop is None:
2585
stop = ZZ(start)
2586
start = ZZ(0)
2587
else:
2588
start = ZZ(start)
2589
stop = ZZ(stop)
2590
2591
if algorithm:
2592
for n in xrange(start, stop):
2593
yield fibonacci(n, algorithm=algorithm)
2594
else:
2595
for n in xrange(start, stop):
2596
yield fibonacci(n)
2597
2598
def fibonacci_xrange(start, stop=None, algorithm='pari'):
2599
r"""
2600
Returns an iterator over all of the Fibonacci numbers in the given
2601
range, including ``f_n = start`` up to, but not
2602
including, ``f_n = stop``.
2603
2604
EXAMPLES::
2605
2606
sage: fibs_in_some_range = [i for i in fibonacci_xrange(10^7, 10^8)]
2607
sage: len(fibs_in_some_range)
2608
4
2609
sage: fibs_in_some_range
2610
[14930352, 24157817, 39088169, 63245986]
2611
2612
::
2613
2614
sage: fibs = [i for i in fibonacci_xrange(10, 100)]
2615
sage: fibs
2616
[13, 21, 34, 55, 89]
2617
2618
::
2619
2620
sage: list(fibonacci_xrange(13, 34))
2621
[13, 21]
2622
2623
A solution to the second Project Euler problem::
2624
2625
sage: sum([i for i in fibonacci_xrange(10^6) if is_even(i)])
2626
1089154
2627
2628
.. seealso::
2629
2630
:func:`fibonacci_sequence`
2631
2632
AUTHORS:
2633
2634
- Bobby Moretti
2635
"""
2636
if stop is None:
2637
stop = ZZ(start)
2638
start = ZZ(0)
2639
else:
2640
start = ZZ(start)
2641
stop = ZZ(stop)
2642
2643
# iterate until we've gotten high enough
2644
fn = 0
2645
n = 0
2646
while fn < start:
2647
n += 1
2648
fn = fibonacci(n)
2649
2650
while True:
2651
fn = fibonacci(n)
2652
n += 1
2653
if fn < stop:
2654
yield fn
2655
else:
2656
return
2657
2658
def bernoulli_polynomial(x, n):
2659
r"""
2660
Return the nth Bernoulli polynomial evaluated at x.
2661
2662
The generating function for the Bernoulli polynomials is
2663
2664
.. math::
2665
2666
\frac{t e^{xt}}{e^t-1}= \sum_{n=0}^\infty B_n(x) \frac{t^n}{n!},
2667
2668
and they are given directly by
2669
2670
.. math::
2671
2672
B_n(x) = \sum_{i=0}^n \binom{n}{i}B_{n-i}x^i.
2673
2674
One has `B_n(x) = - n\zeta(1 - n,x)`, where
2675
`\zeta(s,x)` is the Hurwitz zeta function. Thus, in a
2676
certain sense, the Hurwitz zeta function generalizes the
2677
Bernoulli polynomials to non-integer values of n.
2678
2679
EXAMPLES::
2680
2681
sage: y = QQ['y'].0
2682
sage: bernoulli_polynomial(y, 5)
2683
y^5 - 5/2*y^4 + 5/3*y^3 - 1/6*y
2684
sage: bernoulli_polynomial(y, 5)(12)
2685
199870
2686
sage: bernoulli_polynomial(12, 5)
2687
199870
2688
sage: bernoulli_polynomial(y^2 + 1, 5)
2689
y^10 + 5/2*y^8 + 5/3*y^6 - 1/6*y^2
2690
sage: P.<t> = ZZ[]
2691
sage: p = bernoulli_polynomial(t, 6)
2692
sage: p.parent()
2693
Univariate Polynomial Ring in t over Rational Field
2694
2695
We verify an instance of the formula which is the origin of
2696
the Bernoulli polynomials (and numbers)::
2697
2698
sage: power_sum = sum(k^4 for k in range(10))
2699
sage: 5*power_sum == bernoulli_polynomial(10, 5) - bernoulli(5)
2700
True
2701
2702
2703
REFERENCES:
2704
2705
- http://en.wikipedia.org/wiki/Bernoulli_polynomials
2706
"""
2707
try:
2708
n = ZZ(n)
2709
if n < 0:
2710
raise TypeError
2711
except TypeError:
2712
raise ValueError, "The second argument must be a non-negative integer"
2713
2714
if n == 0:
2715
return ZZ(1)
2716
2717
if n == 1:
2718
return x - ZZ(1)/2
2719
2720
k = n.mod(2)
2721
coeffs = [0]*k + sum(([binomial(n, i)*bernoulli(n-i), 0]
2722
for i in range(k, n+1, 2)), [])
2723
coeffs[-3] = -n/2
2724
2725
if isinstance(x, Polynomial):
2726
try:
2727
return x.parent()(coeffs)(x)
2728
except TypeError:
2729
pass
2730
2731
x2 = x*x
2732
xi = x**k
2733
s = 0
2734
for i in range(k, n-1, 2):
2735
s += coeffs[i]*xi
2736
t = xi
2737
xi *= x2
2738
s += xi - t*x*n/2
2739
return s
2740
2741