\\ PARI Functions for Combinatorial Species
\\ ----------------------------------------
\\ Author: Andrew Howroyd, Dec 2020
\\ This script or parts thereof may be freely used for private, educational, open source or commercial purposes.
\\ No waranty of any kind. Use at own risk.

\\ This is not intended to be a library.
\\     Source code will not be updated with improvements and bug fixes. 
\\     Functions do not check arguments for correct usage. Garbage in garbage out.
\\
\\ If you use any of the functions to enumerate new sequences, you should independently check initial terms by
\\ another method.  You are 100% responsible. Do not blame the script or author for wrong results / bugs.
\\
\\ Please feel free to adapt any lines for your own personal collection. There is no need to acknowledge
\\ source when using. (Once you copy characters from this file they are basically yours)



\\ How to use:
\\ Just cut and paste everything from the top down to the start of examples into a (preferably fresh) instance
\\ of PARI-GP. All non-executable text is in the form of comments.
\\ Also, give your PARI environment more memory by calling 'allocatemem' a few times.


\\ References:
\\ The methods presented here are for practical enumeration purposes.
\\ In comments ITSS means Bergeron/Labelle/Leroux reference.
\\ Wikipedia:
\\   https://en.wikipedia.org/wiki/Combinatorial_species
\\   https://en.wikipedia.org/wiki/Cycle_index
\\   https://en.wikipedia.org/wiki/P%C3%B3lya_enumeration_theorem
\\   https://en.wikipedia.org/wiki/Burnside%27s_lemma
\\ References of a more mathematical nature: (available online)
\\   F. Bergeron, G. Labelle, and P. Leroux, "Introduction to the Theory of Species of Structures", Nov 2013.
\\   F. Harary and E. M. Palmer, Graphical Enumeration, Academic Press, NY, 1973.
\\ There is also a published book by the first authors:
\\   F. Bergeron, G. Labelle and P. Leroux, Combinatorial Species and Tree-Like Structures, Camb. 1998.
\\ Other background reading on combinatorial classes: (available online)
\\   P. Flajolet and R. Sedgewick, Analytic Combinatorics, 2009.
\\ 
\\ See also PARI functions for Polya Enumeration (https://oeis.org/A056391/a056391.txt)
\\ This is in part derived from that, but is geared to handling power series of cycle indexes, which
\\ gives many new capabilitues. (In particular the ability to determine the cycle index of objects
\\ that have a complex automorphism group such as trees and connected graphs).

\\ Representation in PARI
\\ Combinatorial species are too abstract to represent directly in PARI.
\\ Instead the cycle index series is used, which is technically only a property of a combinatorial species.
\\ For example, the cycle index probably contains too little information for 'functorial composition' (2.4 in ITSS).
\\ (In this script the term 'species' is used to mean the cycle index series of a species)
\\ Roughly speaking a cycle index series is a power series (in the variable x here) of cycle indexes.
\\ A cycle index is a multivariate polynomial (in the variables s1, s2, ... here).
\\ The sum of the exponents of s1, s2 ... in any monomial must not exceed the exponent of x.
\\ In a sense, a cycle index series is a power series with coefficients over all partitions of the exponent.
\\ A typical cycle index series looks like: (this one for the symmetric groups up to 3 terms)
\\   1 + s1*x + (1/2*s1^2 + 1/2*s2)*x^2 + (1/6*s1^3 + 1/2*s2*s1 + 1/3*s3)*x^3 + O(x^4)
\\ The initial constant term 1 is not required. Some operations are easier with and some without.
\\ A typical series can use a lot of memory!
\\ Most methods also support bivaritate series in x and y. This allows for weighted enumeration.
\\ Arbitrary multivariate series will not work. 
\\ A cycle index series represents a 'combinatorial class' [Flajolet and Sedgewick], like sets, trees, graphs, etc.
\\ The cycle index carries enough information for both labeled and unlabeled enumeration.
\\

\\ Basic operations:
\\ Important basic operation are handled natively (as power series)
\\ A + B: is the disjoint union of two combinatorial classes A and B.
\\ A - B: removes the objects enumerated by B from A.
\\ c*A: with c a constant forms the class of c copies of A.
\\ A * B: is the composition of two combinatorial classes A and B. The composition is essentially
\\   the class consisting of ordered pairs (a,b) with a in A and b in B.
\\ A / B: division is the inverse operation to multiplication. An example use is:
\\   Gcr = Gr / G meaning there is a bijection between a rooted graph (Gr) and a pair consisting 
\\   of a connected rooted graph (Gcr) and an unrooted graph (G). Although division often appears
\\   in program scripts and formula, it is usually easier to comprehend its meaning in terms of multiplication. 
\\ 1/(1 - A): is the class of sequences of zero or more objects of A. It is assumed that A does not
\\   have any object of size 0. (Any constant term should first be removed).
\\ x*sv(1): is the class of a singleton. In most texts this is denoted Z or E. Most of the methods require
\\   a power series representation so in practice you must use x*sv(1) + O(x^n) where n is at least 2.

\\ Other operations are implemented below in PARI.
\\ A(B): Substitution implemented by sSubstOp(A,B); sSolve is inverse.
\\ A X B: Cartesian product.
\\ ^A: Pointing.
\\ Others are exponential and log transforms. The exponential transform is just a short hand
\\ for a substitution with the first argument being sets (symmetric group series) - in other words it
\\ gives the combinatorial class of sets of the second argument.


\\ PARI performance problem (bug) and work around:
\\ Versions of PARI up to at least 2.13.0 contain a critical performance problem in the handling of multivariate
\\ series. See https://pari.math.u-bordeaux.fr/cgi-bin/bugreport.cgi?bug=2264 for status.
\\ Fortunately, there is a good work around (discovered accidentally by author), albeit slightly messy.
\\ The work around which is mostly encapsulated in 'xsubstvec' (a replacement for substvec) is to add a term with
\\ a floating point coeficient at the end of the series.

\\  ************  START OF CODE ************


\\ Hack for bug in PARI - approximate replacement for substvec.
\\ Assumes s is a power series.
\\ Returned series will need to be truncated, left to caller.
xsubstvec(s, v, w)={my(n=serprec(s,x)); substvec(0.0*x^n + O(x*x^n) + Pol(s), v, w)}

\\ Gets n-th subscripted variable.
\\ sv(1) => s1, sv(2) => s2, etc.
sv(n)={eval(Str("'s",n))}

\\ Gets the monomial (in s1, s2, ...) corresponding to a partition.
\\ sMonomial([1,1,3,5]) => s1^2*s3*s5.
sMonomial(v)={prod(k=1, #v, sv(v[k]))}

\\ Cycle index series for identity group. Species of linear objects.
\\ linearSeries(3) => 1 + s1*x + s1^2*x^2 + s1^3*x^3 + O(x^4).
linearSeries(n)={1/(1-x*sv(1)) + O(x*x^n)}

\\ Cycle index of cyclic group of order n. Z(C_n).
\\ cycGroupCycleIndex(6) => 1/6*s1^6 + (1/6*s2^3 + (1/3*s3^2 + 1/3*s6)).
cycGroupCycleIndex(n)={sumdiv(n, d, eulerphi(d)*sv(d)^(n/d))/n}

\\ Cycle index series of the cyclic groups up to n. Includes inital 1.
\\ cycGroupSeries(2) => 1 + s1*x + (1/2*s1^2 + 1/2*s2)*x^2 + O(x^3).
cycGroupSeries(n)={1 + sum(k=1, n, cycGroupCycleIndex(k)*x^k) + O(x*x^n)}

\\ Cycle index of dihedral group of order n. Z(D_n).
dihedralGroupCycleIndex(n)={(cycGroupCycleIndex(n) + if(n%2, sv(1)*sv(2)^(n\2), sv(2)^(n/2-1)*(sv(1)^2 + sv(2))/2))/2}
\\ Cycle index series of dihedral groups up to n.
dihedralGroupSeries(n)={1 + sum(k=1, n, dihedralGroupCycleIndex(k)*x^k) + O(x*x^n)}

\\ Cycle index series of symmetric groups up to n. Species of sets.
symGroupSeries(n)={my(v=vector(1+n)); v[1] = 1; for(n=1, #v-1, v[n+1] = sum(k=1, n, sv(k)*v[1+n-k])/n); Ser(v)}
\\ Cycle index of symmetric group of order n.
symGroupCycleIndex(n)={polcoef(symGroupSeries(n),n)}

\\ Cycle index series for strongly normal labelings.
stronglyNormalSeries(n)={my(p=symGroupSeries(n)); 1/prod(k=1, n, 1 - polcoef(p,k)*x^k + O(x*x^n))}

\\ Operations
\\ In power series p, substitutes x with x^e, y with y^e and s_i with s_{i*e}.
\\ sRaise(4*s1*y^2*x + s2*s1^2*x^4 + O(x^5), 3) => 4*s3*y^6*x^3 + s6*s3^2*x^12 + O(x^15)
sRaise(p,e)={my(n=serprec(p,x)-1); O(x^(n*e+e)) + xsubstvec(p, concat([x, 'y], vector(n, i, sv(i))), concat([x^e, 'y^e], vector(n, i, sv(i*e))))}

\\ In a cycle index, substitutes s_i with s_{i*e} for i in 1..n.
sRaiseCI(p, n, e)={polcoef(xsubstvec(p+O(x), vector(n,i,sv(i)), vector(n,i,sv(i*e))), 0)}

\\ In power series p, substitutes s_i with f(i)
sFuncSubst(p, f)={my(n=serprec(p,x)-1); O(x*x^n) + xsubstvec(p, vector(n,i,sv(i)), vector(n,i,f(i)))}

\\ apparantly called plethystic substitution
sApplyCI(p, pn, q, qn)={my(qv=vector(qn, k, sv(k))); substvec(p, vector(pn, i, sv(i)), vector(pn,i, substvec(q, qv, vector(qn, k, sv(i*k)))))}



\\ Species substitution binary operation.
\\ Returns species A(B).
\\ See 2.2.1 in ITSS.
\\ Effectively replaces the unit objects of A with objects of B.
\\ This method allows A or B to be polynomials instead of power series. If both are polynomials then the result is a polynomial.
\\ Otherwise a power series is returned. In many cases, the returned series precision is currently suboptimal.
sSubstOp(A, B)={
  B -= polcoef(B,0);
  my(m=serprec(A,x), n=serprec(B,x), ispol=0);
  if(m==oo,
     m=poldegree(A,x); if(n==oo, n=m*poldegree(B); ispol=1, n--),
     m--; if(n==oo, n=m, n=min(m, n-1)));
  if(m > n, m=n);
  my(vars=vector(n, i, sv(i)));
  my(f=vector(m, k, xsubstvec(B + O(x*x^(n\k)), concat([x, 'y], vars[1..n\k]), concat([x^k, 'y^k], vector(n\k, i, vars[i*k])))/x^k ));
  \\ my(r=O(x*x^n) + xsubstvec(A + O(x*x^n), vars[1..m], f));
  \\ Faster version of above (in v2.11.1 of PARI) - diy outer loop
  my(r=polcoef(A,0) + sum(i=1, m, my(t=polcoef(A,i,x), X=O(x^(n-i+1)));
     x^i*substvec(t + X, vars[1..i], vector(i,j,f[j] + X))
  ));
  if(ispol, Pol(r), r);
}

\\ Inverse of sSubstOp.
\\ Returns species A such that Z=A(B).
sSolve(Z, B)={
  B -= polcoef(B,0);
  my(m=serprec(Z,x)-1, n=serprec(B,x));
  if(n==oo, n=m, n=min(m, n-1));
  if(m > n, m=n);
  my(vars=vector(n, i, sv(i)));
  my(f=vector(m, k, xsubstvec(B + O(x*x^(n\k)), concat([x, 'y], vars[1..n\k]), concat([x^k, 'y^k], vector(n\k, i, vars[i*k])))/x^k ));
  my(b1=polcoef(polcoef(B,1,x), 1, sv(1)));
  my(r=polcoef(Z,0,x));
  for(i=1, n, my(t=polcoef(Z,i,x)/b1); r += t*x^i;
     \\      Z-=x^i*xsubstvec(t + O(x^(n-i+1)), vars[1..i], f[1..i])
     \\ Faster version of above (in v2.11.1 of PARI)
     my(X=O(x^(n-i+1))); Z-=x^i*substvec(t + X, vars[1..i], vector(i,j,f[j] + X))
  );
  r + O(x*x^n);
}

\\ Cycle index series reversion.
\\ Returns species A such that X = A(B), where X is the unit species.
sRevert(B)={
  sSolve(x*sv(1) + O(x^serprec(B,x)), B)
}

\\ Exponential transform. Gives species of sets of A.
\\ Shorthand for sSubstOp with first argument symmetric series
sExp(A)={
  my(n=serprec(A,x)-1, vars=concat([x, 'y], vector(n,i,sv(i))));
  A -= polcoef(A, 0);
  exp( O(x*x^n) + sum(k=1, n, xsubstvec(A + O(x*x^(n\k)), vars[1..2+n\k], concat([x^k, 'y^k], vector(n\k, i, vars[2+i*k])))/k ))
}

\\ Logarithm transform (inverse of exponential).
sLog(A)={
  my(n=serprec(A,x)-1, vars=concat([x, 'y], vector(n,i,sv(i))));
  A = log(1 - polcoef(A,0) + A);
  O(x*x^n) + sum(k=1, n, moebius(k)*xsubstvec(A + O(x*x^(n\k)), vars[1..2+n\k], concat([x^k, 'y^k], vector(n\k, i, vars[2+i*k])))/k )
}

\\ In a previous version, sExp/sLog were called sEulerT/sInvEulerT
sEulerT=sExp;
sInvEulerT=sLog;

\\ Weigh transform: distinct collections (set) of elements of a combinatorial class.
\\ Warning - this is probably not a legitimate species operation. Use with extreme caution and expect wrong results!
sWeighT(p)={
  my(n=serprec(p,x)-1, vars=concat([x, 'y], vector(n,i,sv(i)))); 
  p -= polcoef(p, 0);
  exp( O(x*x^n) + sum(k=1, n, (-1)^(k-1)*xsubstvec(p + O(x*x^(n\k)), vars[1..2+n\k], concat([x^k, 'y^k], vector(n\k, i, vars[2+i*k])))/k ))
}

\\ Inverse Weigh transform.
sInvWeighT(p)={
  my(n=serprec(p,x)-1, vars=concat([x, 'y], vector(n,i,sv(i))));
  p = log(1 - polcoef(p,0) + p);
  O(x*x^n) + sum(k=1, n, my(b=valuation(k,2)); moebius(k/2^b)*2^max(0,b-1)*xsubstvec(p + O(x*x^(n\k)), vars[1..2+n\k], concat([x^k, 'y^k], vector(n\k, i, vars[2+i*k])))/k )
}


\\ Pointed combinatorial class (rooted objects from non-rooted objects)
\\ Return species ^A
\\ See 2.3.1 in ITSS.
sPoint(A)={my(s1=sv(1)); s1*deriv(A, s1)}

\\ Cartesian product - often called superposition.
\\ Returns species A X B.
\\ See 2.3.2 in ITSS.
sCartProd(A,B)={
  my(n=min(serprec(A,x),serprec(B,x))-1);
  my(vars=vector(n, i, sv(i)));
  my(recurse(k, r, p, q)=if(p && q, if(k>r, p*q, my(v=vars[k], b=x^k); sum(i=0, min(poldegree(p,v), poldegree(q,v)), (b*v*k)^i*i!*self()(k+1, r-i*k, polcoef(p,i,v)/b^i, polcoef(q,i,v)/b^i)))));
  recurse(1, n, Pol(A), Pol(B)) + O(x*x^n);
}

\\ Cartesian product of A with itself e times.
\\ Equivalent to calling sCartProd with itself e times.
sCartPower(A, e)={
  my(n=serprec(A,x)-1);
  my(vars=vector(n, i, sv(i)));
  my(recurse(k, r, p)=if(p, if(k>r, p^e, my(v=vars[k], b=x^k); sum(i=0, poldegree(p,v), (b*v*k^(e-1))^i*(i!)^(e-1)*self()(k+1, r-i*k, polcoef(p,i,v)/b^i)))));
  recurse(1, n, Pol(A)) + O(x*x^n);  
}

\\ Cycle index from a permutation
\\ PermCycleIndex([2,3,1,4]) => s_3*s_1.
PermCycleIndex(p)={my(r=1);for(i=1,#p,my(m=1,t=p[i]);while(t>i,t=p[t];m++);if(t==i,r*=sv(m)));r}
\\ Cycle index from a group of permutations
GroupCycleIndex(group)={sum(i=1,#group,PermCycleIndex(group[i]))/#group}

\\ Converts a partition into a polynomial in x. Each term t in the input list is made into x^t.
\\ PartPoly([1,1,1,5],x) => x^5 + 3*x
PartPoly(p,x)={sum(i=1, #p, x^p[i])}

\\ The number of permutations corresponding to a given partition. The partition is given in polynomial form.
\\ PermMultiplicityByPart(PartPoly([1,1,1,5], x)) => 1344
PermMultiplicityByPart(partPoly) = {my(m=1,s=0); for(i=0, poldegree(partPoly), my(c=polcoef(partPoly,i)); if(c, s+=c*i; m*=c!*i^c)); s!/m}

\\ Same as PermMultiplicity(PartPoly(p))
\\ Included here as a helper (used in many graph sequences).
permcount(v) = {my(m=1, s=0, k=0, t); for(i=1, #v, t=v[i]; k=if(i>1&&t==v[i-1], k+1, 1); m*=t*k; s+=t); s!/m}


QPoly(n)={sumdiv(n,d,d*sv(d))}


\\ Computes the number of non equivalent colorings where all color permutations are allowed. (color group is symmetric group)
StructsByCycleIndex(n, cycleIndex, colors)={
    my(sv = vector(n, i, sv(i)));
    my(qp = vector(n, i, QPoly(i)));
    my(total = 0);
    forpart(p=colors, my(px=PartPoly(p,x)); total += PermMultiplicityByPart(px) * substvec(cycleIndex, sv, substvec(qp, sv, vector(n, i, polcoef(px,i)) ) ) );
    total / colors!;
}

\\ Gets the generating series for unlabeled objects.
\\ Replace s_i in the cycle index series with 1: Z[species; 1, 1, 1, ... ]
\\ If 'k' is given as the parameter then that is used in place of 1.
OgfSeries(species, k=1)={sFuncSubst(species, i->k)}

\\ Gets the exponential generating series for labeled objects.
EgfSeries(species)={sFuncSubst(species, i->i==1)}

\\ Numbers of labeled object as a vector.
NumLabeledObjsSeq(species)={Vec(serlaplace(EgfSeries(species)))}

\\ Number of unlabeled objects as a vector.
NumUnlabeledObjsSeq(species)={Vec(OgfSeries(species))}

\\ Number of colorings of the unlabeled objects using at most k colors (default k=2).
ColoringsSeq(species, k=2)={Vec(OgfSeries(species, k))}

\\ Number of colorings of the unlabeled objects using exactly k colors (default k=2).
ColoringsExactlySeq(species, k=2)={my(y=sv(1)); p=OgfSeries(species, y); Vec(sum(j=0, k, (-1)^(k-j)*subst(p, y, j)*binomial(k,j)), 1-serprec(species,x))}

\\ Number of colorings of the unlabeled objects by number of colors used. Returns matrix.
\\ The first column will be number of unlabeled objects, and main diagonal number of labeled objects.
ColoringsTriangle(species)={my(y=sv(1), n=serprec(species,x)-1, p=OgfSeries(species, y)-polcoef(species,0), v=vector(n, k, Vec(subst(p, y, k))~)); Mat(vector(n, k, sum(i=1, k, (-1)^(k-i)*binomial(k, i)*v[i])))}

\\ Number of labelings of the unlabeled objects covering an initial interval of positive integers. Row sums of ColoringsTriangle.
NormalLabelingsSeq(species)={my(y=sv(1), n=serprec(species,x)-1, p=OgfSeries(species, y)); Vec(sum(k=0, n, subst(p,y,k)*sum(r=k, n, binomial(r, k)*(-1)^(r-k)))) }

\\ Number of strongly normal labelings of the objects.
\\ Strongly normal means labels have weakly decreasing frequency.
StronglyNormalLabelingsSeq(species)={NumUnlabeledObjsSeq(sCartProd(species, stronglyNormalSeries(serprec(species,x)-1)))}

\\ Number of inequivalent colorings of the unlabeled objects (up to permutation of colors)
InequivalentColoringsSeq(species, colors=oo)={
   my( n=serprec(species,x)-1,
       sv = vector(n, i, sv(i)),
       qp = sum(i=1, n, QPoly(i)*x^i) + O(x*x^n));
    my(total = O(x*x^n));
    colors=min(colors,n);
    forpart(p=colors, my(px=PartPoly(p,x)); total += PermMultiplicityByPart(px) * xsubstvec(species, sv, Vec(O(x*x^n)+xsubstvec(qp, sv, vector(n, i, polcoef(px,i)) ), -n) ) );
    Vec(total / colors!)
}

\\ Number of inequivalent colorings of the unlabeled objects using exactly k colors.
InequivalentColoringsExactlySeq(species, colors=2)={InequivalentColoringsSeq(species, colors) - if(colors>1, InequivalentColoringsSeq(species, colors-1))}

\\ Number of inequivalent coloring of the unlabeled objects by number of colors used. Returns matrix.
\\ The first column and the main diagonal will be the number of unlabeled objects.
InequivalentColoringsTriangle(species)={my(v=vector(serprec(species,x)-1, k, InequivalentColoringsSeq(species,k)~)); Mat(vector(#v, k, v[k]-if(k>1,v[k-1])) )}


\\ ********************************************************************************************
\\ **********************                END OF CODE                  *************************
\\ ********************************************************************************************
\\ The rest of the file is examples.

\\ Examples of sequences that have been extended using the above functions.
\\ These mostly relate to unlabeled graphs and inequivalent colorings and non-isomorphic objects.

\\ A002218: Number of unlabeled nonseparable (or 2-connected) graphs (or blocks) with n nodes.
\\ A004115: Number of unlabeled rooted nonseparable graphs with n nodes.
\\ A007145: Number of rooted bridgeless graphs with n nodes.
\\ A101389: Number of n-vertex unlabeled oriented graphs without endpoints.
\\ A259115: Number of unrooted binary ordered tanglegrams of size n.
\\ A318226: Number of inequivalent leaf-colorings of rooted trees with n nodes.
\\ A318228: Number of inequivalent leaf-colorings of planted achiral trees with n nodes.
\\ A318230: Number of inequivalent leaf-colorings of binary rooted trees with 2n + 1 nodes.
\\ A318231: Number of inequivalent leaf-colorings of series-reduced rooted trees with n nodes.
\\ A320173: Number of inequivalent colorings of series-reduced balanced rooted trees with n leaves.
\\ A330465: Number of non-isomorphic series-reduced rooted trees whose leaves are multisets with a total of n elements.
\\ A330467: Number of series-reduced rooted trees whose leaves are multisets whose multiset union is a strongly normal multiset of size n.
\\ A330470: Number of non-isomorphic series/singleton-reduced rooted trees on a multiset of size n.
\\ A339233: Number of inequivalent colorings of oriented series-parallel networks with n colored elements.
\\ A339287: Number of inequivalent colorings of unoriented series-parallel networks with n colored elements.
\\ A339645: T(n,k) is the number of inequivalent colorings of series-reduced rooted trees with n leaves of exactly k colors.


\\ ============================================================
\\ Example: Species of sets
symGroupSeries(3)
\\ => 1 + s1*x + (1/2*s1^2 + 1/2*s2)*x^2 + (1/6*s1^3 + 1/2*s2*s1 + 1/3*s3)*x^3 + O(x^4)

\\ Number of labeled sets of size n.
NumLabeledObjsSeq(symGroupSeries(10))
\\ => [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]

\\ Number of unlabeled sets of size n.
NumUnlabeledObjsSeq(symGroupSeries(10))
\\ => [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]

\\ Number of inequivalent ways to color a set of size n.
InequivalentColoringsSeq(symGroupSeries(10))
\\ => [1, 1, 2, 3, 5, 7, 11, 15, 22, 30, 42] = A000041 (number of partitions)

\\ Number of inequivalent colorings of an n-set by number of colors used.
InequivalentColoringsTriangle(symGroupSeries(10))
\\ => A008284 (number of partitions of n into k positive parts)

\\ Unlabeled sets of lists are integer partitions
NumUnlabeledObjsSeq(sExp(linearSeries(12)))
\\ => [1, 1, 2, 3, 5, 7, 11, 15, 22, 30, 42, 56, 77] = A000041

\\ The labeled equivalant is:
NumLabeledObjsSeq(sExp(linearSeries(10)))
\\ => [1, 1, 3, 13, 73, 501, 4051, 37633, 394353, 4596553, 58941091] = A000262

\\ Unlabeled sets of sets also give integer partitions
NumUnlabeledObjsSeq(sExp(symGroupSeries(12)))
\\ => [1, 1, 2, 3, 5, 7, 11, 15, 22, 30, 42, 56, 77] = A000041

\\ The labeled equivalent is:
NumLabeledObjsSeq(sExp(symGroupSeries(10)))
\\ => [1, 1, 2, 5, 15, 52, 203, 877, 4140, 21147, 115975] = A000110

\\ Sets of sets of lists are partitions of partitions
NumUnlabeledObjsSeq(sExp(sExp(linearSeries(12))))
\\ => [1, 1, 3, 6, 14, 27, 58, 111, 223, 424, 817, 1527, 2870] = A001970

\\ List of sets
NumLabeledObjsSeq(1/(2-symGroupSeries(10)))
\\ => [1, 1, 3, 13, 75, 541, 4683, 47293, 545835, 7087261, 102247563] = A000670 (Fubini numbers)

\\ Sets of distinct sets are partitions into distinct parts
NumUnlabeledObjsSeq(sWeighT(linearSeries(15)))
\\ => [1, 1, 1, 2, 2, 3, 4, 5, 6, 8, 10, 12, 15, 18, 22, 27] = A000009

\\ A permutation can be considered to be a collection of labeled cycles
NumLabeledObjsSeq(sExp(cycGroupSeries(10)))
\\ => [1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800] = A000142

\\ Removing cycles of length 1 gives derangements
NumLabeledObjsSeq(sExp(cycGroupSeries(10) - x*sv(1)))
\\ => [1, 0, 1, 2, 9, 44, 265, 1854, 14833, 133496, 1334961] = A000166

\\ Number of unlabeled 2-regular graphs with loops forbidden
NumUnlabeledObjsSeq(sExp(cycGroupSeries(12) - x*sv(1)))
\\ => [1, 0, 1, 1, 2, 2, 4, 4, 7, 8, 12, 14, 21] = A002865

\\ Set partitions up to rotation and reflection
InequivalentColoringsSeq(dihedralGroupSeries(12)-1)
\\ => [1, 2, 3, 7, 12, 37, 93, 354, 1350, 6351, 31950, 179307] = A084708

\\ Inequivalent colorings of bracelets by number of colors used
InequivalentColoringsTriangle(dihedralGroupSeries(10)-1)
\\ => A152176

\\ The Cartesian product implements superposition
\\ Inequivalent colorings are a superposition with a set of sets (each set is a color)
NumUnlabeledObjsSeq(sCartProd(dihedralGroupSeries(12)-1, sExp(symGroupSeries(12))))
\\ => [1, 2, 3, 7, 12, 37, 93, 354, 1350, 6351, 31950, 179307] = A084708

\\ The superposition of 3 matchings
my(M=sExp(x^2*symGroupCycleIndex(2) + O(x^17))); NumUnlabeledObjsSeq(sCartProd(M,sCartProd(M,M)))
\\ => [1, 0, 1, 0, 5, 0, 16, 0, 86, 0, 448, 0, 3580, 0, 34981, 0, 448628] = aerated A002830

\\ The superposition of 3 matchings written another way
NumUnlabeledObjsSeq(sCartPower(sExp(x^2*symGroupCycleIndex(2) + O(x^17)), 3))
\\ => [1, 0, 1, 0, 5, 0, 16, 0, 86, 0, 448, 0, 3580, 0, 34981, 0, 448628] = aerated A002830

\\ Achiral patterns
my(N=12); InequivalentColoringsSeq(sum(k=0, N, sv(1)^(k%2)*sv(2)^(k\2)*x^k) + O(x*x^N))
\\ => [1, 1, 2, 3, 7, 12, 31, 59, 164, 339, 999, 2210, 6841] = A080107

\\ Unoriented patterns
my(N=12); InequivalentColoringsSeq(sum(k=0, N, (sv(1)^k + sv(1)^(k%2)*sv(2)^(k\2))*x^k/2) + O(x*x^N))
\\ => [1, 1, 2, 4, 11, 32, 117, 468, 2152, 10743, 58487, 340390, 2110219] = A103293


\\ ============================================================
\\ Examples on species of trees.
\\ The species of rooted tree can be constructed iteratively.
rootedTreesSeries(n)={my(p = sv(1)*x + O(x^2)); for(n=2, n, p = sv(1)*x*sExp(p)); p}
rootedTreesSeries(4)
\\ s1*x + s1^2*x^2 + (3/2*s1^3 + 1/2*s2*s1)*x^3 + (8/3*s1^4 + s2*s1^2 + 1/3*s3*s1)*x^4 + O(x^5)

\\ Number of labeled rooted trees
NumLabeledObjsSeq(rootedTreesSeries(10))
\\ => [1, 2, 9, 64, 625, 7776, 117649, 2097152, 43046721, 1000000000] = A000169

\\ Number of unlabeled rooted trees
NumUnlabeledObjsSeq(rootedTreesSeries(10))
\\ => [1, 1, 2, 4, 9, 20, 48, 115, 286, 719] = A000081

\\ The following is probably the most important example to understand (everything)
\\ A rooted tree can also be defined by the equation: T = x*sv(1)*sExp(T), or rewitten x*sv(1) = T/sExp(T).
\\ In labeled terms serreverse(x/exp(x)) gives e.g.f. of A000169. This is well known.
\\ This doesn't translate to unlabeled o.g.f because EulerT isn't a function (so the normal method is iteration as shown above).
\\ But, with series of cycle indices the property is restored.
rootedTreesSeriesAlt(n)={sRevert(x*sv(1)/sExp(x*sv(1) + O(x^n)))}
NumUnlabeledObjsSeq(rootedTreesSeriesAlt(10))
\\ => [1, 1, 2, 4, 9, 20, 48, 115, 286, 719] = A000081 (Amazing!)

\\ Pairs of rooted trees
NumUnlabeledObjsSeq(rootedTreesSeries(10)^2)
\\ => [1, 2, 5, 12, 30, 74, 188, 478, 1235, 3214] = A000106

\\ Unordered pairs of rooted trees
NumUnlabeledObjsSeq(sSubstOp(x^2*symGroupCycleIndex(2), rootedTreesSeries(12)))
\\ => [1, 1, 3, 6, 16, 37, 96, 239, 622, 1607, 4235] = A027852

\\ Linear forests of rooted trees
NumUnlabeledObjsSeq(1/(1-rootedTreesSeries(10))-1)
\\ => [1, 2, 5, 13, 35, 95, 262, 727, 2033, 5714] = A000107

\\ Number of rooted trees with nodes of 2 colors
ColoringsSeq(rootedTreesSeries(10),2)
\\ => [2, 4, 14, 52, 214, 916, 4116, 18996, 89894, 433196] = A038055

\\ Number of rooted trees with colored nodes of exactly 2 colors
ColoringsExactlySeq(rootedTreesSeries(10),2)
\\ => [0, 2, 10, 44, 196, 876, 4020, 18766, 89322, 431758] = A339642

\\ Number of rooted trees with integer labeled nodes covering initial interval of positive integers
NormalLabelingsSeq(rootedTreesSeries(10))
\\ => [1, 3, 21, 214, 3004, 53696, 1169220, 30017582, 887835091, 29728120594] = A339644

\\ Number of rooted trees with colored nodes using k colors 
ColoringsTriangle(rootedTreesSeries(8))
\\ => A141610

\\ Similar triangle with the root not colored
ColoringsTriangle(rootedTreesSeries(10)/(x*sv(1))-1)
\\ => A256064

\\ A directed cycle of trees is a connected function
my(N=10); NumUnlabeledObjsSeq(sSubstOp(cycGroupSeries(N)-1, rootedTreesSeries(N)))
\\ => [1, 2, 4, 9, 20, 51, 125, 329, 862, 2311] = A002861

\\ On labeled nodes
my(N=10); NumLabeledObjsSeq(sSubstOp(cycGroupSeries(N)-1, rootedTreesSeries(N)))
\\ => [1, 3, 17, 142, 1569, 21576, 355081, 6805296, 148869153, 3660215680] = A001865

\\ An undirected cycle of trees is graph with a single cycle (>=1 allowed)
my(N=10); NumUnlabeledObjsSeq(sSubstOp(dihedralGroupSeries(N)-1, rootedTreesSeries(N)))
\\ => [1, 2, 4, 9, 20, 49, 118, 300, 765, 1998] = A068051

\\ Number of connected unicyclic graphs
my(N=12); NumUnlabeledObjsSeq(sSubstOp(dihedralGroupSeries(N)-Pol(dihedralGroupSeries(2)), rootedTreesSeries(N)))
\\ => [1, 2, 5, 13, 33, 89, 240, 657, 1806, 5026] = A001429

\\ Number of labeled connected unicyclic graphs
my(N=10); NumLabeledObjsSeq(sSubstOp(dihedralGroupSeries(N)-Pol(dihedralGroupSeries(2)), rootedTreesSeries(N)))
\\ => [1, 15, 222, 3660, 68295, 1436568, 33779340, 880107840] = A057500

\\ Number of unlabeled unrooted trees (applying Otter's theorem)
NumUnlabeledObjsSeq(sSubstOp(x*sv(1) + x^2*symGroupCycleIndex(2) - (x*sv(1))^2, rootedTreesSeries(15)))
\\ => [1, 1, 1, 2, 3, 6, 11, 23, 47, 106, 235, 551, 1301, 3159, 7741] = A000055

\\ Number of labeled unrooted trees (applying Otter's theorem)
NumLabeledObjsSeq(sSubstOp(x*sv(1) + x^2*symGroupCycleIndex(2) - (x*sv(1))^2, rootedTreesSeries(10)))
\\ => [1, 1, 3, 16, 125, 1296, 16807, 262144, 4782969, 100000000] = A000272

\\ ============================================================
\\ Examples on species of graph.
\\ The following is an adaptation of code in A000088.
edges(v) = {sum(i=2, #v, sum(j=1, i-1, gcd(v[i], v[j]))) + sum(i=1, #v, v[i]\2)}
graphsCycleIndex(n)={my(s=0); forpart(p=n, s+=permcount(p)*2^edges(p)*sMonomial(p)); s/n!}
graphsSeries(n)={sum(k=0, n, graphsCycleIndex(k)*x^k) + O(x*x^n)}

graphsSeries(3)
\\ => 1 + s1*x + (s1^2 + s2)*x^2 + (4/3*s1^3 + 2*s2*s1 + 2/3*s3)*x^3 + O(x^4)

\\ Number of labeled graphs.
NumLabeledObjsSeq(graphsSeries(8))
\\ => [1, 1, 2, 8, 64, 1024, 32768, 2097152, 268435456] = A006125

\\ Number of unlabeled graphs.
NumUnlabeledObjsSeq(graphsSeries(10))
\\ => [1, 1, 2, 4, 11, 34, 156, 1044, 12346, 274668, 12005168] = A000088

\\ Number of connected labeled graphs.
NumLabeledObjsSeq(1 + sLog(graphsSeries(8)))
\\ => [1, 1, 1, 4, 38, 728, 26704, 1866256, 251548592] = A001187

\\ Number of connected unlabeled graphs.
NumUnlabeledObjsSeq(1 + sLog(graphsSeries(10)))
\\ => [1, 1, 1, 2, 6, 21, 112, 853, 11117, 261080, 11716571] = A001349

\\ sLog is the inverse transform of sExp 
my(p=graphsSeries(10)); sExp(sLog(p)) - p
\\ =>  O(x^11)

\\ sExp is the same as sSubstOp with the species of symmetric group as the first argument.
my(N=10, p=sLog(graphsSeries(N))); NumUnlabeledObjsSeq(sSubstOp(symGroupSeries(N), p))
\\ => [1, 1, 2, 4, 11, 34, 156, 1044, 12346, 274668, 12005168] = A000088

\\ Number of unlabeled graphs with distinct components.
NumUnlabeledObjsSeq(sWeighT(sLog(graphsSeries(10))))
\\ => [1, 1, 1, 3, 8, 29, 142, 1005, 12173, 273582, 11992634] = A207828

\\ Number of unlabeled rooted graphs.
NumUnlabeledObjsSeq(sPoint(graphsSeries(10)))
\\ => [1, 2, 6, 20, 90, 544, 5096, 79264, 2208612, 113743760] = A000666

\\ Number of graphs rooted at two noninterchangeable vertices.
\\ Double pointing doesn't work (unless both roots may coincide)
\\ (see implementation of sPoint - need to do both derivatives first) 
my(s1=sv(1)); NumUnlabeledObjsSeq(s1^2*deriv(deriv(graphsSeries(10), s1), s1))
\\ => [2, 8, 40, 240, 1992, 24416, 483040, 16343872, 982635280] = A304070

\\ Number of graphs rooted at two indistinguishable vertices.
\\ Use superposition to put the roots into one set.
NumUnlabeledObjsSeq(sCartProd(graphsSeries(10), x^2*symGroupCycleIndex(2)*symGroupSeries(8)))
\\ => [2, 6, 28, 148, 1144, 13128, 250240, 8295664, 494367376] = A303829

\\ Number of labeled connected rooted graphs.
NumLabeledObjsSeq(sPoint(sLog(graphsSeries(8))))
\\ => [1, 2, 12, 152, 3640, 160224, 13063792, 2012388736] = A053549

\\ Number of labeled connected rooted graphs, by another route.
NumLabeledObjsSeq(sPoint(graphsSeries(8))/graphsSeries(8))
\\ => [1, 2, 12, 152, 3640, 160224, 13063792, 2012388736] = A053549

\\ Number of unlabeled connected rooted graphs.
NumUnlabeledObjsSeq(sPoint(sLog(graphsSeries(10))))
\\ => [1, 1, 3, 11, 58, 407, 4306, 72489, 2111013, 111172234] = A126100


\\ ============================================================