jerzy.karczmarczuk at info.unicaen.fr jerzy.karczmarczuk at info.unicaen.fr
Thu Nov 8 10:30:01 EST 2007

```Henning Thielemann writes:

> And now we have much Haskell code for one sequence to be submitted to the
> Online Encyclopedia of Integer Sequences!

Is there anything really there in Haskell?...

Well, if you are interested in something more "venerable" than rabbits, why
not try the sequence which gives the number of rooted trees:

T_n = 0,1,1,2,4,9,20,48,115,286,719,1842,4766,12486,32973,87811,...

Such combinatorial coefficients are sometimes useful in practice.
(Even in theoretical physics...)
Their generating function T(x) = Sum_{n >= 1} T_n*x^n
fulfils the Polya relation:

T(x) = x*exp(T(x) + T(x^2)/2 + T(x^3)/3 + T(x^4)/4+...+ T(x^m)/m + ...)

from which you see clearly that you don't see anything clearly.

For T_n there is a recurrence:

T_(n+1) = (1/n) * sum_{k=1..n} ( sum_{d|k} d*T_(d) ) * T_(n-k+1).

(where d|k means the sum over the divisors of k)
whose advantage is that is really so awfully ugly, that it looks
professional... And it seems that the formulae in Maple, etc., on the
page of Sloane, are based essentially on that.

==

I tried to code the stuff *directly* from the Polya identity. I'll give
the solution below, but, perhaps, try yourself? It is not a one-liner,
though, and uses a small package for formal series manipulation. Recall
that if a series u_0 + u_1*x + u_2*x^2 + ... is represented by a list
[u_0, u_1,...], then the sum and the difference are just zipWith (+/-).

The multiplication is: (u0:uq)*v@(v0:vq) = u0*v0 : u0*>vq + v*uq
and the division:
(u0:uq) / v@(v0:vq) = let w0 = u0/v0
in  w0:(uq - w0*>vq)/v
where (*>) = map . (*), multiplies a series by a scalar.

Now, we have here an exponential, so it seems that a Floating instance of
the series is necessary. However, in order to prevent some people from
pointing out that this would break down, let's define a *restricted*
exponential, exp0(s) for s such that s_0=0. Then, if s is rational, the
exponential will also be rational. We shall use the ideology described
well in Knuth. If w=exp(s), then w' = w*s', and
w = (exp s_0) + INTEG w*s'
which is a co-recursive relation for w. In Haskell:

intgs = nt 1 where nt n = n : nt (n+1)
sDif (_:sq) = zipWith (*) sq intgs   -- differentiation of a series
sInt c s = c : zipWith (/) s intgs   -- integration, with the int. constant

Don't ask me why not [1 ..]. This forces an instance of Enum, I don't
want to hear about. Now, the series exponential:

exp0 u@(0:_) = w where w = sInt 1 (sDif u * w)  -- otherwise it fails.

and finally, the Polya formula. Notice the following:

a. The T(x) starts with x*..., so the zeroth coeff is zero. The restricted
exponential should work.
b. T(x^k) puts some zeros between the elements of T(x), that's all. This
is done by the function (\x->x : replicate...) below.
c. The infinite sum of T(x^m)/m is *co-recursively doable!!*

So, finally,   tau = 0 : psi where

psi = exp0(foldr1 (\x s->0:x+s)           -- this (+) is for series.
(map (\m->(1/fromInteger m)*>
concatMap (\x->x : replicate (fromInteger(m-1)) 0) psi) intgs))

So, after all, if the auxiliary series functions are considered to be a part
of our "standard" environment, we got a one-liner, although the length of
this line is substantial.

Anybody tries to shorten this? If you don't, I'll put it on my grave stone,
but hopefully not yet tomorrow.

Jerzy Karczmarczuk

```