\\ ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
\\ Strongly Connected Digraphs
\\ Author: Andrew Howroyd, Jan 2022
\\ ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

\\ If intending to compute more than a few terms then give PARI more memory!
\\ This can be done with 'default'. For example:
\\ default(parisizemax, 10^9)

\\ INTRODUCTION
\\ ------------
\\ Strong graphs were first enumerated by Liskovets. The method used here
\\ includes significant (and ingeneous) improvements by Wright and R W Robinson
\\ (and/or perhaps Read). [Warning: This is not intended as an accurate or complete
\\ history]
\\ The primary referene paper here is:
\\ R. W. Robinson, Counting digraphs with restrictions on the strong components,
\\ Combinatorics and Graph Theory '95 (T.-H. Ku, ed.), World Scientific, Singapore (1995), 343-354.
\\ Although this paper deals with the labeled case and the additional complexity of
\\ a source and sink, it contains the method used here as a simple equation:
\\
\\ In particular eqn (18) states:
\\    Ggf(exp(-C(x))) * Ggf(G(x)) = 1,
\\ where C(x) is the e.g.f. of any set of strongly connected digraphs
\\ and G(x) is the e.g.f. is the set of graphs whose strongly components are C(x)
\\ and Ggf(), which Robertsons denotes by a triangle, is the operator
\\ that transforms an e.g.f. into a g.g.f. (graphical generating function). 
\\ Although the equation is written for the labeled case represented by
\\ their e.g.f.'s, it applies equaly to combinatorial classes represented by cycle index series. 
\\
\\ Some particular cases:
\\     Ggf(exp(-Singleton)) * Ggf(AcyclicDigraphs) = 1,
\\     Ggf(exp(-StronglyConnectedDigraphs)) * Ggf(Digraphs) = 1,
\\     Ggf(exp(-StronglyConnectedOrientedDiagraphs)) * Ggf(OrientedDigraphs) = 1.
\\
\\ The 2nd equation can be rewritten:
\\     StronglyConnectedDigraphs = -log(UnGcf(1/Ggf(Digraphs))).
\\
\\ The following PARI which is commented out here quickly demonstrates the power of this
\\ equation for the labeled case:
\\ The Ggf operator just divides the terms of a series by 2^(k*(k-1)/2) and 
\\ UnGgf does the oppossite.
/*
Ggf(p)={my(n=serprec(p,x)-1); serconvol(p, sum(k=0, n, x^k/2^(k*(k-1)/2), O(x*x^n)))}
UnGgf(p)={my(n=serprec(p,x)-1); serconvol(p, sum(k=0, n, x^k*2^(k*(k-1)/2), O(x*x^n)))}
StrongComponents(graphEgf)={-log(UnGgf(1/Ggf(graphEgf)))}
DigraphEgf(n)={sum(k=0, n, 2^(k*(k-1))*x^k/k!, O(x*x^n) )}

\\ sequence A053763 - number of digraphs.
Vec(serlaplace(DigraphEgf(12)))
\\ sequence A003030 - number of strongly connected digraphs.
Vec(serlaplace(StrongComponents(DigraphEgf(12))))
*/
\\
\\ The unlabeled case follows the same lines except there is no
\\ equivalent to a Ggf (to my knowledge) so 'UnGgf(1/Ggf(graphEgf))'
\\ is performed as a single operation 'InvGgfCIData'. Instead
\\ of working with e.g.f's a cycle index series is required.
\\ (See A339645 for some PARI code on cycle index series represented as PARI power series)
\\ Here we use a more efficient array representation because the all
\\ important operation has no quick implementation using standard power series functions.

\\ START OF ACTUAL CODE.

