Paraffins: Difference between revisions
Line 976: | Line 976: | ||
for (b = 1, 4-B, |
for (b = 1, 4-B, |
||
if ((s = S + b * n) < P, |
if ((s = S + b * n) < P, |
||
c = R[n+1] * C * prod(i = 1, b-1, (R[n+1]+i) / (i+1)) |
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 (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))))) |
if (B+b < 4, R[s+1] += c; i = n; while (i--, self()(B+b, i, c, s, l))))) |
Revision as of 14:41, 16 December 2016
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:
C
Can't show the tree shapes; count only. <lang c>#include <stdio.h>
- define MAX_N 33 /* max number of tree nodes */
- 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;
- 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>
- include <stdio.h>
- include <stdlib.h>
- define MAX_BRANCH 4
- 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
<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
<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
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;
- state: [unrooted, ra]
- 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
- 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;
- 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
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
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
<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
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
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
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
<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.
<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)
- 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>
- 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
<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
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
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