Ticket #14542: trac_14542_cycle_index_arithmetic_product.patch

File trac_14542_cycle_index_arithmetic_product.patch, 4.9 KB (added by agd, 9 years ago)

Patch to implement arithmetic product of cycle indices

  • sage/combinat/species/generating_series.py

    # HG changeset patch
    # User Andrew Gainer-Dewar <andrew dot gainer dot dewar at gmail.com>
    # Date 1367877237 18000
    # Node ID bd8b4321be5a6615e6e5680449af6d75fe8a5868
    # Parent  9c903bae39b955af5d4e5fa1a33bdc9e4c21b7bd
    Implement arithmetic product of cycle index series
    
    diff --git a/sage/combinat/species/generating_series.py b/sage/combinat/species/generating_series.py
    a b class CycleIndexSeries(LazyPowerSeries): 
    554554            yield res
    555555            n += 1
    556556
     557    def arithmetic_product(self, g):
     558        """
     559        Returns the arithmetic product of ``self`` with ``g``.
     560       
     561        For species `M` and `N` such that `M[\\varnothing] = N[\\varnothing] = \\varnothing`,
     562        their arithmetic product is the species `M \\boxdot N` of "`M`-assemblies of cloned `N`-structures".
     563        This operation is defined and several examples are given in [MM]_.
     564
     565        The cycle index series for `M \\boxdot N` can be computed in terms of the component series `Z_M` and `Z_N`,
     566        as implemented in this method.
     567       
     568        EXAMPLES:
     569
     570        For `C` the species of (oriented) cycles and `L_{+}` the species of nonempty linear orders, `C \\boxdot L_{+}` corresponds
     571        to the species of "regular octopuses"; a `(C \\boxdot L_{+})`-structure is a cycle of some length, each of whose elements
     572        is an ordered list of a length which is consistent for all the lists in the structure. ::
     573
     574            sage: C = species.CycleSpecies().cycle_index_series()
     575            sage: Lplus = species.LinearOrderSpecies(min=1).cycle_index_series()
     576            sage: RegularOctopuses = C.arithmetic_product(Lplus)
     577            sage: RegOctSpeciesSeq = RegularOctopuses.generating_series().counts(8)
     578            sage: RegOctSpeciesSeq
     579            [0, 1, 3, 8, 42, 144, 1440, 5760]
     580
     581        It is shown in [MM]_ that the exponential generating function for regular octopuses satisfies
     582        `(C \\boxdot L_{+}) (x) = \\sum_{n \geq 1} \\sigma (n) (n - 1)! \\frac{x^{n}}{n!}` (where `\\sigma (n)` is
     583        the sum of the divisors of `n`). ::
     584
     585            sage: RegOctDirectSeq = [0] + [sum(divisors(i))*factorial(i-1) for i in range(1,8)]
     586            sage: RegOctDirectSeq == RegOctSpeciesSeq
     587            True           
     588       
     589        AUTHORS:
     590
     591        - Andrew Gainer-Dewar (2013)
     592
     593        REFERENCES:
     594       
     595        .. [MM] M. Maia and M. Mendez. "On the arithmetic product of combinatorial species." Discrete Mathematics, vol. 308, issue 23, 2008, pp. 5407-5427.
     596          http://arxiv.org/abs/math/0503436.
     597       
     598        """
     599        from sage.combinat.partition import Partition, Partitions
     600        from sage.combinat.species.stream import Stream, _integers_from
     601        from sage.rings.arith import gcd, lcm, divisors
     602        from itertools import product, repeat, chain
     603
     604        p = self.base_ring()
     605
     606        assert self.coefficient(0) == p.zero()
     607        assert g.coefficient(0) == p.zero()
     608
     609        # We first define an operation `\\boxtimes` on partitions as in Lemma 2.1 of [MM]_.
     610        def arith_prod_of_partitions (l1, l2):
     611            # Given two partitions `l_1` and `l_2`, we construct a new partition `l_1 \\boxtimes l_2` by
     612            # the following procedure: each pair of parts `a \\in l_1` and `b \\in l_2` contributes a part
     613            # `\\lcm (a, b)` to `l_1 \\boxtimes l_2` with multiplicity `\\gcm (l_1, l_2)``.
     614            term_iterable = chain.from_iterable( repeat(lcm(pair), times=gcd(pair)) for pair in product(l1, l2) )
     615            term_list = sorted(term_iterable, reverse=True)
     616            res = Partition(term_list)
     617            return res
     618
     619        # We then extend this to an operation on symmetric functions as per eq. 52 of [MM]_.
     620        def arith_prod_sf (x, y):
     621            ap_sf_wrapper = lambda l1, l2: p(arith_prod_of_partitions(l1, l2))
     622            return p._apply_multi_module_morphism(x, y, ap_sf_wrapper)
     623
     624        # Sage stores cycle index series by degree.
     625        # Thus, to compute the arithmetic product `Z_M \\boxdot Z_N` it is useful
     626        # to compute all terms of a given degree `n` at once.
     627        def arith_prod_coeff (n):
     628            if n == 0:
     629                res = p.zero()
     630            else:
     631                index_set = ((d, n/d) for d in divisors(n))
     632                res = sum(arith_prod_sf(self.coefficient(i), g.coefficient(j)) for i,j in index_set)
     633
     634            # Build a list which has res in the nth slot and 0's before and after
     635            # to feed to sum_generator
     636            res_in_seq = [p.zero()]*n + [res, p.zero()]
     637           
     638            return self.parent(res_in_seq)
     639
     640        # Finally, we use the sum_generator method to assemble these results into a single
     641        # LazyPowerSeries object.
     642        return self.parent().sum_generator(arith_prod_coeff(n) for n in _integers_from(0))
     643
    557644    def _cycle_type(self, s):
    558645        """
    559646        EXAMPLES::