"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:
  422             # 1 is not found
  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
  661         for k, sign in penta:
  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