## "Fossies" - the Fresh Open Source Software Archive

### Member "NZMATH-1.2.0/nzmath/combinatorial.py" (19 Nov 2012, 18741 Bytes) of package /linux/misc/old/NZMATH-1.2.0.tar.gz:

As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) Python source code syntax highlighting (style: standard) with prefixed line numbers. Alternatively you can here view or download the uninterpreted source code file. For more information about "combinatorial.py" see the Fossies "Dox" file reference documentation.

    1 """
2 Combinatorial functions
3 """
4
5 import bisect
6 import itertools
7 from nzmath.rational import Integer, Rational
8
9
10 def binomial(n, m):
11     """
12     The binomial coefficient.
13     binomial(n, m) returns n ! / ((n - m) ! * m !).
14
15     n must be a positive integer and m must be a non-negative integer.
16     For convinience, binomial(n, n+i) = 0 for positive i, and
17     binomial(0,0) = 1.
18
19     In other cases, it raises an exception.
20     """
21     if not isinstance(n, (int, long)):
22         raise TypeError("integer is expected, %s detected." % n.__class__)
23     if not isinstance(m, (int, long)):
24         raise TypeError("integer is expected, %s detected." % m.__class__)
25     if n == m >= 0 or m == 0 and n > 0:
26         return 1
27     if n <= 0:
28         raise ValueError("non-positive number: %d" % n)
29     if m < 0:
30         raise ValueError("negative number: %d" % m)
31     if n < m:
32         return 0
33     if m*2 > n:
34         m = n - m
35     retval = n
36     for i in range(1, m):
37         retval *= (n - i)
38         retval /= (i + 1)
39     return Integer(retval)
40
41 def factorial(n):
42     """
43     Return n! for non negative integer n.
44     """
45     if not isinstance(n, (int, long)):
46         raise TypeError("integer is expected, %s detected." % n.__class__)
47     elif n < 0:
48         raise ValueError("argument must not be a negative integer.")
49     elif n < 1500:
50         # naive loop is best for small n.
51         result = 1
52         for i in xrange(2, n + 1):
53             result *= i
54         return Integer(result)
55     # The following algorithm keeps temporary results rather small as
56     # possible in order to make the function run faster than the naive
57     # loop.
58     l = range(1, n + 1)
59     while len(l) > 1:
60         for i in xrange(len(l) >> 1):
61             l[i] *= l.pop()
62     return Integer(l.pop())
63
64 def bernoulli(n):
65     """
66     Return n-th Bernoulli number.
67     """
68     if n != 1 and n & 1:
69         return Integer(0)
70     B = {0:Integer(1),
71          1:Rational(-1, 2)}
72     for i in range(2, n + 1, 2):
73         a = B[0] + (i + 1) * B[1]
74         for j in range(2, i, 2):
75             a += binomial(i + 1, j) * B[j]
76         B[i] = -a / (i + 1)
77     return B[n]
78
79 def catalan(n):
80     """
81     Return n-th Catalan number.
82     """
83     return binomial(2*n, n) // (n+1)
84
85 def euler(n):
86     """
87     Return n-th Euler number.
88     """
89     if n & 1:
90         return Integer(0)
91     E = {0:Integer(1)}
92     for i in range(2, n + 1, 2):
93         E[i] = sum([-binomial(i, j) * E[j] for j in range(0, i, 2)])
94     return E[n]
95
96 def fallingfactorial(n, m):
97     """
98     Return the falling factorial; n to the m falling, i.e. n(n-1)..(n-m+1).
99
100     For Example:
101     >>> fallingfactorial(7, 3)
102     210
103     """
104     r = 1
105     for i in range(n, n-m, -1):
106         r *= i
107     return r
108
109 def risingfactorial(n, m):
110     """
111     Return the rising factorial; n to the m rising, i.e. n(n+1)..(n+m-1).
112
113     For example:
114     >>> risingfactorial(7, 3)
115     504
116     """
117     r = 1
118     for i in range(n, n+m):
119         r *= i
120     return r
121
122 def multinomial(n, parts):
123     """
124     Return multinomial coefficient.
125
126     parts MUST be a sequence of natural numbers and n==sum(parts) holds.
127     """
128     if n != sum(parts):
129         raise ValueError("sum of parts must be equal to n.")
130     for part in parts:
131         if not isinstance(part, (int, long)) or part < 0:
132             raise ValueError("parts must be a sequence of natural numbers.")
133     # binomial case
134     if len(parts) == 2:
135         return binomial(n, parts[0])
136     # other cases
137     result = factorial(n)
138     for part in parts:
139         if part >= 2: # 0! = 1! = 1 are negligible
140             result //= factorial(part)
141     return result
142
143 def stirling1(n, m):
144     """
145     Stirling number of the first kind.
146
147     Let s denote the Stirling number, (x)_n falling factorial, then
148       (x)_n = \sum_{i=0}^{n} s(n, i) * x**i
149     and s satisfies the recurrence relation:
150       s(n, m) = s(n-1, m-1) - (n-1)*s(n-1, m)
151     """
152     if n == m:
153         return 1
154     elif n < m or n*m < 0:
155         return 0
156     elif m == 0:
157         return 0
158     elif m == 1:
159         if n & 1:
160             sign = 1
161         else:
162             sign = -1
163         return sign * factorial(n - 1)
164     elif m == n - 1:
165         return -binomial(n, 2)
166     else:
167         # recursively calculate by s(n, m) = s(n-1, m-1) - (n-1)*S(n-1, m)
168         slist = (1,) # S(1, 1) = 1
169         for i in range(1, n):
170             l, u = max(1, m - n + i + 1), min(i + 2, m + 1)
171             if l == 1 and len(slist) < n - l:
172                 nlist = [-i * slist[0]]
173             else:
174                 nlist = []
175             for j in range(len(slist) - 1):
176                 nlist.append(slist[j] - i * slist[j + 1])
177             if len(slist) <= u - l:
178                 nlist.append(slist[-1])
179             slist = tuple(nlist)
180         return slist[0]
181
182 def stirling2(n, m):
183     """
184     Stirling number of the second kind.
185
186     Let S denote the Stirling number, (x)_i falling factorial, then:
187       x**n = \sum_{i=0}^{n} S(n, i) * (x)_i
188     S satisfies:
189       S(n, m) = S(n-1, m-1) + m*S(n-1, m)
190     """
191     if n == m:
192         return 1
193     elif n < m or n * m <= 0:
194         return 0
195     elif n < 0:
196         return stirling2_negative(-m, -n)
197     elif m == 1:
198         return 1
199     elif m == 2:
200         return 2**(n - 1) - 1
201     elif m == n - 1:
202         return n * m >> 1
203     else:
204         r = 0
205         b = m
206         l = 1
207         while 1:
208             if l & 1:
209                 r -= b * l**n
210             else:
211                 r += b * l**n
212             if l == m:
213                 break
214             b = (b * (m - l)) // (l + 1)
215             l += 1
216         if m & 1:
217             r = -r
218         return r // factorial(m)
219
220 def stirling2_negative(n, m):
221     """
222     Stiring number of the second kind extended to negative numbers.
223
224     Let S# denote the extended Stirling number, S the original, then
225     S#(n, m) = S(-m, -n) and extended by the recurrence relation:
226       S#(n, m) = S#(n-1, m-1) + (n-1)*S#(n-1, m)
227     """
228     if n == m:
229         return 1
230     elif n < m or n * m <= 0:
231         return 0
232     elif n < 0:
233         return stirling2(-m, -n)
234     elif m == 1:
235         return factorial(n - 1)
236     else: # return S(n-1,m-1)+(n-1)*S(n-1,m)
237         slist = [1]
238         i = 1
239         while i < n:
240             l, u = max(1, m - n + i + 1), min(i + 2, m + 1)
241             if len(slist) < u - l:
242                 nlist = [i * slist[0]]
243             else:
244                 nlist = []
245             for j in range(len(slist) - 1):
246                 nlist.append(slist[j] + i * slist[j + 1])
247             if len(slist) <= u - l:
248                 nlist.append(slist[-1])
249             slist = nlist
250             i += 1
251         return slist[-1]
252
253 def bell(n):
254     """
255     Bell number.
256
257     The Bell number b is defined by:
258       b(n) = \sum_{i=0}^{n} S(n, i)
259     where S denotes Stirling number of the second kind.
260     """
261     return sum([stirling2(n, i) for i in range(n + 1)])
262
263
264 # generators
265 def combination_index_generator(n, m):
266     """
267     Generate indices of m elment subsets of n element set.
268
269     The number of generated indices is binomial(n, m).
270
271     For example,
272     combinationIndexGenerator(5,3) generates the following lists:
273         [0, 1, 2]
274         [0, 1, 3]
275         [0, 1, 4]
276         [0, 2, 3]
277         [0, 2, 4]
278         [0, 3, 4]
279         [1, 2, 3]
280         [1, 2, 4]
281         [1, 3, 4]
282         [2, 3, 4]
283     """
284     assert n >= m > 0
285     idx = range(m)
286     while True:
287         yield list(idx)
288         for i in range(1, m+1):
289             if idx[-i] != n-i:
290                 idx[-i] += 1
291                 for j in range(i-1, 0, -1):
292                     idx[-j] = idx[-j-1] + 1
293                 break
294         else:
295             raise StopIteration
296
297
298 def permutation_generator(n):
299     """
300     Generate all permutations of n elements as lists.
301
302     The number of generated lists is n!, so be careful to use big n.
303
304     For example,
305     permutationGenerator(3) generates the following lists:
306         [0, 1, 2]
307         [0, 2, 1]
308         [1, 0, 2]
309         [1, 2, 0]
310         [2, 0, 1]
311         [2, 1, 0]
312     """
313     # traverse tree by depth first
314     perm = range(n)
315     unused = []
316     while True:
317         # leaf is reached, thus yield the value.
318         yield list(perm)
319
320         # track back until node with subtree yet to be traversed
321         last = n - 1
322         unused.append(perm[-1])
323         while last and perm[last - 1] > unused[-1]:
324             last -= 1
325             unused.append(perm[last])
326
327         # exhausted
328         if not last:
329             break
330
331         # assert unused == sorted(unused)
332         # replace with just bigger than perm[last - 1]
333         index = bisect.bisect(unused, perm[last - 1])
334         unused[index], perm[last - 1] = perm[last - 1], unused[index]
335         # replace remaining part
336         perm[last:] = unused
337         del unused[:]
338
339
340 def dyck_word_generator(n, alphabet=(0, 1)):
341     """
342     Generate all Dyck words of length 2*n as tuples.
343
344     The Dyck words are words on a two character alphabet.
345     The number of each character in a word is equal,
346     and the number of the second character never exceeds the first
347     in any initial parts of the word.
348
349     The number of generated words is the n-th Catalan number.
350
351     The alphabet is {0, 1} by default, but you can pass it into the
352     optional argument 'alphabet'.
353
354     For example,
355     >>> for word in dyck_word_generator(3, alphabet=("(", ")")):
356     ...     print "".join(word)
357     ...
358     ()()()
359     ()(())
360     (())()
361     (()())
362     ((()))
363     >>>
364     """
365     if n == 0:
366         yield ()
367         raise StopIteration()
368     else:
369         # memo[i] is a list of Dyck words of length 2*i
370         memo = [[()]]
371         alpha, beta = (alphabet[0],), (alphabet[1],)
372
373     # prepare up to n-1
374     for i in range(1, n):
375         memo.append([])
376         for j in range(i):
377             for p in memo[j]:
378                 prefix = alpha + p + beta
379                 for q in memo[i - j - 1]:
380                     memo[i].append(prefix + q)
381
382     # output
383     for j in range(n):
384         for p in memo[j]:
385             prefix = alpha + p + beta
386             for q in memo[n - j - 1]:
387                 yield prefix + q
388
389
390 # partition related functions and classes
391 def partition_generator(n, maxi=None):
392     """
393     Generate partitions of n.
394     If maxi is given, then parts are limited to at most maxi.
395     """
396     if not n:
397         yield ()
398         raise StopIteration
399     if maxi is None or maxi > n:
400         maxi = n
401     partition = [maxi]
402     rest = n - maxi
403     while True:
404         key = partition[-1]
405         while rest >= key:
406             partition.append(key)
407             rest -= key
408         if rest:
409             partition.append(rest)
410
411         yield tuple(partition)
412
413         try:
414             # wind up all 1's.
415             first_one = partition.index(1)
416             if not first_one:
417                 # partition was [1]*n means all partitions have been generated.
418                 raise StopIteration
419             rest = len(partition) - first_one
420             del partition[first_one:]
421         except ValueError:
423             rest = 0
424
425         partition[-1] -= 1
426         rest += 1
427
428
429 class PartitionDriver(object):
430     """
431     PartitionDriver knows how to construct partitions.
432
433     This class can generate all partitions without any additional
434     condition, yet it should be usable as a base class for subclasses
435     to generate conditioned partitions, such as limited maximum, with
436     odd parts, with distinct parts, etc.
437     """
438     def __init__(self):
439         """
440         Initialize the class.
441
442         Public attributes are partition and rest, though they are not
443         expected to be modified directly.
444         """
445         self.partition = []
446         self.rest = 0
447
448     def set_parameters(self, integer):
449         """
450         Set parameters for generate:
451           integer: the number to be partitioned
452
453         The method signature may differ class to class.
454         """
455         self.partition = []
456         self.rest = integer
457
458     def ispreconditioned(self):
459         """
460         Check whether preconditions are satisfied or not.
461
462         For example, partition into evens needs rest be even.
463         """
464         return True
465
466     def pull(self):
467         """
468         Pull the maximum part obtainable from rest to partition.
469
470         When you override the method, please try to make it precise.
471         If it is unavoidable to do some failing attempt, raising a
472         LookupError is the signal to tell that no parts can be pulled.
473         """
474         if self.partition:
475             if self.rest >= self.partition[-1]:
476                 part = self.partition[-1]
477             else:
478                 part = self.rest
479         else:
480             part = self.rest
481
482         self.partition.append(part)
483         self.rest -= part
484
485     def backtrack(self):
486         """
487         Back track the history from either broken or complete partition.
488         """
489         try:
490             if self.partition[-1] == 1:
491                 # wind up all 1's.
492                 first_one = self.partition.index(1)
493                 self.rest = len(self.partition) - first_one
494                 del self.partition[first_one:]
495             else:
496                 self.rest = 0
497         except Exception:
498             print self.partition
499             print self.rest
500
501     def push(self):
502         """
503         On each end of backtrack, the tail of partition has to be
504         decreased and the decrement is cancelled by an increment of
505         rest.
506         """
507         self.partition[-1] -= 1
508         self.rest += 1
509
510     def shouldterminate(self):
511         """
512         Knowing that the generation should terminate, return True to
513         stop the loop.
514         """
515         return not self.partition
516
517     def generate(self, *args):
518         """
519         Generate all partitions of integer
520
521         Variable arguments *args is used for inheritance.
522         """
523         self.set_parameters(*args)
524         if not self.ispreconditioned():
525             raise StopIteration
526
527         while True:
528             while self.rest:
529                 try:
530                     self.pull()
531                 except LookupError:
532                     break
533
534             if not self.rest:
535                 yield tuple(self.partition)
536
537             self.backtrack()
538             if self.shouldterminate():
539                 break
540             self.push()
541
542
543 class LimitedMaximumPartitionDriver(PartitionDriver):
544     """
545     Only limit the muximum of parts.
546     """
547     def __init__(self):
548         PartitionDriver.__init__(self)
549         self.limit = 0
550
551     def set_parameters(self, integer, limit):
552         """
553         Set parameters for generate:
554           integer: the number to be partitioned
555           limit: the maximum of parts
556         """
557         self.partition = []
558         self.rest = integer
559         self.limit = limit
560
561     def pull(self):
562         """
563         Pull the maximum part obtainable from rest to partition.
564         """
565         if self.partition:
566             if self.rest >= self.partition[-1]:
567                 part = self.partition[-1]
568             else:
569                 part = self.rest
570         elif self.rest >= self.limit:
571             part = self.limit
572         else:
573             part = self.rest
574
575         self.partition.append(part)
576         self.rest -= part
577
578
579 class OddPartitionDriver(PartitionDriver):
580     """
581     All parts are odd.
582     """
583     def __init__(self):
584         PartitionDriver.__init__(self)
585
586     def pull(self):
587         """
588         Pull the maximum odd part obtainable from rest to partition.
589         """
590         if self.partition and self.rest >= self.partition[-1]:
591             part = self.partition[-1]
592         elif self.rest & 1:
593             part = self.rest
594         else:
595             part = self.rest - 1
596
597         self.partition.append(part)
598         self.rest -= part
599
600     def push(self):
601         """
602         On each end of backtrack, the tail of partition has to be
603         decreased and the decrement is cancelled by an increment of
604         rest keeping oddity of parts.
605         """
606         self.partition[-1] -= 2
607         self.rest += 2
608
609
610 class OddMaximumPartitionDriver(LimitedMaximumPartitionDriver, OddPartitionDriver):
611     def __init__(self):
612         # Calling LimitedMaximum suffices, since extra work isn't needed for Odd
613         LimitedMaximumPartitionDriver.__init__(self)
614
615     def set_parameters(self, integer, limit):
616         """
617         Set parameters for generate:
618           integer: the number to be partitioned
619           limit: the maximum of parts
620         """
621         if not (limit & 1):
622             limit -= 1
623         LimitedMaximumPartitionDriver.set_parameters(self, integer, limit)
624
625     def pull(self):
626         """
627         Pull the maximum odd part obtainable from rest to partition.
628         """
629         if self.partition and self.rest >= self.partition[-1]:
630             part = self.partition[-1]
631         elif self.rest >= self.limit:
632             part = self.limit
633         elif self.rest & 1:
634             part = self.rest
635         else:
636             part = self.rest - 1
637
638         self.partition.append(part)
639         self.rest -= part
640
641
642 def partition_into_odd_generator(n, maxi=None):
643     if maxi is None:
644         driver = OddPartitionDriver()
645         return driver.generate(n)
646     else:
647         driver = OddMaximumPartitionDriver()
648         return driver.generate(n, maxi)
649
650
651 def partition_numbers_upto(n):
652     """
653     Return the partition numbers for 0 to '''n''' (inclusive).
654     """
655     penta = list(itertools.izip(
656         itertools.takewhile(lambda k: k <= n, _pentagonal()),
657         itertools.cycle((1, 1, -1, -1))))
658     p = [1]
659     for i in range(1, n + 1):
660         s = 0
662             if k > i:
663                 break
664             s += sign * p[i - k]
665         p.append(s)
666     return p
667
668 def _pentagonal():
669     """
670     Generates pentagonal and skew pentagonal numbers.
671     (1, 2, 5, 7, 12, 15, ...)
672     """
673     j = 1
674     while True:
675         yield j*(3*j - 1) >> 1
676         yield j*(3*j + 1) >> 1
677         j += 1
678
679 def partition_number(n):
680     """
681     Return the partition number for '''n'''.
682     """
683     return partition_numbers_upto(n)[-1]
684
685 def partition_conjugate(partition):
686     """
687     Return the conjugate partition of 'partition'.
688
689     For example:
690     >>> partition_conjugate((5, 3, 1))
691     (3, 2, 2, 1, 1)
692     """
693     conj = []
694     last_index = len(partition) - 1
695     for i, part in enumerate(partition):
696         if i == last_index:
697             conj.extend([i + 1] * part)
698         elif part != partition[i + 1]:
699             conj.extend([i + 1] * (part - partition[i + 1]))
700     conj.reverse()
701     return tuple(conj)
702
703
704 # aliases
705 combinationIndexGenerator = combination_index_generator
706 partitionGenerator = partition_generator
707 permutationGenerator = permutation_generator
708 DyckWordGenerator = dyck_word_generator