1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
|
## howmanytrees.R (2022-10-10)
## Calculate Numbers of Phylogenetic Trees
## Copyright 2004-2022 Emmanuel Paradis
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
LargeNumber <- function(a, b)
{
c <- b * log10(a)
n <- floor(c)
x <- 10^(c - n)
structure(c(x = x, n = n), class = "LargeNumber")
}
print.LargeNumber <- function(x, latex = FALSE, digits = 1, ...)
{
if (latex) {
cat("$\\approx ", round(x["x"], digits), " \\times 10^{", x["n"], "}$\n", sep = "")
} else {
cat("approximately ", x["x"], " * 10^", x["n"], "\n", sep = "")
}
}
howmanytrees <- function(n, rooted = TRUE, binary = TRUE,
labeled = TRUE, detail = FALSE)
{
if (!labeled && !(rooted & binary))
stop("can compute number of unlabeled trees only for rooted binary cases.")
if (n < 3) N <- 1 else {
if (labeled) {
if (!rooted) n <- n - 1
if (binary) {
if (n < 152) {
N <- prod(seq(1, 2*n - 3, by = 2)) # double factorial
} else {
ldfac <- lfactorial(2 * n - 3) - (n - 2) * log(2) - lfactorial(n - 2)
N <- LargeNumber(exp(1), ldfac)
}
}
else {
N <- matrix(0, n, n - 1)
N[1:n, 1] <- 1
for (i in 3:n)
for (j in 2:(i - 1))
N[i, j] <- (i + j - 2)*N[i - 1, j - 1] + j*N[i - 1, j]
if (detail) {
rownames(N) <- 1:n
colnames(N) <- 1:(n - 1)
} else N <- sum(N[n, ])
}
} else {
N <- numeric(n)
N[1] <- 1
for (i in 2:n) {
if (i %% 2) {
im1 <- i - 1L
x <- N[1:(im1 / 2)]
y <- N[im1:((i + 1) / 2)]
} else {
ion2 <- i / 2
x <- N[1:ion2]
y <- N[(i - 1):ion2]
ny <- length(y)
y[ny] <- (y[ny] + 1) / 2
}
N[i] <- sum(x * y)
}
if (detail) names(N) <- 1:n else N <- N[n]
}
}
N
}
|
