Paraffins: Difference between revisions

From Rosetta Code
Content added Content deleted
m (added an image of a paraffins: isopentane (common name).)
Line 964: Line 964:
249: 5814271898167303040368103945830220447130073898083466852225709084407144308593691069932064987528870826155297
249: 5814271898167303040368103945830220447130073898083466852225709084407144308593691069932064987528870826155297
250: 16206624309085062837751018464745815688226709117091506494175397665527493805947344857313038875654104100026504</pre>
250: 16206624309085062837751018464745815688226709117091506494175397665527493805947344857313038875654104100026504</pre>

=={{header|PARI/GP}}==

This function is for recent PARI/GP:
<lang parigp>paraffin(p) =
{
local (P = p+1, R, U = R = Vec([1,1], P));

for (n = 1, p,
((B,n,C,S,l=n) -> my(b,c,i,s);
for (b = 1, 4-B,
if ((s = S + b * n) < P,
c = R[n+1] * C * prod(i = 1, b-1, (R[n+1]+i) / (i+1)) ;
if (l+l < s, U[s+1] += c);
if (B+b < 4, R[s+1] += c; i = n; while (i--, self()(B+b, i, c, s, l)))))
)(0,n,1,1);
if (n % 2,, U[n+1] += R[n/2+1] * (R[n/2+1]+1) / 2);
print([n, U[n+1]]))
}</lang>


Code for older version of PARI/GP < 2.9:
<lang parigp>iso(B,n,C,S,l=n) =
{
my(b,c,i,s);

for (b = 1, 4-B,
if ((s = S + b * n) < P,
c = R[n+1] * C * prod(i = 1, b-1, (R[n+1]+i) / (i+1));
if (l+l < s, U[s+1] += c);
if (B+b < 4, R[s+1] += c; i = n; while (i--, iso(B+b, i, c, s, l)))))
}

paraffin(p) =
{
local (P = p+1, R, U = R = Vec([1,1], P));

for (n = 1, p, iso(0, n, 1, 1);
if (n % 2,, U[n+1] += R[n/2+1] * (R[n/2+1]+1) / 2);
print([n, U[n+1]]))
}</lang>


Output:<pre>paraffin(32)

[1, 1]
[2, 1]
[3, 1]
[4, 2]
[5, 3]
[6, 5]
[7, 9]
[8, 18]
[9, 35]
[10, 75]
[11, 159]
[12, 355]
[13, 802]
[14, 1858]
[15, 4347]
[16, 10359]
[17, 24894]
[18, 60523]
[19, 148284]
[20, 366319]
[21, 910726]
[22, 2278658]
[23, 5731580]
[24, 14490245]
[25, 36797588]
[26, 93839412]
[27, 240215803]
[28, 617105614]
[29, 1590507121]
[30, 4111846763]
[31, 10660307791]
[32, 27711253769]
</pre>


=={{header|Pascal}}==
=={{header|Pascal}}==

Revision as of 14:40, 16 December 2016

Task
Paraffins
You are encouraged to solve this task according to the task description, using any language you may know.

This organic chemistry task is essentially to implement a tree enumeration algorithm.


Task

Enumerate, without repetitions and in order of increasing size, all possible paraffin molecules (also known as alkanes).


Paraffins are built up using only carbon atoms, which has four bonds, and hydrogen, which has one bond. All bonds for each atom must be used, so it is easiest to think of an alkane as linked carbon atoms forming the "backbone" structure, with adding hydrogen atoms linking the remaining unused bonds.

In a paraffin, one is allowed neither double bonds (two bonds between the same pair of atoms), nor cycles of linked carbons. So all paraffins with   n   carbon atoms share the empirical formula   CnH2n+2

But for all   n ≥ 4   there are several distinct molecules ("isomers") with the same formula but different structures.

The number of isomers rises rather rapidly when   n   increases.