\\ Concatenates two sorted small vectors, preserving order.
\\ Used exclusively for combining two partitions.
\\ For example: Merge([1,1,4], [3,6]) => [1,1,3,4,6]
Merge(p1, p2)={
  my(v=vectorsmall(#p1+#p2), i=1, j=1, k=1);
  while(i<=#p1 && j<=#p2, if(p1[i]<=p2[j], v[k]=p1[i]; i++, v[k]=p2[j]; j++); k++);
  while(i<=#p1, v[k]=p1[i]; i++; k++);
  while(j<=#p2, v[k]=p2[j]; j++; k++); 
  v;
}

\\ Creates a map of partition to rank.
\\ The rank of the partition is the order it appears in the generation order of paritions.
\\ This is used to support matrices indexed by partition.
\\ Although there are more memory efficient ways of doing this taking advantage of
\\ the particular ordering of partiions, this is a small part of the overall memory footprint.
BuildPartitionIndex(n)={
   my(M=Map(), ix=1);
   forpart(p = n, mapput(M, p, ix); ix++);
   M;
}

\\ Gives the number of permutations of a particular type where the type is a partition.
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}


\\ Functions for constructing cycle index series for various kinds of graph.
\\ The parameter v is a partition (of the cycle type)
\\ The parameter yf is to allow edge data to be captured if required.
DigraphEdges(v, yf) = {prod(i=2, #v, prod(j=1, i-1, my(g=gcd(v[i], v[j])); yf(v[i]*v[j]/g)^(2*g))) * prod(i=1, #v, yf(v[i])^(v[i]-1))}
DigraphWithLoopEdges(v, yf) = {prod(i=2, #v, prod(j=1, i-1, my(g=gcd(v[i], v[j])); yf(v[i]*v[j]/g)^(2*g))) * prod(i=1, #v, yf(v[i])^v[i])}
OrientedGraphEdges(v, yf) = {prod(i=2, #v, prod(j=1, i-1, my(g=gcd(v[i], v[j])); yf(v[i]*v[j]/g)^g )) * prod(i=1, #v, my(c=v[i]); yf(c)^((c-1)\2))}

\\ One of the above can be passed to this function to actually create the cycle index series.
\\ The cycle index series is represented as vector of vectors.
GraphCIData(n, edgeFunc, yf=e->2)={
  my(results=vector(n+1));
  for(n=0, n,
     my(v=vector(numbpart(n)), xi=0);
     forpart(p=n, xi++;
        v[xi] = permcount(p)*edgeFunc(p, yf);
     );
     results[1+n] = v;
  );
  results;
}

\\ This is the unlabeled equivalent of InvGgf(1/Ggf( graphEgfSeries ))
\\ It is its own inverse. This is the core algorithm for finding the number of strong graphs.
\\ Compare with 'DagCycleIndexData' function given in A122078.
InvGgfCIData(blocks,yf=e->2)={
  my(n=#blocks-1, results=vector(n+1));
  results[1] = [1];
  for(n=1, n,
    my(v=vector(numbpart(n)), pm=BuildPartitionIndex(n));
    for(i=1, n,
      my(j=n-i, uj=results[j+1], ui=blocks[i+1], b=binomial(n,i), Q=vector(i), xj=0);
      forpart(pj=j, xj++;
        for(s=1, #Q, Q[s]=prod(t=1, #pj, my(g=gcd(s, pj[t])); yf(s*pj[t]/g)^g));
        my(xi=0);
        forpart(pi=i, xi++;
           my(col=mapget(pm, Merge(pi,pj)));
           v[col] -= uj[xj]*b*ui[xi]*prod(t=1, #pi, Q[pi[t]]);
        ); \\ pi
      ); \\ pj
    ); \\ i
    results[n+1] = v;
  ); \\ n
  results;
}

\\ Helper to convert a cycle index series stored as a vector of vectors into an o.g.f. counting the unlabeled objects.
UnlabeledOgf(data)={
  sum(n=0, #data-1, vecsum(data[n+1])*x^n/n!, O(x^#data));
}


\\ The only other operation that is required is the combinatorial logarithm. A cycle index series implementation of this is given
\\ in A339645 (functions for combinatorial species). Here we cheat, by first converting to an o.g.f. and then taking
\\ the inverse euler transform which is the o.g.f. equivalent of combinatorial logorithm. 

\\ Inverse Euler transforms (single variable and multivariate)
\\ The single variate version takes a vector and returns a vector of same length.
InvEulerT(v)={my(p=log(1+x*Ser(v))); dirdiv(vector(#v,n,polcoef(p,n)), vector(#v,n,1/n))}
\\ The multivariate version works with power series. The primary variable must be x.
InvEulerMTS(p)={my(n=serprec(p,x)-1, q=log(p), vars=variables(p)); sum(i=1, n, moebius(i)*substvec(q + O(x*x^(n\i)), vars, apply(v->v^i,vars))/i)}


\\ Sequences for strongly connected oriented graphs.
A350489seq(n)={concat([1], -InvEulerT(Vec(-1 + UnlabeledOgf(InvGgfCIData(GraphCIData(n, OrientedGraphEdges, e->3))))))}
A350750rows(n)={my(v=Vec(-InvEulerMTS(UnlabeledOgf(InvGgfCIData(GraphCIData(n, OrientedGraphEdges, e->1+2*y^e), e->1+y^e))))); vector(#v, n, Vecrev(v[n], 1+n*(n-1)/2)) }
A350751seq(n)={my(g=-InvEulerMTS(UnlabeledOgf(InvGgfCIData(GraphCIData(n, OrientedGraphEdges, e->1+2*y^e + O(y*y^n)), e->1+y^e+O(y*y^n))))); Vec(sum(i=1, n, polcoef(g,i,x)), -n)}

\\ Sequences for strongly connected digraphs.
A035512seq(n)={concat([1], -InvEulerT(Vec(-1 + UnlabeledOgf(InvGgfCIData(GraphCIData(n, DigraphEdges))))))}
A057276rows(n)={[Vecrev(p) | p<-Vec(-InvEulerMTS(UnlabeledOgf(InvGgfCIData(GraphCIData(n, DigraphEdges, e->1+y^e), e->1+y^e))))]}
A350753rows(n)={my(yf(e)=1+y^e+O(y*y^n)); my(g=-InvEulerMTS(UnlabeledOgf(InvGgfCIData(GraphCIData(n, DigraphEdges, yf), yf)))); vector(n, i, Vec(O(x*x^i)+polcoef(g,i-1,y), -i)) }
A350752seq(n)={my(yf(e)=1+y^e+O(y*y^n)); my(g=-InvEulerMTS(UnlabeledOgf(InvGgfCIData(GraphCIData(n, DigraphEdges, yf), yf)))); Vec(sum(i=1, n, polcoef(g,i,x)), -n)}

\\ Sequences for strongly connected digraphs with multiple edges and loops.
A139622rows(n)={my(yf(e)=1/(1-y^e)+O(y*y^n)); my(g=-InvEulerMTS(UnlabeledOgf(InvGgfCIData(GraphCIData(n, DigraphWithLoopEdges, yf), yf)))); vector(n, i, Vec(O(x*x^i)+polcoef(g,i,y), -i)) }
A139627seq(n)={my(A=A139622rows(n)); concat([1], vector(#A, n, vecsum(A[n])))}

\\ The following is not used here, but converts a cycle index series represented as a vector of vectors
\\ into a power series. Depends on sMonomial function defined in A339645 PARI script.
ToCycleIndexSeries(data)={
  sum(n=0, #data-1,
    my(row=data[n+1], t=0, xi=0);
    forpart(pi=n, xi++; t += row[xi]*sMonomial(pi) );
    x^n*t/n!,
    O(x^#data)
  );
}


\\ End