\\ ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
\\ A361953: Unlabeled weakly graded posets.
\\ Author: Andrew Howroyd, Apr 2023
\\ ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
\\ The labeled version of A361952 is A361950.
\\ The method here is at a high level based on the labeled enumeration.
\\ Instead of e.g.f's, a mixture of cycle index series and o.g.f's are 
\\ used. The cycle index series introduce some complexity.
\\ The reference for the labeled enumeration is:
\\ D. A. Klarner, The number of graded partially ordered sets, J. Combin. Theory, 6 (1969), 12-19.
\\


\\ 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}


\\ This and C(n) compute A361952 as a vector of column g.f.
\\ src is a cycle index series represented as a vector of vectors.
\\ Each entry in src is a vector of length numbpart(n-1) containing the cycle-index coefficients
\\ in partition order. The cycle index series coefficients are power series in x counting the part
\\ of the poset that it is not covered by the top layer. The coefficients are scaled to be integers.
\\ The return value is also a similar cycle index series.
AddLayer(src,x='x)={
  my(n=#src-1, results=vector(n+1));
  for(j=0, n,
    my(v=vector(numbpart(j)));
    for(i=0, n-j,
      my(ui=src[i+1], Q=vector(i), xj=0, A=O(x^(n-j-i+1)));
      forpart(pj=j, xj++;
        for(s=1, #Q, Q[s]=2^sum(t=1, #pj, gcd(s, pj[t])));
        my(xi=0, t=0);
        forpart(pi=i, xi++;
           t += (ui[xi] + A)*prod(t=1, #pi, Q[pi[t]]);
        ); \\ pi
        v[xj] += t * x^i * permcount(pj)/i!
      ); \\ pj
    ); \\ i
    results[j+1] = v;
  ); \\ j
  results;
}

\\ Computaton of A361952, n rows (elements), m columns (height).
C(n,m=n)={
  my(data=vector(n+1, i, vector(numbpart(i-1))), c=vector(m+1)); 
  c[1] = 1 + O(x*x^n); data[1][1] = c[1];
  for(h=1, m, 
     data = AddLayer(data);
     c[h+1] = sum(i=0, n, x^i*vecsum(data[i+1])/i!);
  );
  c;
}

\\ Inverse Euler transforms
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)}

\\ The sequences
A361952tabl(n,m=n) = {my(c=C(n,m)); Mat(vector(#c, i, Col(c[i])))}
A361952col(k,n) = {Vec(C(n,k)[k+1])}
A361952row(k,n) = {[polcoef(p, k) | p<-C(k,n)]} 

A361953tabl(n)={my(c=C(n), b=vector(n+1, h, c[h]/c[max(h-1, 1)])); Mat(vector(n+1, h, Col(b[h]-if(h>1, b[h-1]), -n-1)))}
A361953col(k,n)={my(c=C(n+k-1,k)); Vec(c[k+1]/c[max(k,1)] - if(k>0, c[k]/c[max(k-1,1)]))}
A362953partsums(k,n)={my(c=C(n,k)); Vec(c[k+1]/c[max(k,1)])}
A361920seq(n)={my(c=C(n)); Vec(c[n+1]/c[n])}

A361954tabl(n)={my(c=C(n), b=vector(n+1, h, InvEulerMTS(c[h]/c[max(h-1, 1)]))); Mat(vector(n, h, Col(b[h+1]-b[h], -n)))}
A361954col(k,n)={my(c=C(n+k-1,k)); Vec(InvEulerMTS(c[k+1]/c[max(k,1)]) - if(k>0, InvEulerMTS(c[k]/c[max(k-1,1)])))}
A362954partsums(k,n)={my(c=C(n,k)); Vec(1 + InvEulerMTS(c[k+1]/c[max(k,1)]))}
A361955seq(n)={my(c=C(n)); Vec(1 + InvEulerMTS(c[n+1]/c[n]))}

\\ End