/* A132860 method:
The basic idea is to group the trees according to their coefficients mod k,
then build a polynomial to count each class based on the initial value.

Unfortunately, mod k does not give enough information to determine the later
coefficients. Theorem 4 from [1] states that the coefficients mod k*A007947(k)
are uniquely determined by the kth root mod A007947(k).  We need enough leading
terms mod k*A007947(k) to fix the rest mod k (except the last term, which does
not matter).  The number of leading terms needed is ceil(n/A020639(k)).

The counting polynomials are built backwards from the last term.
Write f(x,m) for the number of trees fixed up to term m, which is x.
Within the mod k part, f(x,m-1) = sum(i=1..x, f(i*k-p,m)) for appropriate p.
Within the leading terms the situation is only slightly more complicated.
Finally we evaluate f(k,2).

A132860_simple method:
This implementation is meant as a double-check since it is very slow.
It builds the trees directly using the rule that potential values differ by k.
Two levels of optimisation are available which cut out the last one or two loops.

[1] N. Heninger, E. M. Rains and N. J. A. Sloane, 
<a href="http://arXiv.org/abs/math.NT/0509316">
On the Integrality of n-th Roots of Generating Functions</a>,
J. Combinatorial Theory, Series A, 113 (2006), 1732-1745.

Martin Fuller
*/

IndefiniteSum(X)=
{
  if (#IndefiniteSum_Matrix < #X,
    IndefiniteSum_Matrix =
      concat(IndefiniteSum_Matrix, vector(#X-#IndefiniteSum_Matrix));
  );
  if (type(IndefiniteSum_Matrix[#X]) <> "t_MAT",
    IndefiniteSum_Matrix[#X] = matrix(#X,#X, i,j, 
      if(j<=i, (-1)^(i-j)*binomial(#X+1-j,#X-i), 0)
    )^-1;
  );
  Pol(concat(IndefiniteSum_Matrix[#X] * Col(X), [0]))
}
IndefiniteSum_Matrix = [];
addhelp(IndefiniteSum, "IndefiniteSum(X):compute the sum of polynomial X starting from 1.");

A007947(n)=local(p);p=factor(n)[,1];prod(i=1,#p,p[i])
A020639(n)=if(n==1,1,factor(n)[1,1])

A132860(n,k)=
{
  local(result, mu, lead, a, count, p, q);
  if (k == 1, return(1));
  mu = k*A007947(k);
  lead = ceil(n/A020639(k));
  forvec (b=vector(n, i, if(i<=2, [1,1], if(i<=lead, [0,mu/k-1], [0,0]))),
    a = Vec(Ser(b)^k);
    count = 1;
    forstep(i=n, max(lead+1,3), -1,
      p = (-a[i]) % k;
      count = IndefiniteSum(subst(count, 'x, 'x*k - p));
    );
    forstep(i=lead, 3, -1,
      p = (-a[i]) % mu;
      q = a[i-1] % (mu/k);
      count = IndefiniteSum(subst(count, 'x, 'x*mu - p));
      count = subst(count, 'x, ('x - q) * k / mu + (q*k + p >= mu));
    );
    result += subst(count, 'x, k);
  );
  result
}

A132860_simple(n,k,optimise=0,v=[1])=
{
  if(n==#v, if(default(debug), print(v)); 1,
  if(optimise>=1 & n==#v+1, v[#v],
    sum(i=0, k-1, if(denominator(polcoeff(Ser(concat(v,-i))^(1/k), #v)) == 1,
      if(optimise>=2 & n==#v+2, k*v[#v]*(v[#v]+1)/2 - v[#v]*i,
        sum(j=1, v[#v], A132860_simple(n,k,optimise,concat(v,j*k-i))))))))
}