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