# HG changeset patch
# User Andrew Gainer-Dewar
# 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/sage/combinat/species/generating_series.py
+++ b/sage/combinat/species/generating_series.py
@@ -554,6 +554,93 @@ class CycleIndexSeries(LazyPowerSeries):
yield res
n += 1
+ def arithmetic_product(self, g):
+ """
+ Returns the arithmetic product of ``self`` with ``g``.
+
+ For species `M` and `N` such that `M[\\varnothing] = N[\\varnothing] = \\varnothing`,
+ their arithmetic product is the species `M \\boxdot N` of "`M`-assemblies of cloned `N`-structures".
+ This operation is defined and several examples are given in [MM]_.
+
+ The cycle index series for `M \\boxdot N` can be computed in terms of the component series `Z_M` and `Z_N`,
+ as implemented in this method.
+
+ EXAMPLES:
+
+ For `C` the species of (oriented) cycles and `L_{+}` the species of nonempty linear orders, `C \\boxdot L_{+}` corresponds
+ to the species of "regular octopuses"; a `(C \\boxdot L_{+})`-structure is a cycle of some length, each of whose elements
+ is an ordered list of a length which is consistent for all the lists in the structure. ::
+
+ sage: C = species.CycleSpecies().cycle_index_series()
+ sage: Lplus = species.LinearOrderSpecies(min=1).cycle_index_series()
+ sage: RegularOctopuses = C.arithmetic_product(Lplus)
+ sage: RegOctSpeciesSeq = RegularOctopuses.generating_series().counts(8)
+ sage: RegOctSpeciesSeq
+ [0, 1, 3, 8, 42, 144, 1440, 5760]
+
+ It is shown in [MM]_ that the exponential generating function for regular octopuses satisfies
+ `(C \\boxdot L_{+}) (x) = \\sum_{n \geq 1} \\sigma (n) (n - 1)! \\frac{x^{n}}{n!}` (where `\\sigma (n)` is
+ the sum of the divisors of `n`). ::
+
+ sage: RegOctDirectSeq = [0] + [sum(divisors(i))*factorial(i-1) for i in range(1,8)]
+ sage: RegOctDirectSeq == RegOctSpeciesSeq
+ True
+
+ AUTHORS:
+
+ - Andrew Gainer-Dewar (2013)
+
+ REFERENCES:
+
+ .. [MM] M. Maia and M. Mendez. "On the arithmetic product of combinatorial species." Discrete Mathematics, vol. 308, issue 23, 2008, pp. 5407-5427.
+ http://arxiv.org/abs/math/0503436.
+
+ """
+ from sage.combinat.partition import Partition, Partitions
+ from sage.combinat.species.stream import Stream, _integers_from
+ from sage.rings.arith import gcd, lcm, divisors
+ from itertools import product, repeat, chain
+
+ p = self.base_ring()
+
+ assert self.coefficient(0) == p.zero()
+ assert g.coefficient(0) == p.zero()
+
+ # We first define an operation `\\boxtimes` on partitions as in Lemma 2.1 of [MM]_.
+ def arith_prod_of_partitions (l1, l2):
+ # Given two partitions `l_1` and `l_2`, we construct a new partition `l_1 \\boxtimes l_2` by
+ # the following procedure: each pair of parts `a \\in l_1` and `b \\in l_2` contributes a part
+ # `\\lcm (a, b)` to `l_1 \\boxtimes l_2` with multiplicity `\\gcm (l_1, l_2)``.
+ term_iterable = chain.from_iterable( repeat(lcm(pair), times=gcd(pair)) for pair in product(l1, l2) )
+ term_list = sorted(term_iterable, reverse=True)
+ res = Partition(term_list)
+ return res
+
+ # We then extend this to an operation on symmetric functions as per eq. 52 of [MM]_.
+ def arith_prod_sf (x, y):
+ ap_sf_wrapper = lambda l1, l2: p(arith_prod_of_partitions(l1, l2))
+ return p._apply_multi_module_morphism(x, y, ap_sf_wrapper)
+
+ # Sage stores cycle index series by degree.
+ # Thus, to compute the arithmetic product `Z_M \\boxdot Z_N` it is useful
+ # to compute all terms of a given degree `n` at once.
+ def arith_prod_coeff (n):
+ if n == 0:
+ res = p.zero()
+ else:
+ index_set = ((d, n/d) for d in divisors(n))
+ res = sum(arith_prod_sf(self.coefficient(i), g.coefficient(j)) for i,j in index_set)
+
+ # Build a list which has res in the nth slot and 0's before and after
+ # to feed to sum_generator
+ res_in_seq = [p.zero()]*n + [res, p.zero()]
+
+ return self.parent(res_in_seq)
+
+ # Finally, we use the sum_generator method to assemble these results into a single
+ # LazyPowerSeries object.
+ return self.parent().sum_generator(arith_prod_coeff(n) for n in _integers_from(0))
+
def _cycle_type(self, s):
"""
EXAMPLES::