In counting isomers it should be borne in mind that the four bond positions on a given carbon atom can be freely interchanged and bonds rotated (including 3-D "out of the paper" rotations when it's being observed on a flat diagram), so rotations or re-orientations of parts of the molecule (without breaking bonds) do not give different isomers. So what seem at first to be different molecules may in fact turn out to be different orientations of the same molecule.


Example

With   n = 3   there is only one way of linking the carbons despite the different orientations the molecule can be drawn;   and with   n = 4   there are two configurations:

  • a   straight chain:     (CH3)(CH2)(CH2)(CH3)
  • a branched chain:     (CH3)(CH(CH3))(CH3)


Due to bond rotations, it doesn't matter which direction the branch points in.

The phenomenon of "stereo-isomerism" (a molecule being different from its mirror image due to the actual 3-D arrangement of bonds) is ignored for the purpose of this task.

The input is the number   n   of carbon atoms of a molecule (for instance 17).

The output is how many different different paraffins there are with   n   carbon atoms (for instance   24,894   if   n = 17).

The sequence of those results is visible in the Sloane encyclopedia. The sequence is (the index starts from zero, and represents the number of carbon atoms): 

1, 1, 1, 1, 2, 3, 5, 9, 18, 35, 75, 159, 355, 802, 1858, 4347, 10359,
24894, 60523, 148284, 366319, 910726, 2278658, 5731580, 14490245,
36797588, 93839412, 240215803, 617105614, 1590507121, 4111846763,
10660307791, 27711253769, ...


Extra credit

Show the paraffins in some way.

A flat 1D representation, with arrays or lists is enough, for instance:

<lang haskell>*Main> all_paraffins 1 

[CCP H H H H]
  • Main> all_paraffins 2
[BCP (C H H H) (C H H H)]
  • Main> all_paraffins 3
[CCP H H (C H H H) (C H H H)]
  • Main> all_paraffins 4
[BCP (C H H (C H H H)) (C H H (C H H H)),CCP H (C H H H) (C H H H)
(C H H H)]
  • Main> all_paraffins 5
[CCP H H (C H H (C H H H)) (C H H (C H H H)),CCP H (C H H H)
(C H H H) (C H H (C H H H)),CCP (C H H H) (C H H H) (C H H H)
(C H H H)]
  • Main> all_paraffins 6
[BCP (C H H (C H H (C H H H))) (C H H (C H H (C H H H))),BCP
(C H H (C H H (C H H H))) (C H (C H H H) (C H H H)),BCP (C H
(C H H H) (C H H H)) (C H (C H H H) (C H H H)),CCP H (C H H H)
(C H H (C H H H)) (C H H (C H H H)),CCP (C H H H) (C H H H)
(C H H H) (C H H (C H H H))]</lang>

Showing a basic 2D ASCII-art representation of the paraffins is better; for instance (molecule names aren't necessary): 

 Methane         Ethane              Propane              Isobutane
                 
    H             H   H             H   H   H             H   H   H
    |             |   |             |   |   |             |   |   |
H - C - H     H - C - C - H     H - C - C - C - H     H - C - C - C - H
    |             |   |             |   |   |             |   |   |
    H             H   H             H   H   H             H   |   H
                                                              |
                                                          H - C - H
                                                              |
                                                              H


Links
  • A paper that explains the problem and its solution in a functional language:

http://www.cs.wright.edu/~tkprasad/courses/cs776/paraffins-turner.pdf

  • A Haskell implementation:

https://github.com/ghc/nofib/blob/master/imaginary/paraffins/Main.hs

  • A Scheme implementation:

http://www.ccs.neu.edu/home/will/Twobit/src/paraffins.scm

  • A Fortress implementation:

http://java.net/projects/projectfortress/sources/sources/content/ProjectFortress/demos/turnersParaffins0.fss?rev=3005

C

Can't show the tree shapes; count only. <lang c>#include <stdio.h>

  1. define MAX_N 33 /* max number of tree nodes */
  2. define BRANCH 4 /* max number of edges a single node can have */

/* The basic idea: a paraffin molecule can be thought as a simple tree 

  with each node being a carbon atom.  Counting molecules is thus the
  problem of counting free (unrooted) trees of given number of nodes.

  An unrooted tree needs to be uniquely represented, so we need a way
  to cannonicalize equivalent free trees.  For that, we need to first
  define the cannonical form of rooted trees.  Since rooted trees can
  be constructed by a root node and up to BRANCH rooted subtrees that
  are arranged in some definite order, we can define it thusly:
    * Given the root of a tree, the weight of each of its branches is
      the number of nodes contained in that branch;
    * A cannonical rooted tree would have its direct subtrees ordered
      in descending order by weight;
    * In case multiple subtrees are the same weight, they are ordered
      by some unstated, but definite, order (this code doesn't really
      care what the ordering is; it only counts the number of choices
      in such a case, not enumerating individual trees.)

  A rooted tree of N nodes can then be constructed by adding smaller,
  cannonical rooted trees to a root node, such that:
    * Each subtree has fewer than BRANCH branches (since it must have
      an empty slot for an edge to connect to the new root);
    * Weight of those subtrees added later are no higher than earlier
      ones;
    * Their weight total N-1.
  A rooted tree so constructed would be itself cannonical.

  For an unrooted tree, we can define the radius of any of its nodes:
  it's the maximum weight of any of the subtrees if this node is used
  as the root.  A node is the center of a tree if it has the smallest
  radius among all the nodes.  A tree can have either one or two such
  centers; if two, they must be adjacent (cf. Knuth, tAoCP 2.3.4.4).

  An important fact is that, a node in a tree is its sole center, IFF
  its radius times 2 is no greater than the sum of the weights of all
  branches (ibid).  While we are making rooted trees, we can add such
  trees encountered to the count of cannonical unrooted trees.

  A bi-centered unrooted tree with N nodes can be made by joining two
  trees, each with N/2 nodes and fewer than BRANCH subtrees, at root.
  The pair must be ordered in aforementioned implicit way so that the
  product is cannonical. */

typedef unsigned long long xint;

  1. define FMT "llu"

xint rooted[MAX_N] = {1, 1, 0}; xint unrooted[MAX_N] = {1, 1, 0};

/* choose k out of m possible values; chosen values may repeat, but the 

  ordering of them does not matter.  It's binomial(m + k - 1, k) */

xint choose(xint m, xint k) { xint i, r;

if (k == 1) return m; for (r = m, i = 1; i < k; i++) r = r * (m + i) / (i + 1); return r; }

/* constructing rooted trees of BR branches at root, with at most 

  N radius, and SUM nodes in the partial tree already built. It's
  recursive, and CNT and L carry down the number of combinations
  and the tree radius already encountered. */

void tree(xint br, xint n, xint cnt, xint sum, xint l) { xint b, c, m, s;

for (b = br + 1; b <= BRANCH; b++) { s = sum + (b - br) * n; if (s >= MAX_N) return;

/* First B of BR branches are all of weight n; the rest are at most of weight N-1 */ c = choose(rooted[n], b - br) * cnt;

/* This partial tree is singly centered as is */ if (l * 2 < s) unrooted[s] += c;

/* Trees saturate at root can't be used as building blocks for larger trees, so forget them */ if (b == BRANCH) return; rooted[s] += c;

/* Build the rest of the branches */ for (m = n; --m; ) tree(b, m, c, s, l); } }

void bicenter(int s) { if (s & 1) return;

/* Pick two of the half-size building blocks, allowing repetition. */ unrooted[s] += rooted[s/2] * (rooted[s/2] + 1) / 2; }

int main() { xint n; for (n = 1; n < MAX_N; n++) { tree(0, n, 1, 1, n); bicenter(n); printf("%"FMT": %"FMT"\n", n, unrooted[n]); }

return 0; }</lang>Same idea, with GMP, and done somewhat differently:<lang c>#include <gmp.h>

  1. include <stdio.h>
  2. include <stdlib.h>
  1. define MAX_BRANCH 4
  2. define MAX_N 500

mpz_t bcache[MAX_N + 1]; mpz_t ucache[MAX_N + 1]; mpz_t *rcache[MAX_N + 1][MAX_BRANCH + 1];

mpz_t tmp1, tmp2; void choose(mpz_t r, mpz_t m, int k) { int i; mpz_set(r, m);

mpz_add_ui(tmp1, m, 1); for (i = 1; i < k; ) { mpz_mul(r, r, tmp1); mpz_divexact_ui(r, r, ++i);

if (i >= k) break; mpz_add_ui(tmp1, tmp1, 1); } }

mpz_t rtmp1, rtmp2; void calc_rooted(mpz_t res, int n, int b, int r) { mpz_set_ui(res, 0);

if (n == 1 && b == 0 && r == 0) { mpz_set_ui(res, 1); return; } else if (n <= b || n <= r || n == 1 || b == 0 || r == 0) return;

int b1, r1; for (b1 = 1; b1 <= b && r * b1 < n; b1++) { choose(rtmp1, bcache[r], b1); mpz_set_ui(rtmp2, 0); for (r1 = 0; r1 < r && r1 + r * b1 < n; r1++) mpz_add(rtmp2, rtmp2, rcache[n - r * b1][b - b1][r1]);

mpz_addmul(res, rtmp1, rtmp2); } }

void calc_first_branch(int n) { int b, r; mpz_init_set_ui(bcache[n], 0);

for (b = 0; b < MAX_BRANCH; b++) for (r = 0; r < n; r++) mpz_add(bcache[n], bcache[n], rcache[n][b][r]); }

void calc_unrooted(int n) { int b, r;

for (b = 0; b <= MAX_BRANCH; b++) { mpz_t *p = malloc(sizeof(mpz_t) * n); rcache[n][b] = p; for (r = 0; r < n; r++) { mpz_init(p[r]); calc_rooted(p[r], n, b, r); } }

calc_first_branch(n);

mpz_init_set_ui(ucache[n], 0); for (r = 0; r * 2 < n; r++) for (b = 0; b <= MAX_BRANCH; b++) mpz_add(ucache[n], ucache[n], rcache[n][b][r]);

if (!(n & 1)) { mpz_add_ui(rtmp1, bcache[n/2], 1); mpz_mul(rtmp1, rtmp1, bcache[n/2]); mpz_divexact_ui(rtmp1, rtmp1, 2); mpz_add(ucache[n], ucache[n], rtmp1); } }

void init(void) { mpz_init(tmp1), mpz_init(tmp2); mpz_init(rtmp1), mpz_init(rtmp2); }

int main(void) { int i;

init();

for (i = 0; i <= MAX_N; i++) { calc_unrooted(i); gmp_printf("%d: %Zd\n", i, ucache[i]); }

return 0; }</lang>

D

Translation of: Go

<lang d>import std.stdio, std.bigint;

enum uint nMax = 250; enum uint nBranches = 4;

__gshared BigInt[nMax + 1] rooted = [1.BigInt, 1.BigInt /*...*/], 

                          unrooted = [1.BigInt, 1.BigInt /*...*/];

void tree(in uint br, in uint n, in uint l, in uint inSum, 

         in BigInt cnt) nothrow {
   __gshared static BigInt[nBranches] c;

   uint sum = inSum;
   foreach (immutable b; br + 1 .. nBranches + 1) {
       sum += n;
       if (sum > nMax || (l * 2 >= sum && b >= nBranches))
           return;
       if (b == br + 1) {
           c[br] = rooted[n] * cnt;
       } else {
           c[br] *= rooted[n] + b - br - 1;
           c[br] /= b - br;
       }
       if (l * 2 < sum)
           unrooted[sum] += c[br];
       if (b < nBranches)
           rooted[sum] += c[br];
       foreach_reverse (immutable m; 1 .. n)
           tree(b, m, l, sum, c[br]);
   }

}

void bicenter(in uint s) nothrow { 

   if ((s & 1) == 0)
       unrooted[s] += rooted[s / 2] * (rooted[s / 2] + 1) / 2;

}

void main() { 

   foreach (immutable n; 1 .. nMax + 1) {
       tree(0, n, n, 1, 1.BigInt);
       n.bicenter;
       writeln(n, ": ", unrooted[n]);
   }

}</lang>

Output:
1: 1
2: 1
3: 1
4: 2
5: 3
6: 5
7: 9
8: 18
9: 35
10: 75
11: 159
12: 355
13: 802
14: 1858
15: 4347
16: 10359
17: 24894
18: 60523
19: 148284
20: 366319
21: 910726
...
247: 748434113260252449609376828666343341456610378512725586135955779939163320693444008584504390382026391947130
248: 2086006351917005252913566124773054331962205157167696706926185063169623907656246841866717933958839366769700
249: 5814271898167303040368103945830220447130073898083466852225709084407144308593691069932064987528870826155297
250: 16206624309085062837751018464745815688226709117091506494175397665527493805947344857313038875654104100026504

Run-time with nMax = 250 is about 3.6 seconds (about twice the Go entry using go1.2).

Go

Translation of: C

<lang go>package main

import ( 

   "fmt"
   "math/big"

)

const branches = 4 const nMax = 500

var rooted, unrooted [nMax + 1]big.Int var c [branches]big.Int var tmp = new(big.Int) var one = big.NewInt(1)

func tree(br, n, l, sum int, cnt *big.Int) { 

   for b := br + 1; b <= branches; b++ {
       sum += n
       if sum > nMax {
           return
       }
       if l*2 >= sum && b >= branches {
           return
       }
       if b == br+1 {
           c[br].Mul(&rooted[n], cnt)
       } else {
           tmp.Add(&rooted[n], tmp.SetInt64(int64(b-br-1)))
           c[br].Mul(&c[br], tmp)
           c[br].Div(&c[br], tmp.SetInt64(int64(b-br)))
       }
       if l*2 < sum {
           unrooted[sum].Add(&unrooted[sum], &c[br])
       }
       if b < branches {
           rooted[sum].Add(&rooted[sum], &c[br])
       }
       for m := n - 1; m > 0; m-- {
           tree(b, m, l, sum, &c[br])
       }
   }

}

func bicenter(s int) { 

   if s&1 == 0 {
       tmp.Rsh(tmp.Mul(&rooted[s/2], tmp.Add(&rooted[s/2], one)), 1)
       unrooted[s].Add(&unrooted[s], tmp)
   }

}

func main() { 

   rooted[0].SetInt64(1)
   rooted[1].SetInt64(1)
   unrooted[0].SetInt64(1)
   unrooted[1].SetInt64(1)
   for n := 1; n <= nMax; n++ {
       tree(0, n, n, 1, big.NewInt(1))
       bicenter(n)
       fmt.Printf("%d: %d\n", n, &unrooted[n])
   }

}</lang> Output: (trimmed) 

1: 1
2: 1
3: 1
4: 2
5: 3
6: 5
7: 9
8: 18
9: 35
10: 75
11: 159
12: 355
13: 802
14: 1858
15: 4347
16: 10359
17: 24894
18: 60523
19: 148284
20: 366319
21: 910726
22: 2278658
23: 5731580
24: 14490245
25: 36797588
26: 93839412
27: 240215803
28: 617105614
29: 1590507121
30: 4111846763
31: 10660307791
32: 27711253769
33: 72214088660
34: 188626236139
35: 493782952902
...
499: 2494778503569518415974916458171385008434293725770646826239156067
242244971575380601975205121412154447546983792046199118690890088809311
353933980135961153586789990128720078992516401152141409189797471635981
1774035382416687
500: 6988880697796831865200756887232350988314555919591933763737656855
676489655063916537419419883456606489718804407295350416444314716722189
676957495404000563480243902403405244758705372196760892838840473514056
5376686372508743

Haskell

Using formula from OEIS page, similar to the Mathematica entry below: <lang haskell>-- polynomial utils a `nmul` n = map (*n) a a `ndiv` n = map (`div` n) a

instance (Integral a) => Num [a] where 

 (+) = zipWith (+)
 negate = map negate
 a * b = foldr f undefined b where
   f x z = (a `nmul` x) + (0 : z)
 abs _ = undefined
 signum _ = undefined
 fromInteger n = fromInteger n : repeat 0

-- replace x in polynomial with x^n repl a n = concatMap (: replicate (n-1) 0) a

-- S2: (a^2 + b)/2 cycleIndexS2 a b = (a*a + b)`ndiv` 2

-- S4: (a^4 + 6 a^2 b + 8 a c + 3 b^2 + 6 d) / 24 cycleIndexS4 a b c d = ((a ^ 4) + (a ^ 2 * b) `nmul` 6 + (a * c) `nmul` 8 + (b ^ 2) `nmul` 3 + d `nmul` 6) `ndiv` 24


a598 = x1 -- A000598: A(x) = 1 + (1/6)*x*(A(x)^3 + 3*A(x)*A(x^2) + 2*A(x^3)) x1 = 1 : ((x1^3) + ((x2*x1)`nmul` 3) + (x3`nmul`2)) `ndiv` 6 x2 = x1`repl`2 x3 = x1`repl`3 x4 = x1`repl`4

-- A000678 = x CycleIndex(S4, A000598(x)) a678 = 0 : cycleIndexS4 x1 x2 x3 x4

-- A000599 = CycleIndex(S2, A000598(x) - 1) a599 = cycleIndexS2 (0 : tail x1) (0 : tail x2)

-- A000602 = A000678(x) - A000599(x) + A000599(x^2) a602 = a678 - a599 + x2

main = mapM_ print $ take 200 $ zip [0 ..] a602</lang>

Counting trees with some fairly primitive caching: <lang haskell>import Data.Array

choose :: Integer -> Int -> Integer choose m k = let kk = toInteger k in (product [m..m+kk-1]) `div` (product [1..kk])

max_branches = 4 max_nodes = 200

bcache = listArray (0, max_nodes) [sum[rcache!n!b!r | r <- [0..n], b <- [0..max_branches-1]] | n <- [0..max_nodes]] build_block = (bcache !)

rcache = listArray (0,max_nodes) [arr_b i | i <- [0..max_nodes]] where arr_b n = listArray(0,max_branches) [arr_r b n | b <- [0..max_branches]] arr_r b n = listArray(0,n) [rooted n b r | r <- [0..n]]

rooted 1 0 0 = 1 rooted 1 _ _ = 0 rooted _ 0 _ = 0 rooted _ _ 0 = 0 rooted n b r | (n <= b) || (n <= r) = 0 | otherwise = sum [(firsts b1) * (rests b1) | b1 <- [1..b], r * b1 < n] where firsts = choose (build_block r) rests bb = sum [rcache!(n-r*bb)!(b - bb)!r1 | r1 <- [0..r-1], r1 < (n-r*bb)]

unrooted n = unicenter + bycenter where unicenter = sum [ rcache!n!b!r | b <- [0..max_branches], r <-[0..n], r * 2 < n] bycenter| odd n = 0 | otherwise = x * (x + 1) `div` 2 where x = build_block (n `div` 2)

main = mapM_ print $ map (\x->(x, unrooted x)) [1..max_nodes]</lang>

Output:
(1,1)
(2,1)
(3,1)
(4,2)
(5,3)
(6,5)
(7,9)
(8,18)
(9,35)
(10,75)
(11,159)
(12,355)
(13,802)
(14,1858)
...
(199,339176261988518728096836182493660862745709169352281541101577697702699073887422989905)
(200,943043328799038505167332910595466006794464252841664581909549826351576307818857723954)

J

The following code is an interpretation of the Haskell program listed in the links above.

<lang j>part3=: ;@((<@([(],.(-+/"1))],.]+i.@(]-~1+<.@-:@-))"0 i.@>:@<.@%&3))

part4=: 3 :0 ij=.; (,.]+i.@:(]-~1+[:<.3%~y-]))&.> i.1+<.y%4 (,.y - +/"1) ; (<@(],"1 0 <.@-:@(y-[) (] + i.@>:@-) {:@] >. (>.-:y)-[)~+/)"1 ij )

c0=: */@:{ c1=: 13 :'(*-:@(*>:))/y{~}:x' c2=: 13 :'(*-:@(*>:))~/y{~}.x' c3=: 13 :'3!2+y{~{.x'

radGenN=: [:;[:(],[:+/c0`c1`c2`c3@.(#.@(}.=}:)@[)"1)&.>/(<1x),~part3&.>@ i.@-

bcpGenN=: [: , 0 ,.~ -:@(*>:)@({~i.)

c11=: 13 :'*/(y{~0 1{x), -:(*>:)y{~{:x' c12=: 13 :'*/(y{~0 3{x), -:(*>:)y{~2{x' c13=: 13 :'*/(y{~{.x) , 3!2+ y{~{: x' c14=: 13 :'*/(y{~_2{.x), -:(*>:)y{~{.x' c15=: 13 :'*/ -:(*>:) y{~0 3{x' c16=: 13 :'*/(y{~{:x) , 3!2+ y{~{. x' c17=: 13 :'4!3+y{~{.x'

cassl=: c0`c11`c12`c13`c14`c15`c16`c17

ccpGenN=: 4 :0 if. 0=y do. i.0 return. end. y{.2({.,0,}.) 0,+/@:(x cassl@.(#.@(}.=}:)@[)"1~[)@:part4"0 [1-.~i.y-1 )

NofParaff=: {. radGenN ((ccpGenN +:) + bcpGenN ) 2&|+<.@-:</lang> Output: <lang j> 6 6 $ NofParaff 36 

        1           1           1           1            2            3
        5           9          18          35           75          159
      355         802        1858        4347        10359        24894
    60523      148284      366319      910726      2278658      5731580
 14490245    36797588    93839412   240215803    617105614   1590507121

4111846763 10660307791 27711253769 72214088660 188626236139 493782952902</lang>

Java

Translation of C via Go and D <lang java>import java.math.BigInteger; import java.util.Arrays;

class Test { 

   final static int nMax = 250;
   final static int nBranches = 4;

   static BigInteger[] rooted = new BigInteger[nMax + 1];
   static BigInteger[] unrooted = new BigInteger[nMax + 1];
   static BigInteger[] c = new BigInteger[nBranches];

   static void tree(int br, int n, int l, int inSum, BigInteger cnt) {
       int sum = inSum;
       for (int b = br + 1; b <= nBranches; b++) {
           sum += n;

           if (sum > nMax || (l * 2 >= sum && b >= nBranches))
               return;

           BigInteger tmp = rooted[n];
           if (b == br + 1) {
               c[br] = tmp.multiply(cnt);
           } else {
               c[br] = c[br].multiply(tmp.add(BigInteger.valueOf(b - br - 1)));
               c[br] = c[br].divide(BigInteger.valueOf(b - br));
           }

           if (l * 2 < sum)
               unrooted[sum] = unrooted[sum].add(c[br]);

           if (b < nBranches)
               rooted[sum] = rooted[sum].add(c[br]);

           for (int m = n - 1; m > 0; m--)
               tree(b, m, l, sum, c[br]);
       }
   }

   static void bicenter(int s) {
       if ((s & 1) == 0) {
           BigInteger tmp = rooted[s / 2];
           tmp = tmp.add(BigInteger.ONE).multiply(rooted[s / 2]);
           unrooted[s] = unrooted[s].add(tmp.shiftRight(1));
       }
   }

   public static void main(String[] args) {
       Arrays.fill(rooted, BigInteger.ZERO);
       Arrays.fill(unrooted, BigInteger.ZERO);
       rooted[0] = rooted[1] = BigInteger.ONE;
       unrooted[0] = unrooted[1] = BigInteger.ONE;

       for (int n = 1; n <= nMax; n++) {
           tree(0, n, n, 1, BigInteger.ONE);
           bicenter(n);
           System.out.printf("%d: %s%n", n, unrooted[n]);
       }
   }

}</lang> 

1: 1
2: 1
3: 1
4: 2
5: 3
6: 5
7: 9
8: 18
9: 35
10: 75
(...)
248: 2086006351917005252913566124773054331962205157167696706926185063169623907656246841866717933958839366769700
249: 5814271898167303040368103945830220447130073898083466852225709084407144308593691069932064987528870826155297
250: 16206624309085062837751018464745815688226709117091506494175397665527493805947344857313038875654104100026504

jq

Works with: jq version 1.4

The following version is based on the C/python/ruby implementations.

Currently jq uses IEEE 754 64-bit numbers. Large integers are approximated by floats and therefore the results generated by the program presented here are only precise for n up to and including 45.

<lang jq>def MAX_N: 500; # imprecision begins at 46 def BRANCH: 4;

  1. state: [unrooted, ra]
  2. tree(br; n; l; sum; cnt) where initially: l=n, sum=1 and cnt=1

def tree(br; n; l; sum; cnt): 

 # The inner function is used to implement the range(b+1; BRANCH) loop
 # as there are early exits.
 # On completion, _tree returns [unrooted, ra]
 def _tree: # state [ (b, c, sum),  (unrooted, ra)]
   if length != 5 then error("_tree input has length \(length)") else . end
   | .[0] as $b | .[1] as $c | .[2] as $sum | .[3] as $unrooted | .[4] as $ra
   | if $b > BRANCH then [$unrooted, $ra]
     else
       ($sum + n) as $sum
       | if $sum >= MAX_N or
            # prevent unneeded long math
            ( l * 2 >= $sum and $b >= BRANCH) then [$unrooted, $ra]                      # return
         else (if $b == br + 1 then $ra[n] * cnt 
               else ($c * ($ra[n] + (($b - br - 1)))) / ($b - br) | floor
               end) as $c
         | (if l * 2 < $sum then ($unrooted | .[$sum] += $c)
            else $unrooted end) as $unrooted
         | if $b >= BRANCH then [$b+1, $c, $sum, $unrooted, $ra] | _tree                 # next
           else  [$unrooted, ($ra | .[$sum] += $c) ]
           | reduce range(1; n) as $m (.;  tree($b; $m; l; $sum; $c)) 
           | ([$b + 1, $c, $sum] + .) | _tree 
           end
         end
     end
 ;

 # start by incrementing b, and prepending values for (b,c,sum)
 ([br+1, cnt, sum] + .)  | _tree
  1. input and output: [unrooted, ra]

def bicenter(s): 

 if s % 2 == 1 then .
 else
    .[1][s / 2] as $aux
    | .[0][s] += ($aux * ($aux + 1)) / 2 # 2 divides odd*even
 end

def array(n;init): [][n-1] = init | map(init);

def ra: array( MAX_N; 0) | .[0] = 1 | .[1] = 1;

def unrooted: ra;

  1. See below for a simpler implementation using "foreach"

def paraffins: 

 # range(1; MAX_N)
 def _paraffins(n):
   if n >= MAX_N then empty
   else tree(0; n; n; 1; 1) | bicenter(n)
   | [n, .[0][n]],  # output
     _paraffins(n+1)
   end;
 [unrooted, ra] | _paraffins(1)

paraffins</lang> Output (trimmed): 

$ jq -M -n -c -f paraffins.jq
[1,1]
[2,1]
[3,1]
[4,2]
[5,3]
[6,5]
[7,9]
[8,18]
[9,35]
[10,75]
[11,159]
[12,355]
[13,802]
[14,1858]
[15,4347]
[16,10359]
[17,24894]
[18,60523]
[19,148284]
[20,366319]
[21,910726]
[22,2278658]
[23,5731580]
[24,14490245]
[25,36797588]
[26,93839412]
[27,240215803]
[28,617105614]
[29,1590507121]
[30,4111846763]
[31,10660307791]
[32,27711253769]
[33,72214088660]
[34,188626236139]
[35,493782952902]
[36,1295297588128]
[37,3404490780161]
[38,8964747474595]
[39,23647478933969]
[40,62481801147341]
[41,165351455535782]
[42,438242894769226]
[43,1163169707886427]
[44,3091461011836856]
[45,8227162372221203]
# The above answer for 45 is the last precisely correct answer -- floating point approximations hereon:
...
[100,5.921072038125814e+39]
...
[499,2.4947785035695037e+217]

Using foreach

Works with: jq version >1.4

The following is a more elegant alternative to paraffins/0 as defined above but requires "foreach": <lang jq>def paraffins: 

 foreach range(1; MAX_N) as $n
   ( [unrooted, ra];
     tree(0; $n; $n; 1; 1) | bicenter($n);
     [$n, .[0][$n]]
   )
</lang>

Mathematica

Works with: Mathematica version 9.0

Using the formula on OEIS. <lang mathematica>G000602[n_] := 

 Block[{x}, 
  x*CycleIndexPolynomial[SymmetricGroup[4], 
       Table[ComposeSeries[#, x^i + O[x]^(n + 1)], {i, 4}]] - 
     CycleIndexPolynomial[SymmetricGroup[2], 
      Table[ComposeSeries[# - 1, x^i + O[x]^(n + 1)], {i, 2}]] + 
     ComposeSeries[#, x^2 + O[x]^(n + 1)] &@
   Fold[Series[
      1 + x/6 (#1^3 + 3 #1 ComposeSeries[#1, x^2 + O[x]^#2] + 
          2 ComposeSeries[#1, x^3 + O[x]^#2]), {x, 0, #2}] &, 
    1 + O[x], Range[n + 1]]];

A000602[n_] := SeriesCoefficient[G000602[n], n]; A000602List[n_] := CoefficientList[G000602[n], x]; Grid@Transpose@{Range[0, 200], A000602List@200}</lang>

Output:
0	1
1	1
2	1
3	1
4	2
5	3
6	5
7	9
8	18
9	35
10	75
11	159
12	355
13	802
...
199	339176261988518728096836182493660862745709169352281541101577697702699073887422989905
200	943043328799038505167332910595466006794464252841664581909549826351576307818857723954

Nim

Translation of: C

<lang nim>import bigints

const 

 nMax: int32 = 250
 nBranches = 4

var rooted, unrooted: array[nMax + 1, BigInt] rooted[0..1] = [1.initBigInt, 1.initBigInt] unrooted[0..1] = [1.initBigInt, 1.initBigInt] for i in 2 .. nMax: 

 rooted[i] = 0.initBigInt
 unrooted[i] = 0.initBigInt

proc choose(m, k): BigInt = 

 result = m
 if k == 1: return
 for i in 1 .. < k:
   result = result * (m + i) div (i + 1)

proc tree(br, n, l, sum: int32, cnt: BigInt) = 

 var s: int32 = 0
 for b in br + 1 .. nBranches:
   s = sum + (b - br) * n
   if s > nMax: return

   let c = choose(rooted[n], b - br) * cnt

   if l * 2 < s: unrooted[s] += c
   if b == nBranches: return
   rooted[s] += c
   for m in countdown(n-1, 1):
     tree b, m, l, s, c

proc bicenter(s) = 

 if (s and 1) == 0:
   unrooted[s] += rooted[s div 2] * (rooted[s div 2] + 1) div 2

for n in 1 .. nMax: 

 tree 0, n, n, 1, 1.initBigInt
 n.bicenter
 echo n, ": ", unrooted[n]</lang>

Output: 

1: 1
2: 1
3: 1
4: 2
5: 3
6: 5
7: 9
8: 18
9: 35
10: 75
11: 159
12: 355
13: 802
14: 1858
15: 4347
16: 10359
17: 24894
18: 60523
19: 148284
20: 366319
21: 910726
...
249: 5814271898167303040368103945830220447130073898083466852225709084407144308593691069932064987528870826155297
250: 16206624309085062837751018464745815688226709117091506494175397665527493805947344857313038875654104100026504

PARI/GP

This function is for recent PARI/GP: <lang parigp>paraffin(p) = { 

 local (P = p+1, R, U = R = Vec([1,1], P));

 for (n = 1, p,
   ((B,n,C,S,l=n) -> my(b,c,i,s);
     for (b = 1, 4-B,
       if ((s = S + b * n) < P,
         c = R[n+1] * C * prod(i = 1, b-1, (R[n+1]+i) / (i+1)) ;
         if (l+l < s, U[s+1] += c);
         if (B+b < 4, R[s+1] += c; i = n; while (i--, self()(B+b, i, c, s, l)))))
   )(0,n,1,1);
   if (n % 2,, U[n+1] += R[n/2+1] * (R[n/2+1]+1) / 2);
   print([n, U[n+1]]))

}</lang>


Code for older version of PARI/GP < 2.9: <lang parigp>iso(B,n,C,S,l=n) = { 

 my(b,c,i,s);

 for (b = 1, 4-B,
   if ((s = S + b * n) < P,
     c = R[n+1] * C * prod(i = 1, b-1, (R[n+1]+i) / (i+1));
     if (l+l < s, U[s+1] += c);
     if (B+b < 4, R[s+1] += c; i = n; while (i--, iso(B+b, i, c, s, l)))))

}

paraffin(p) = { 

 local (P = p+1, R, U = R = Vec([1,1], P));

 for (n = 1, p, iso(0, n, 1, 1);
   if (n % 2,, U[n+1] += R[n/2+1] * (R[n/2+1]+1) / 2);
   print([n, U[n+1]]))

}</lang>


Output:

paraffin(32)

[1, 1]
[2, 1]
[3, 1]
[4, 2]
[5, 3]
[6, 5]
[7, 9]
[8, 18]
[9, 35]
[10, 75]
[11, 159]
[12, 355]
[13, 802]
[14, 1858]
[15, 4347]
[16, 10359]
[17, 24894]
[18, 60523]
[19, 148284]
[20, 366319]
[21, 910726]
[22, 2278658]
[23, 5731580]
[24, 14490245]
[25, 36797588]
[26, 93839412]
[27, 240215803]
[28, 617105614]
[29, 1590507121]
[30, 4111846763]
[31, 10660307791]
[32, 27711253769]

Pascal

Works with: Free_Pascal
Library: GMP

Conversion of the C example: <lang pascal>Program Paraffins;

uses 

 gmp;

const 

 max_n = 500;
 branch = 4;

var 

 rooted, unrooted: array [0 .. max_n-1] of mpz_t;
 c: array [0 .. branch-1] of mpz_t;
 cnt, tmp: mpz_t;
 n: integer;
 fmt: pchar;
 sum: integer;

procedure tree(br, n, l: integer; sum: integer; cnt: mpz_t); 

 var
   b, m: integer;
 begin
   for b := br + 1 to branch do
   begin
     sum := sum + n;
     if sum >= max_n then

exit; 

     (* prevent unneeded long math *)
     if (l * 2 >= sum) and (b >= branch) then

exit; 

     if b = (br + 1) then

mpz_mul(c[br], rooted[n], cnt) 

     else
     begin

mpz_add_ui(tmp, rooted[n], b - br - 1); mpz_mul(c[br], c[br], tmp); mpz_divexact_ui(c[br], c[br], b - br); 

     end;

     if l * 2 < sum then

mpz_add(unrooted[sum], unrooted[sum], c[br]); 

     if b < branch then
     begin

mpz_add(rooted[sum], rooted[sum], c[br]); for m := n-1 downto 1 do tree(b, m, l, sum, c[br]); 

     end;
   end;
 end;

procedure bicenter(s: integer); begin 

 if odd(s) then
   exit;
 mpz_add_ui(tmp, rooted[s div 2], 1);
 mpz_mul(tmp, rooted[s div 2], tmp);
 mpz_tdiv_q_2exp(tmp, tmp, 1);

 mpz_add(unrooted[s], unrooted[s], tmp);

end;

begin 

 for n := 0 to 1 do
 begin
   mpz_init_set_ui(rooted[n], 1);
   mpz_init_set_ui(unrooted[n], 1);
 end;
 for n := 2 to max_n-1 do
 begin
   mpz_init_set_ui(rooted[n], 0);
   mpz_init_set_ui(unrooted[n], 0);
 end;
 for n := 0 to BRANCH-1 do
   mpz_init(c[n]);

 mpz_init(tmp);

 mpz_init_set_ui(cnt, 1);
 sum := 1;
 for n := 1 to MAX_N do
 begin
   tree(0, n, n, sum, cnt);
   bicenter(n);
   mp_printf('%d: %Zd'+chr(13)+chr(10), n, @unrooted[n]);
 end;

end.</lang> Output (trimmed): 

1: 1
2: 1
3: 1
4: 2
5: 3
6: 5
7: 9
8: 18
9: 35
10: 75
11: 159
12: 355
13: 802
14: 1858
15: 4347
16: 10359
17: 24894
18: 60523
19: 148284
20: 366319
21: 910726
22: 2278658
23: 5731580
24: 14490245
25: 36797588
26: 93839412
27: 240215803
28: 617105614
29: 1590507121
30: 4111846763
31: 10660307791
32: 27711253769
33: 72214088660
34: 188626236139
35: 493782952902
..
499: 2494778503569518415974916458171385008434293725770646826239156067242244971575380601975205121412154447546983792
     0461991186908900888093113539339801359611535867899901287200789925164011521414091897974716359811774035382416687
500: 35027241765694350953134643689731689510165490186734465975716452976264824623941928748229034890456068034873862265
     02245470029705798383457667757879016306516643706616463708365586217812176623950823312143595975577679489832553344

Perl

Translation of: C

This is using Math::GMPz for best performance. Math::GMP works almost as well. Math::BigInt is in core and only 9 times slower for this task. <lang perl>use Math::GMPz;

my $nmax = 250; my $nbranches = 4;

my @rooted = map { Math::GMPz->new($_) } 1,1,(0) x $nmax; my @unrooted = map { Math::GMPz->new($_) } 1,1,(0) x $nmax; my @c = map { Math::GMPz->new(0) } 0 .. $nbranches-1;

sub tree { 

 my($br, $n, $l, $sum, $cnt) = @_;
 for my $b ($br+1 .. $nbranches) {
   $sum += $n;
   return if $sum > $nmax || ($l*2 >= $sum && $b >= $nbranches);
   if ($b == $br+1) {
     $c[$br] = $rooted[$n] * $cnt;
   } else {
     $c[$br] *= $rooted[$n] + $b - $br - 1;
     $c[$br] /= $b - $br;
   }
   $unrooted[$sum] += $c[$br] if $l*2 < $sum;
   return if $b >= $nbranches;
   $rooted[$sum] += $c[$br];
   for my $m (reverse 1 .. $n-1) {
     next if $sum+$m > $nmax;
     tree($b, $m, $l, $sum, $c[$br]);
   }
 }

}

sub bicenter { 

 my $s = shift;
 $unrooted[$s] += $rooted[$s/2] * ($rooted[$s/2]+1) / 2  unless $s & 1;

}

for my $n (1 .. $nmax) { 

 tree(0, $n, $n, 1, Math::GMPz->new(1));
 bicenter($n);
 print "$n: $unrooted[$n]\n";

}</lang>

Output identical to C GMP example (truncated to 250).

Perl 6

Works with: rakudo version 2016.04

Counting only, same algorithm as the C solution with some refactorings.

Note how lexical scoping — rather than global variables or repeated arguments — is used to pass down information to subroutines.

<lang perl6>sub count-unrooted-trees(Int $max-branches, Int $max-weight) { 

   my @rooted   = flat 1,1,0 xx $max-weight - 1;
   my @unrooted = flat 1,1,0 xx $max-weight - 1;

   sub count-trees-with-centroid(Int $radius) {
       sub add-branches(
           Int $branches,        # number of branches to add
           Int $w,               # weight of heaviest branch to add
           Int $weight  is copy, # accumulated weight of tree
           Int $choices is copy, # number of choices so far
       ) {
           $choices *= @rooted[$w];
           for 1 .. $branches -> $b {
               ($weight += $w) <= $max-weight or last;
               @unrooted[$weight] += $choices if $weight > 2*$radius;
               if $b < $branches {
                   @rooted[$weight] += $choices;
                   add-branches($branches - $b, $_, $weight, $choices) for 1 ..^ $w;
                   $choices = $choices * (@rooted[$w] + $b) div ($b + 1);
               }
           }
       }
       add-branches($max-branches, $radius, 1, 1);
   }

   sub count-trees-with-bicentroid(Int $weight) {
       if $weight %% 2 {
           my \halfs = @rooted[$weight div 2];
           @unrooted[$weight] += (halfs * (halfs + 1)) div 2;
       }
   }

   gather {
       take 1;
       for 1 .. $max-weight {
           count-trees-with-centroid($_);
           count-trees-with-bicentroid($_);
           take @unrooted[$_];
       }
   }

}

my constant N = 100; my @paraffins = count-unrooted-trees(4, N); say .fmt('%3d'), ': ', @paraffins[$_] for flat 1 .. 30, N;</lang>

Output:
  1: 1
  2: 1
  3: 1
  4: 2
  5: 3
  6: 5
  7: 9
  8: 18
  9: 35
 10: 75
 11: 159
 12: 355
 13: 802
 14: 1858
 15: 4347
 16: 10359
 17: 24894
 18: 60523
 19: 148284
 20: 366319
 21: 910726
 22: 2278658
 23: 5731580
 24: 14490245
 25: 36797588
 26: 93839412
 27: 240215803
 28: 617105614
 29: 1590507121
 30: 4111846763
100: 5921072038125809849884993369103538010139

Pike

Translation of: Python

<lang Pike>int MAX_N = 300; int BRANCH = 4;

array ra = allocate(MAX_N); array unrooted = allocate(MAX_N);

void tree(int br, int n, int l, int sum, int cnt) { 

   int c;
   for (int b = br + 1; b < BRANCH + 1; b++)
   {
       sum += n;
       if (sum >= MAX_N)
           return;

       // prevent unneeded long math
       if (l * 2 >= sum && b >= BRANCH)
           return;

       if (b == br + 1)
       {
           c = ra[n] * cnt;
       }
       else
       {
           c = c * (ra[n] + (b - br - 1)) / (b - br);
       }

       if (l * 2 < sum)
           unrooted[sum] += c;

       if (b < BRANCH)
       {
           ra[sum] += c;
           for (int m=1; m < n; m++)
           {
               tree(b, m, l, sum, c);
           }
       }
   }

}

void bicenter(int s) { 

   if (!(s & 1))
   {
       int aux = ra[s / 2];
       unrooted[s] += aux * (aux + 1) / 2;
   }

}


void main() { 

   ra[0] = ra[1] = unrooted[0] = unrooted[1] = 1;

   for (int n = 1; n < MAX_N; n++)
   {
       tree(0, n, n, 1, 1);
       bicenter(n);
       write("%d: %d\n", n, unrooted[n]);
   }

}</lang>

Python

This version only counts different paraffins. The multi-precision integers of Python avoid overflows.

Translation of: C

<lang python>try: 

   import psyco
   psyco.full()

except ImportError: 

   pass

MAX_N = 300 BRANCH = 4

ra = [0] * MAX_N unrooted = [0] * MAX_N

def tree(br, n, l, sum = 1, cnt = 1): 

   global ra, unrooted, MAX_N, BRANCH
   for b in xrange(br + 1, BRANCH + 1):
       sum += n
       if sum >= MAX_N:
           return

       # prevent unneeded long math
       if l * 2 >= sum and b >= BRANCH:
           return

       if b == br + 1:
           c = ra[n] * cnt
       else:
           c = c * (ra[n] + (b - br - 1)) / (b - br)

       if l * 2 < sum:
           unrooted[sum] += c

       if b < BRANCH:
           ra[sum] += c;
           for m in range(1, n):
               tree(b, m, l, sum, c)

def bicenter(s): 

   global ra, unrooted
   if not (s & 1):
       aux = ra[s / 2]
       unrooted[s] += aux * (aux + 1) / 2


def main(): 

   global ra, unrooted, MAX_N
   ra[0] = ra[1] = unrooted[0] = unrooted[1] = 1

   for n in xrange(1, MAX_N):
       tree(0, n, n)
       bicenter(n)
       print "%d: %d" % (n, unrooted[n])

main()</lang> Output (newlines added): 

1: 1
2: 1
3: 1
4: 2
5: 3
6: 5
7: 9
8: 18
9: 35
10: 75
11: 159
12: 355
13: 802
14: 1858
15: 4347
16: 10359
17: 24894
18: 60523
19: 148284
20: 366319
21: 910726
22: 2278658
23: 5731580
24: 14490245
25: 36797588
26: 93839412
27: 240215803
28: 617105614
29: 1590507121
30: 4111846763
31: 10660307791
32: 27711253769
33: 72214088660
34: 188626236139
35: 493782952902
36: 1295297588128
...
498: 8905549466780876642396073343654258905917693028466999354614122414
334803973053490695667380301725459866905052509157458138657523790806693
604983304945043446074074823147488033014740909049433466213254508146498
309733349188106
499: 2494778503569518415974916458171385008434293725770646826239156067
242244971575380601975205121412154447546983792046199118690890088809311
353933980135961153586789990128720078992516401152141409189797471635981
1774035382416687

Using generating function

This is almost the same as the one in Formal power series. Compare to the Mathematica and Haskell solutions. <lang python>from itertools import count, chain, tee, islice, cycle from fractions import Fraction from sys import setrecursionlimit setrecursionlimit(5000)

def frac(a,b): return a//b if a%b == 0 else Fraction(a,b)

  1. infinite polynomial class

class Poly: 

   def __init__(self, gen = None):
       self.gen, self.source = (None, gen) if type(gen) is Poly \
           else (gen, None)

   def __iter__(self):
       # We're essentially tee'ing it everytime the iterator
       # is, well, iterated.  This may be excessive.
       return Poly(self)

   def getsource(self):
       if self.gen == None:
           s = self.source
           s.getsource()
           s.gen, self.gen = tee(s.gen, 2)

   def next(self):
       self.getsource()
       return next(self.gen)

   __next__ = next

   # Overload "<<" as stream input operator. Hey, C++ does it.
   def __lshift__(self, a): self.gen = a

   # The other operators are pretty much what one would expect
   def __neg__(self): return Poly(-x for x in self)

   def __sub__(a, b): return a + (-b)

   def __rsub__(a, n):
       a = Poly(a)
       def gen():
           yield(n - next(a))
           for x in a: yield(-x)
       return Poly(gen())

   def __add__(a, b):
       if type(b) is Poly:
           return Poly(x + y for (x,y) in zip(a,b))

       a = Poly(a)
       def gen():
           yield(next(a) + b)
           for x in a: yield(x)

       return Poly(gen())

   def __radd__(a,b):
       return a + b

   def __mul__(a,b):
       if not type(b) is Poly:
           return Poly(x*b for x in a)

       def gen():
           s = Poly(cycle([0]))
           for y in b:
               s += y*a
               yield(next(s))

       return Poly(gen())

   def __rmul__(a,b): return a*b

   def __truediv__(a,b):
       if not type(b) is Poly:
           return Poly(frac(x, b) for x in a)

       a, b = Poly(a), Poly(b)
       def gen():
           r, bb = a,next(b)
           while True:
               aa = next(r)
               q = frac(aa, bb)
               yield(q)
               r -= q*b

       return Poly(gen())

   def repl(self, n):
       def gen():
           for x in self:
               yield(x)
               for i in range(n-1): yield(0)
       return Poly(gen())

   def __pow__(self, n):
       return Poly(self) if n == 1 else self * self**(n-1)

def S2(a,b): return (a*a + b)/2 def S4(a,b,c,d): return a**4/24 + a**2*b/4 + a*c/3 + b**2/8 + d/4

x1 = Poly() x2 = x1.repl(2) x3 = x1.repl(3) x4 = x1.repl(4) x1 << chain([1], (x1**3 + 3*x1*x2 + 2*x3)/6)

a598 = x1 a678 = Poly(chain([0], S4(x1, x2, x3, x4))) a599 = S2(x1 - 1, x2 - 1) a602 = a678 - a599 + x2

for n,x in zip(count(0), islice(a602, 500)): print(n,x)</lang>

Racket

This Scheme solution runs in Racket too:

Or, a direct translation of the C entry: <lang racket>

  1. lang racket

(define MAX_N 33) (define BRANCH 4)

(define rooted (make-vector MAX_N 0)) (define unrooted (make-vector MAX_N 0)) (for ([i 2]) (vector-set! rooted i 1) (vector-set! unrooted i 1))

(define (vector-inc! v i d) (vector-set! v i (+ d (vector-ref v i))))

(define (choose m k) 

 (if (= k 1) m
     (for/fold ([r m]) ([i (in-range 1 k)]) (/ (* r (+ m i)) (add1 i)))))

(define (tree br n cnt sum l) 

 (let/ec return
   (for ([b (in-range (add1 br) (add1 BRANCH))])
     (define s (+ sum (* (- b br) n)))
     (when (>= s MAX_N) (return))
     (define c (* (choose (vector-ref rooted n) (- b br)) cnt))
     (when (< (* l 2) s) (vector-inc! unrooted s c))
     (when (= b BRANCH) (return))
     (vector-inc! rooted s c)
     (for ([m (in-range (sub1 n) 0 -1)]) (tree b m c s l)))))

(define (bicenter s) 

 (when (even? s)
   (vector-inc! unrooted s (* (vector-ref rooted (/ s 2))
                              (add1 (vector-ref rooted (/ s 2)))
                              1/2))))

(for ([n (in-range 1 MAX_N)]) 

 (tree 0 n 1 1 n)
 (bicenter n)
 (printf "~a: ~a\n" n (vector-ref unrooted n)))

</lang>

REXX

(Based, in large part, on the Pascal version.) <lang rexx>/*REXX program enumerates (without repetition) the # of paraffins with N atoms of carbon*/ parse arg nodes . /*obtain optional argument from the CL.*/ if nodes== | nodes=="," then nodes=100 /*Not specified? Then use the default.*/ 

 rooted. =  0;      rooted.0=1;     rooted.1=1  /*define the  base     rooted  numbers.*/

unrooted. = 0; unrooted.0=1; unrooted.1=1 /* " " " unrooted " */ numeric digits max(9,nodes%2) /*this program may use gihugeic numbers*/ w=length(nodes) /*W: used for aligning formatted nodes.*/ say right(0,w) unrooted.0 /*show enumerations of 0 carbon atoms*/ 

                                                /* [↓]  process all nodes (up to NODES)*/
         do C=1  for nodes;         h=C%2       /*C:   is the number of carbon atoms.  */
         call  tree  0, C, C, 1, 1              /* [↓]  if # of carbon atoms is even···*/
         if C//2==0  then unrooted.C=unrooted.C   +   rooted.h * (rooted.h + 1) % 2
         say right(C,w)  unrooted.C             /*display an aligned formatted number. */
         end   /*C*/

exit /*stick a fork in it, we're all done. */ /*──────────────────────────────────────────────────────────────────────────────────────*/ tree: procedure expose rooted. unrooted. nodes #. /*this function is recursive.*/ 

     parse arg br,n,L,sum,cnt;  nm=n-1;     LL=L+L;       brp=br+1
              do b=brp  to 4;   sum=sum+n;  if sum>nodes  then leave
              if b==4      then             if LL>=sum    then leave
              if b==brp    then #.br=rooted.n * cnt
                           else #.br=#.br * (rooted.n + b - brp) % (b - br)
              if LL<sum    then unrooted.sum=unrooted.sum + #.br
              if b==4      then leave
              rooted.sum = rooted.sum + #.br
                 do m=nm  by -1  for nm;   call tree b, m, L, sum, #.br;   end  /*m*/
              end   /*b*/               /* ↑↑↑↑↑↑↑↑↑    recursive invocation of  TREE. */</lang>

output   when using the input of:   260 

  0 1
  1 1
  2 1
  3 1
  4 2
  5 3
  6 5
  7 9
  8 18
  9 35
 10 75
 11 159
 12 355
 13 802
 14 1858
 15 4347
 16 10359
 17 24894
 18 60523
 19 148284
 20 366319
 21 910726
 22 2278658
 23 5731580
 24 14490245
 25 36797588
 26 93839412
 27 240215803
 28 617105614
 29 1590507121
 30 4111846763
 31 10660307791
 32 27711253769
 33 72214088660
 34 188626236139
 35 493782952902
 36 1295297588128
 37 3404490780161
 38 8964747474595
 39 23647478933969
 40 62481801147341
 41 165351455535782
 42 438242894769226
 43 1163169707886427
 44 3091461011836856
 45 8227162372221203
 46 21921834086683418
 47 58481806621987010
 48 156192366474590639
 49 417612400765382272
 50 1117743651746953270
 51 2994664179967370611
 52 8031081780535296591
 53 21557771913572630901
 54 57919180873148437753
 55 155745431857549699124
 56 419149571193411829372
 57 1128939578361332867936
 58 3043043571906827182530
 59 8208615366863753915949
 60 22158734535770411074184
 61 59858097847706865855186
 62 161805725349297357221898
 63 437671691526158936922623
 64 1184616185385310843585573
 65 3208285066181475821271463
 66 8694130712024868414002815
 67 23573796134448175745408811
 68 63955159527348138708694312
 69 173603007393950249896865875
 70 471484798515330363034639871
 71 1281151315764638215613845510
 72 3482965749140691245110434511
 73 9473447386804490449091871124
 74 25779306238954404972323916397
 75 70183211512214096492433058105
 76 191156381393249393027319384769
 77 520874195248906781713044332539
 78 1419908915343952137338409797325
 79 3872282575137005474139119076135
 80 10564476906946675106953415600016
 81 28833609436277333169440806135431
 82 78725585464391037293036629979444
 83 215027809474796675607407513633870
 84 587531723826577193455385789266377
 85 1605913778494711520354663202536756
 86 4391002908093323425994602631972445
 87 12010257907756938974208750945664835
 88 32861295558120887536942123568548502
 89 89940959024891576997396491928932689
 90 246245150242821439632304475956113295
 91 674391606297983432514229725117306224
 92 1847515048012613337782670842346319120
 93 5062818112121161180862827915688625902
 94 13877857529584521384324419956411729295
 95 38051836070803837001309074456088423358
 96 104363664561059273927704242814298678658
 97 286312976836850192359345859166390622180
 98 785684759853087702778573182234297830503
 99 2156596319845084996862701478402986311496
100 5921072038125809849884993369103538010139
101 16260750014333666174953055376699249561110
102 44667063168726812052821334495766769690630
103 122726610195426301690448676841677827340780
104 337281538963751874669853952178948219200633
105 927143441542280244466720172757699607129825
106 2549176520305910764377448963173035784835631
107 7010510656300876673813654064741809461787182
108 19283877336110239907079044091958051661009951
109 53055727810105880736027950213934519705620559
110 146002972524313232817393491844985704938385801
111 401865724190508834753025926637435418813476039
112 1106339625432222709435767174129826811545391101
113 3046369875968510015403046201590835240153395100
114 8389999420170754836800638580300552381250693062
115 23111326593011774543116815302964652139347135182
116 63675155467360117136901070528242608498818046250
117 175467195960062612437322237574246321515725845634
118 483616671898832299071277369263305813784565460114
119 1333167312321418940566764416056977442040495550342
120 3675740183950426011078357941139728051663026172228
121 10136322901774027447848977748665383292736169662267
122 27956983197937526275999613945221497078279509595407
123 77121096978813982358935411851692069578533009193138
124 212778592638033483022781655638827961970402357080215
125 587155794584829621068447884048323985962957796104395
126 1620497362318232091081355117667505915417499978679013
127 4473132502312622821884079561929897829404710575328024
128 12349306792492607803837161096610238756912653878568775
129 34098849774876383478036291434385545792965491914980650
130 94167748474814466028838037996326649316233175269577493
131 260093170948828891650104553710684162327855828421145690
132 718487205759724277833835055443909476145495116155508155
133 1985050220521088907210323840127550973214943015739291120
134 5485110653386099899275645856977233423965042141295771502
135 15158624968755754576600389921905653584106659889930620820
136 41898053824932615440265041900412507427220728337249680527
137 115820822448502452349822520317304132018285539473087897141
138 320211802888589798701825810680319271475504997973219083170
139 885411355238188116465394365370757710295372148438998022826
140 2448550585524918609600214967948504177437555812600018440773
141 6772180336728084537425567078328320492989943976904644119200
142 18732796033402227075307540055538834651333956120072379687678
143 51823958523558404531622255138725304201359354985976024954747
144 143387634030485523461662580179416231260007790242619239696168
145 396775836020295064920040342935953476579230225268967120945252
146 1098070975453594757891511609218085888434254839495604326720679
147 3039251105982158526063018965393900608357891201531016545453671
148 8413041613874240075848233530979949485087059914285837491890647
149 23291051500594069758631194545655502320903778728677961917787017
150 64487285324785805685734825467573942213924157583096655274158296
151 178569541961158786360447600422369518262867694211679827307797522
152 494525085028771691070376002671999818542495469214503552543494392
153 1369671107847363840368349527801907625890550280754690871159384167
154 3793941909035282970536126899217159922244816989321008990552343933
155 10510197366726219419291185594221700080820692107164072063617583537
156 29118988780427095392911694296588496006042150251385271606702314123
157 80683801316548713731547508195369620842564695928190934573122040053
158 223583881196691039626561929582827819978293196688915654006262185620
159 619638153674192054430980737083649826660231642062541457264064352590
160 1717428978037773119953669826811378686605009454833429087085211111817
161 4760601845949152288761851253868434642647478893273231351009713026471
162 13197353449186709568929592710369551730672357721629256541119049584662
163 36589226826424289787166608629764201634396432211204018812233928108927
164 101451975263926040804307438557581821336425438886780992752450611791029
165 281324901033788598583154170205263556814795889090791609969956549076553
166 780181677818281299965193432627955631078491817302716837810578348379410
167 2163828038323756757063639945018570904120396578192324738853110253083851
168 6001899167570139611127915072874671685163847392112466633395193150607161
169 16649191065671323727576273232609293462550308875754570697371028619095529
170 46188686972056579073145176280276791118176099297131728121378351698216964
171 128149137125681447665302588425507023489080631426025859378879991574150361
172 355576383032176188837060897590191000068058505378256432541548821711409736
173 986703913063443346422020725722084251185909113284392827422830038792419867
174 2738275183964917202164682060710234556685852044781624370789938274187242387
175 7599818348156354735525837090092498330135165342551619766604085368593605623
176 21094284799140605474267783653778252494175708547669907184929527663028371844
177 58554660677719531288883019197284429180673377561888244491058170393359945984
178 162552183133868639189244204285356619593212307470997836346642760606493409411
179 451292826786530619879633220482642976940485477290114448603416892241141577694
180 1253019870825476025726441067676567248038950763298814178748038046446512128926
181 3479293084378459187212303139960535018989517537846033787292960498791544468857
182 9661781855977284524013799278118239872342899149756408496918889491272019198160
183 26832197158239597797570968340612728947891256166650480273266227097169558934791
184 74522545727244539603451333395337695567614066525612042720018110555143893455632
185 206990881176753531116559188573370889805581604324744059300494333307748123498957
186 574971719221297425559348824161112797452658996937464320320048053830834065076638
187 1597250942564001477500533605167309927398304330031144648098072721512668593957703
188 4437423571982333312534972159110678450135834859468229274326790786916695731276497
189 12328758711422329105019389982539560951228986597668702145097956175468519348920309
190 34256124585721478074980980873512523523896822875906637442595527046990665266761523
191 95189094589104790904556217884090558824685828617516319665565748564984723369457220
192 264524521940855272106937702820301986845470010150944446148610489622177560655580196
193 735146927788110318878638054407335543366855876665936464594690408421993895145574507
194 2043202995476015462049187462882169976289343296164934404442378707876446055277665852
195 5679076882963913929265887525377096781591407289261632655627568444557125995319535956
196 15786016263625679649343179544010857226174369384245579060916786714528288038933116607
197 43882930188633901470015828734437451959746697520345645823789728504039719238235284266
198 121996306076853365751053531202168307620916572983606780123661900035869303555630148063
199 339176261988518728096836182493660862745709169352281541101577697702699073887422989905
200 943043328799038505167332910595466006794464252841664581909549826351576307818857723954
201 2622195090600379263364346956264279702121691087227848066895838284611152725976467138514
202 7291640328972323818932818268921088199080628040707037288217930491456875016135548131376
203 20277391980621940663950418790370236703345679504035237316723532280155012704421841134349
204 56393014827755686247101145333945229562531368138401437274854961321951321148392478837326
205 156842815530515935964014240194651333844507609987619166694876840912879248319077130581042
206 436244327522179535577207667646065280269833187002466982658692809627210486721781255271000
207 1213446271931548955557154166653292946893343739485414159573064525913091105471901382363618
208 3375488708820014134062868343953409434477577207616345786043331438445002602958236463369496
209 9390265533842684145381903993662706957889355090166833985181893569515242461156430420087174
210 26124257322713604151166532772505893583948958402141893219135619374261694917895951421995216
211 72683304203243676344755903584747211387194000236278332965871882088118015556154063840306823
212 202231949421481999766866699910650758032171534187352358158377153614156835833334273717131223
213 562715711666310319461011612553561742990998466225938804278692337156743603187662151365368333
214 1565857565336512188705390960430387447731996020954944548096469264981742312029607101653583428
215 4357517959671123300838959993742696621943017700847213254885626681813005595459240045099075216
216 12126894898610872886310565416845280412742023065548437498642351657251934710541687754081654026
217 33750741717021238734330907104325257824118779361033241504478398185637818043649623474315265399
218 93937737312335248803818931862078208076074752854934195324106075781489944419493627198038389700
219 261468709433838684317888993242737209511093586216774554328802510011968600660643442350819036063
220 727816668798656458204462998706538701005731200304158665308881919140152438032334192310185418095
221 2026033924729657796058178057776349080819410108838068133330231127278657184350600523917658261950
222 5640189120704586237460028096950040726608416653216674308167097595939322372553028535058020450921
223 15702277858615709125768061794209929881772942003228681661684701314937614138722747596565403037124
224 43717315979005745749656846671283395378656317609179742248155824749456112069949537586095666617258
225 121721128292306933059993974247348540921262802966704236768446493370986884003643919966941048534839
226 338922110694788427241777802405838150613492773013920995745534133385509651242705886920398223181010
227 943745948264384101974654246655344635790924945832998571660171925129987029630345274360528181753602
228 2628036396188665255122485857407491640837938331141774537726852980481211265329884171692375404118006
229 7318608979755404166520379853847121693286549981969246032578382667986026415973227568070229589183702
230 20381982477815712032070175127992226329196107434900428011148438268022423947798435659837662375820784
231 56765548015352873669830186155391542225154245344836703122834457861735047936794014876150374629860524
232 158104270172530515145265641139760077405201789758339337585177180751307186891559972317181979411965521
233 440374885613967273764570923053527071581411261510943039612254074378568326443950283692532685610902647
234 1226652343389215664746827690619448544313348534138985205724729790807207201280976341202710717893074853
235 3416963010634291402556279395279280039897379691659516509445008289967225395946292319651344904483142876
236 9518723877111155596143748502002445531218893550853129590490604506008598318134823819368242604272122748
237 26517750521770746045881213675376639453619783582926617680386483834380296365299612479395492961088242736
238 73877802385952947032921713908581985111328188063663582097893741466810319537930797134111900695946616163
239 205830832103408707745581845015281353687041390648973776794197115951467164317467871635600767010770430175
240 573490073391076513109930570190363749116344599275683632274375821562743183409191580097184195332209714304
241 1597939145083038119284545851972427651672212885904578621056479058274066462851925440771984480190480136917
242 4452595736241139504848505768149348728327738257189209919868237125085182422374731188865584214640572836004
243 12407515928191678703401447333122259693080851830690277274671393962029357182157422967017478923028296493797
244 34576004768296889785887066794910718730985852896707505707076422305798138880427343561380451664648670542260
245 96356944442415066997623733664230869603611716312377642117711752082444867783045964116385053282421589905891
246 268540209617944059776303921971316267806082307732727328919911735252505071315985779867844056236080213000010
247 748434113260252449609376828666343341456610378512725586135955779939163320693444008584504390382026391947130
248 2086006351917005252913566124773054331962205157167696706926185063169623907656246841866717933958839366769700
249 5814271898167303040368103945830220447130073898083466852225709084407144308593691069932064987528870826155297
250 16206624309085062837751018464745815688226709117091506494175397665527493805947344857313038875654104100026504
251 45175937003790110836570547902793926244224423679908111750183970641450850884015931131155869896510573436377669
252 125932844251389164587464842843820559638763899540510199408003184309130630867218196313865482109056537730398481
253 351065342330255763232857455169736646100231492324922604469018864206124242515874849723301032460307149363100441
254 978709648454897177175122229161387383812277957992041737463093257271644176592407853706597679764899372036161233
255 2728579524947965774924389584178729664141559400119961753692868837602435434199845226128003794876379357316697077
256 7607396703016294839448365361938795245280010901284162245456747052211826446860068129026466203061153094215518724
257 21210557506627112926709102814428055407062200441680249553023821619098519349155457094641057751120559709171040766
258 59140439362561173332419239067243021493683835575971725019164864794611034386190038342974191464909448990949160433
259 164904810640488270400543694428054069489544003984942619936396942506713662286505814885161827364845226768614358614
260 459831050079703806398893682379638814715447009656839612265184371710342022860900424760438175669232045681912458051

Ruby

Translation of: Python

<lang ruby>MAX_N = 500 BRANCH = 4

def tree(br, n, l=n, sum=1, cnt=1) 

 for b in br+1 .. BRANCH
   sum += n
   return if sum >= MAX_N
   # prevent unneeded long math
   return if l * 2 >= sum and b >= BRANCH
   if b == br + 1
     c = $ra[n] * cnt
   else
     c = c * ($ra[n] + (b - br - 1)) / (b - br)
   end
   $unrooted[sum] += c if l * 2 < sum
   next if b >= BRANCH
   $ra[sum] += c
   (1...n).each {|m| tree(b, m, l, sum, c)}
 end

end

def bicenter(s) 

 return if s.odd?
 aux = $ra[s / 2]
 $unrooted[s] += aux * (aux + 1) / 2

end

$ra = [0] * MAX_N $unrooted = [0] * MAX_N

$ra[0] = $ra[1] = $unrooted[0] = $unrooted[1] = 1 for n in 1...MAX_N 

 tree(0, n)
 bicenter(n)
 puts "%d: %d" % [n, $unrooted[n]]

end</lang>

Output:
1: 1
2: 1
3: 1
4: 2
5: 3
6: 5
7: 9
8: 18
9: 35
10: 75
11: 159
12: 355
13: 802
14: 1858
15: 4347
16: 10359
17: 24894
18: 60523
19: 148284
20: 366319
21: 910726
22: 2278658
23: 5731580
24: 14490245
25: 36797588
26: 93839412
27: 240215803
28: 617105614
29: 1590507121
30: 4111846763
31: 10660307791
32: 27711253769
33: 72214088660
34: 188626236139
35: 493782952902
36: 1295297588128
...
498: 8905549466780876642396073343654258905917693028466999354614122414334803973053490695667380301725459866905052509157458138657523790806693604983304945043446074074823147488033014740909049433466213254508146498309733349188106
499: 24947785035695184159749164581713850084342937257706468262391560672422449715753806019752051214121544475469837920461991186908900888093113539339801359611535867899901287200789925164011521414091897974716359811774035382416687

Seed7

<lang seed7>$ include "seed7_05.s7i"; 

 include "bigint.s7i";

const integer: max_n is 500; const integer: branch is 4;

var array bigInteger: rooted is max_n times 0_; var array bigInteger: unrooted is max_n times 0_;

const proc: tree (in integer: br, in integer: n, in integer: l, in var integer: sum, in bigInteger: cnt) is func 

 local
   var integer: b is 0;
   var integer: m is 0;
   var bigInteger: c is 0_;
   var bigInteger: diff is 0_;
 begin
   for b range br + 1 to branch do
     sum +:= n;
     if sum > max_n or l * 2 >= sum and b >= branch then
       # Prevent unneeded long math.
       b := branch;
     else
       if b = (br + 1) then
         c := rooted[n] * cnt;
       else
         diff := bigInteger conv (b - br);
         c := c * (rooted[n] + pred(diff)) div diff;
       end if;
       if l * 2 < sum then
         unrooted[sum] +:= c;
       end if;
       if b < branch then
         rooted[sum] +:= c;
         for m range n-1 downto 1 do
           tree(b, m, l, sum, c);
         end for;
       end if;
     end if;
   end for;
 end func;

const proc: bicenter (in integer: s) is func 

 begin
   if not odd(s) then
     unrooted[s] +:= (rooted[s div 2] * succ(rooted[s div 2])) >> 1;
   end if;
 end func;

const proc: main is func 

 local
   var bigInteger: cnt is 1_;
   var integer: n is 0;
   var integer: sum is 1;
 begin
   rooted[1] := 1_;
   unrooted[1] := 1_;
   for n range 1 to max_n do
     tree(0, n, n, sum, cnt);
     bicenter(n);
     writeln(n <& ": " <& unrooted[n]);
   end for;
 end func;</lang>

Output (trimmed): 

1: 1
2: 1
3: 1
4: 2
5: 3
6: 5
7: 9
8: 18
9: 35
10: 75
11: 159
12: 355
13: 802
14: 1858
15: 4347
16: 10359
17: 24894
18: 60523
19: 148284
20: 366319
21: 910726
22: 2278658
23: 5731580
24: 14490245
25: 36797588
...
499: 24947785035695184159749164581713850084342937257706468262391560672422449715753806019752051214121544475469837920461991186908900888093113539339801359611535867899901287200789925164011521414091897974716359811774035382416687
500: 69888806977968318652007568872323509883145559195919337637376568556764896550639165374194198834566064897188044072953504164443147167221896769574954040005634802439024034052447587053721967608928388404735140565376686372508743

Tcl

Translation of: C

Handles arbitrarily large values. <lang tcl>package require Tcl 8.5

set maxN 200 set rooted [lrepeat $maxN 0] lset rooted 0 1; lset rooted 1 1 set unrooted $rooted

proc choose {m k} { 

   if {$k == 1} {

return $m 

   }
   for {set r $m; set i 1} {$i < $k} {incr i} {

set r [expr {$r * ($m+$i) / ($i+1)}] 

   }
   return $r

}

proc tree {br n cnt sum l} { 

   global maxN rooted unrooted
   for {set b [expr {$br+1}]} {$b <= 4} {incr b} {

set s [expr {$sum + ($b-$br) * $n}] if {$s >= $maxN} return set c [expr {[choose [lindex $rooted $n] [expr {$b-$br}]] * $cnt}] if {$l*2 < $s} { lset unrooted $s [expr {[lindex $unrooted $s] + $c}] } if {$b == 4} return lset rooted $s [expr {[lindex $rooted $s] + $c}] for {set m $n} {[incr m -1]} {} { tree $b $m $c $s $l } 

   }

}

proc bicenter {s} { 

   if {$s & 1} return
   global unrooted rooted
   set r [lindex $rooted [expr {$s/2}]]
   lset unrooted $s [expr {[lindex $unrooted $s] + $r*($r+1)/2}]

}

for {set n 1} {$n < $maxN} {incr n} { 

   tree 0 $n 1 1 $n
   bicenter $n
   puts "${n}: [lindex $unrooted $n]"

}</lang>

zkl

Translation of: D
Translation of: Go

Uses GMP for big ints, mostly modified in place. Rather slow. <lang zkl>var BN=Import("zklBigNum");

const nMax=100, nBranches=4;

var rooted =(nMax+1).pump(List.createLong(nMax+1).write,BN.fp(0)), 

   unrooted=(nMax+1).pump(List.createLong(nMax+1).write,BN.fp(0));

rooted[0]=BN(1); rooted[1]=BN(1); unrooted[0]=BN(1); unrooted[1]=BN(1);

fcn tree(br,n,l,inSum,cnt){ 

  var c=(nBranches).pump(List().write,0);  // happens only once

  sum := inSum;
  foreach b in ([br + 1 .. nBranches]){
     sum += n;
     if (sum > nMax or (l * 2 >= sum and b >= nBranches)) return();
     if (b == br + 1) c[br] = rooted[n] * cnt; // -->BigInt
     else{

c[br].mul(rooted[n] + b - br - 1); c[br].div(b - br); 

     }
     if (l * 2 < sum) unrooted[sum].add(c[br]);
     if (b < nBranches) rooted[sum].add(c[br]);
     foreach m in ([n-1 .. 1,-1]) { tree(b, m, l, sum, c[br]); }
  }

}

fcn bicenter(s){ 

  if (s.isEven) unrooted[s].add(rooted[s / 2] * (rooted[s / 2] + 1) / 2);

}

foreach n in ([1 .. nMax]){ 

  tree(0, n, n, 1, BN(1));
  bicenter(n);
  println(n, ": ", unrooted[n]);

}</lang>

Output:
1: 1
2: 1
3: 1
4: 2
5: 3
6: 5
7: 9
8: 18
9: 35
10: 75
...
97: 286312976836850192359345859166390622180
98: 785684759853087702778573182234297830503
99: 2156596319845084996862701478402986311496
100: 5921072038125809849884993369103538010139
Retrieved from "https://rosettacode.org/wiki/Paraffins?oldid=208208"