summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
authorgit <git@karoshi.linehan.me>2018-08-01 14:03:17 +0000
committergit <git@karoshi.linehan.me>2018-08-01 14:03:17 +0000
commitf2dde877f110ebc26458b619fcf06ed02f1a8adc (patch)
tree700048cc61bca39775cc2498c1474ce734ef7571
downloadmqtc-master.tar.gz
mqtc-master.tar.bz2
mqtc-master.zip
Adds filesHEADmaster
-rw-r--r--10x10.txt10
-rw-r--r--18x18.txt18
-rw-r--r--20x20.txt20
-rw-r--r--30x30.txt30
-rw-r--r--Makefile77
-rw-r--r--input.c159
-rw-r--r--input.h6
-rw-r--r--logs.c68
-rw-r--r--logs.h14
-rw-r--r--main.c206
-rw-r--r--prng/alias.c199
-rw-r--r--prng/alias.h23
-rw-r--r--prng/coin.c58
-rw-r--r--prng/coin.h35
-rw-r--r--prng/coin2.c281
-rw-r--r--prng/coin2.h34
-rw-r--r--prng/dice.c40
-rw-r--r--prng/dice.h8
-rw-r--r--prng/mersenne.c168
-rw-r--r--prng/mersenne.h23
-rw-r--r--prng/prng.c46
-rw-r--r--prng/prng.h10
-rw-r--r--prng/tbl/5tbl.c183
-rw-r--r--prng/tbl/5tbl.c.zipbin0 -> 3369 bytes
-rw-r--r--prng/tbl/__MACOSX/._5tbl.cbin0 -> 336 bytes
-rw-r--r--tree/.node_cost.c.swpbin0 -> 20480 bytes
-rw-r--r--tree/node_alloc.c87
-rw-r--r--tree/node_check.c204
-rw-r--r--tree/node_cost.c257
-rw-r--r--tree/node_count.c84
-rw-r--r--tree/node_get.c265
-rw-r--r--tree/node_insert.c272
-rw-r--r--tree/node_mutate.c322
-rw-r--r--tree/node_print.c407
-rw-r--r--tree/node_traverse.c90
-rw-r--r--tree/old/pnode.c31
-rw-r--r--tree/old/pnode.h137
-rw-r--r--tree/old/print.c406
-rw-r--r--tree/old/print.h9
-rw-r--r--tree/old/ptree.c17
-rw-r--r--tree/old/ptree.h73
-rw-r--r--tree/tree_alloc.c157
-rw-r--r--tree/tree_cost.c49
-rw-r--r--tree/tree_mutate.c267
-rw-r--r--tree/ytree.h188
-rw-r--r--util/list.c297
-rw-r--r--util/list.h153
-rw-r--r--util/math.c204
-rw-r--r--util/math.h19
49 files changed, 5711 insertions, 0 deletions
diff --git a/10x10.txt b/10x10.txt
new file mode 100644
index 0000000..f5245fe
--- /dev/null
+++ b/10x10.txt
@@ -0,0 +1,10 @@
+0.000000 0.900000 0.500000 0.800000 0.300000 0.600000 0.800000 0.900000 0.800000 0.400000
+0.900000 0.000000 0.700000 0.600000 0.900000 0.600000 0.400000 0.300000 0.600000 0.800000
+0.500000 0.700000 0.000000 0.600000 0.500000 0.400000 0.600000 0.700000 0.600000 0.400000
+0.800000 0.600000 0.600000 0.000000 0.800000 0.500000 0.500000 0.600000 0.300000 0.700000
+0.300000 0.900000 0.500000 0.800000 0.000000 0.600000 0.800000 0.900000 0.800000 0.400000
+0.600000 0.600000 0.400000 0.500000 0.600000 0.000000 0.500000 0.600000 0.500000 0.500000
+0.800000 0.400000 0.600000 0.500000 0.800000 0.500000 0.000000 0.400000 0.500000 0.700000
+0.900000 0.300000 0.700000 0.600000 0.900000 0.600000 0.400000 0.000000 0.600000 0.800000
+0.800000 0.600000 0.600000 0.300000 0.800000 0.500000 0.500000 0.600000 0.000000 0.700000
+0.400000 0.800000 0.400000 0.700000 0.400000 0.500000 0.700000 0.800000 0.700000 0.000000
diff --git a/18x18.txt b/18x18.txt
new file mode 100644
index 0000000..a6d6e72
--- /dev/null
+++ b/18x18.txt
@@ -0,0 +1,18 @@
+0.000000 0.277778 0.444444 0.444444 0.388889 0.444444 0.388889 0.388889 0.333333 0.444444 0.166667 0.333333 0.500000 0.555556 0.555556 0.388889 0.388889 0.444444
+0.277778 0.000000 0.333333 0.333333 0.388889 0.333333 0.277778 0.388889 0.333333 0.333333 0.277778 0.333333 0.388889 0.444444 0.444444 0.388889 0.388889 0.333333
+0.444444 0.333333 0.000000 0.388889 0.555556 0.388889 0.222222 0.555556 0.500000 0.166667 0.444444 0.500000 0.444444 0.500000 0.500000 0.555556 0.555556 0.388889
+0.444444 0.333333 0.388889 0.000000 0.555556 0.166667 0.333333 0.555556 0.500000 0.388889 0.444444 0.500000 0.333333 0.388889 0.388889 0.555556 0.555556 0.277778
+0.388889 0.388889 0.555556 0.555556 0.000000 0.555556 0.500000 0.388889 0.333333 0.555556 0.388889 0.222222 0.611111 0.666667 0.666667 0.388889 0.166667 0.555556
+0.444444 0.333333 0.388889 0.166667 0.555556 0.000000 0.333333 0.555556 0.500000 0.388889 0.444444 0.500000 0.333333 0.388889 0.388889 0.555556 0.555556 0.277778
+0.388889 0.277778 0.222222 0.333333 0.500000 0.333333 0.000000 0.500000 0.444444 0.222222 0.388889 0.444444 0.388889 0.444444 0.444444 0.500000 0.500000 0.333333
+0.388889 0.388889 0.555556 0.555556 0.388889 0.555556 0.500000 0.000000 0.222222 0.555556 0.388889 0.333333 0.611111 0.666667 0.666667 0.166667 0.388889 0.555556
+0.333333 0.333333 0.500000 0.500000 0.333333 0.500000 0.444444 0.222222 0.000000 0.500000 0.333333 0.277778 0.555556 0.611111 0.611111 0.222222 0.333333 0.500000
+0.444444 0.333333 0.166667 0.388889 0.555556 0.388889 0.222222 0.555556 0.500000 0.000000 0.444444 0.500000 0.444444 0.500000 0.500000 0.555556 0.555556 0.388889
+0.166667 0.277778 0.444444 0.444444 0.388889 0.444444 0.388889 0.388889 0.333333 0.444444 0.000000 0.333333 0.500000 0.555556 0.555556 0.388889 0.388889 0.444444
+0.333333 0.333333 0.500000 0.500000 0.222222 0.500000 0.444444 0.333333 0.277778 0.500000 0.333333 0.000000 0.555556 0.611111 0.611111 0.333333 0.222222 0.500000
+0.500000 0.388889 0.444444 0.333333 0.611111 0.333333 0.388889 0.611111 0.555556 0.444444 0.500000 0.555556 0.000000 0.222222 0.222222 0.611111 0.611111 0.222222
+0.555556 0.444444 0.500000 0.388889 0.666667 0.388889 0.444444 0.666667 0.611111 0.500000 0.555556 0.611111 0.222222 0.000000 0.166667 0.666667 0.666667 0.277778
+0.555556 0.444444 0.500000 0.388889 0.666667 0.388889 0.444444 0.666667 0.611111 0.500000 0.555556 0.611111 0.222222 0.166667 0.000000 0.666667 0.666667 0.277778
+0.388889 0.388889 0.555556 0.555556 0.388889 0.555556 0.500000 0.166667 0.222222 0.555556 0.388889 0.333333 0.611111 0.666667 0.666667 0.000000 0.388889 0.555556
+0.388889 0.388889 0.555556 0.555556 0.166667 0.555556 0.500000 0.388889 0.333333 0.555556 0.388889 0.222222 0.611111 0.666667 0.666667 0.388889 0.000000 0.555556
+0.444444 0.333333 0.388889 0.277778 0.555556 0.277778 0.333333 0.555556 0.500000 0.388889 0.444444 0.500000 0.222222 0.277778 0.277778 0.555556 0.555556 0.000000
diff --git a/20x20.txt b/20x20.txt
new file mode 100644
index 0000000..3b78fac
--- /dev/null
+++ b/20x20.txt
@@ -0,0 +1,20 @@
+0.000000 0.350000 0.450000 0.400000 0.400000 0.200000 0.400000 0.400000 0.350000 0.350000 0.200000 0.450000 0.350000 0.450000 0.350000 0.350000 0.300000 0.250000 0.450000 0.350000
+0.350000 0.000000 0.650000 0.200000 0.600000 0.300000 0.600000 0.600000 0.550000 0.250000 0.400000 0.250000 0.250000 0.650000 0.350000 0.350000 0.300000 0.450000 0.250000 0.550000
+0.450000 0.650000 0.000000 0.700000 0.200000 0.500000 0.400000 0.400000 0.250000 0.650000 0.400000 0.750000 0.650000 0.150000 0.650000 0.650000 0.600000 0.350000 0.750000 0.350000
+0.400000 0.200000 0.700000 0.000000 0.650000 0.350000 0.650000 0.650000 0.600000 0.300000 0.450000 0.200000 0.300000 0.700000 0.400000 0.400000 0.350000 0.500000 0.200000 0.600000
+0.400000 0.600000 0.200000 0.650000 0.000000 0.450000 0.350000 0.350000 0.200000 0.600000 0.350000 0.700000 0.600000 0.200000 0.600000 0.600000 0.550000 0.300000 0.700000 0.300000
+0.200000 0.300000 0.500000 0.350000 0.450000 0.000000 0.450000 0.450000 0.400000 0.300000 0.250000 0.400000 0.300000 0.500000 0.300000 0.300000 0.250000 0.300000 0.400000 0.400000
+0.400000 0.600000 0.400000 0.650000 0.350000 0.450000 0.000000 0.150000 0.300000 0.600000 0.350000 0.700000 0.600000 0.400000 0.600000 0.600000 0.550000 0.300000 0.700000 0.200000
+0.400000 0.600000 0.400000 0.650000 0.350000 0.450000 0.150000 0.000000 0.300000 0.600000 0.350000 0.700000 0.600000 0.400000 0.600000 0.600000 0.550000 0.300000 0.700000 0.200000
+0.350000 0.550000 0.250000 0.600000 0.200000 0.400000 0.300000 0.300000 0.000000 0.550000 0.300000 0.650000 0.550000 0.250000 0.550000 0.550000 0.500000 0.250000 0.650000 0.250000
+0.350000 0.250000 0.650000 0.300000 0.600000 0.300000 0.600000 0.600000 0.550000 0.000000 0.400000 0.350000 0.150000 0.650000 0.350000 0.350000 0.300000 0.450000 0.350000 0.550000
+0.200000 0.400000 0.400000 0.450000 0.350000 0.250000 0.350000 0.350000 0.300000 0.400000 0.000000 0.500000 0.400000 0.400000 0.400000 0.400000 0.350000 0.200000 0.500000 0.300000
+0.450000 0.250000 0.750000 0.200000 0.700000 0.400000 0.700000 0.700000 0.650000 0.350000 0.500000 0.000000 0.350000 0.750000 0.450000 0.450000 0.400000 0.550000 0.150000 0.650000
+0.350000 0.250000 0.650000 0.300000 0.600000 0.300000 0.600000 0.600000 0.550000 0.150000 0.400000 0.350000 0.000000 0.650000 0.350000 0.350000 0.300000 0.450000 0.350000 0.550000
+0.450000 0.650000 0.150000 0.700000 0.200000 0.500000 0.400000 0.400000 0.250000 0.650000 0.400000 0.750000 0.650000 0.000000 0.650000 0.650000 0.600000 0.350000 0.750000 0.350000
+0.350000 0.350000 0.650000 0.400000 0.600000 0.300000 0.600000 0.600000 0.550000 0.350000 0.400000 0.450000 0.350000 0.650000 0.000000 0.150000 0.200000 0.450000 0.450000 0.550000
+0.350000 0.350000 0.650000 0.400000 0.600000 0.300000 0.600000 0.600000 0.550000 0.350000 0.400000 0.450000 0.350000 0.650000 0.150000 0.000000 0.200000 0.450000 0.450000 0.550000
+0.300000 0.300000 0.600000 0.350000 0.550000 0.250000 0.550000 0.550000 0.500000 0.300000 0.350000 0.400000 0.300000 0.600000 0.200000 0.200000 0.000000 0.400000 0.400000 0.500000
+0.250000 0.450000 0.350000 0.500000 0.300000 0.300000 0.300000 0.300000 0.250000 0.450000 0.200000 0.550000 0.450000 0.350000 0.450000 0.450000 0.400000 0.000000 0.550000 0.250000
+0.450000 0.250000 0.750000 0.200000 0.700000 0.400000 0.700000 0.700000 0.650000 0.350000 0.500000 0.150000 0.350000 0.750000 0.450000 0.450000 0.400000 0.550000 0.000000 0.650000
+0.350000 0.550000 0.350000 0.600000 0.300000 0.400000 0.200000 0.200000 0.250000 0.550000 0.300000 0.650000 0.550000 0.350000 0.550000 0.550000 0.500000 0.250000 0.650000 0.000000
diff --git a/30x30.txt b/30x30.txt
new file mode 100644
index 0000000..d0d7525
--- /dev/null
+++ b/30x30.txt
@@ -0,0 +1,30 @@
+0.000000 0.433333 0.400000 0.166667 0.533333 0.166667 0.266667 0.300000 0.233333 0.166667 0.266667 0.200000 0.233333 0.266667 0.266667 0.533333 0.300000 0.166667 0.400000 0.400000 0.266667 0.500000 0.266667 0.233333 0.166667 0.466667 0.300000 0.233333 0.200000 0.333333
+0.433333 0.000000 0.133333 0.433333 0.200000 0.433333 0.400000 0.233333 0.433333 0.500000 0.266667 0.533333 0.433333 0.466667 0.333333 0.200000 0.500000 0.500000 0.200000 0.200000 0.400000 0.166667 0.333333 0.366667 0.500000 0.133333 0.500000 0.433333 0.533333 0.200000
+0.400000 0.133333 0.000000 0.400000 0.233333 0.400000 0.366667 0.200000 0.400000 0.466667 0.233333 0.500000 0.400000 0.433333 0.300000 0.233333 0.466667 0.466667 0.166667 0.166667 0.366667 0.200000 0.300000 0.333333 0.466667 0.166667 0.466667 0.400000 0.500000 0.166667
+0.166667 0.433333 0.400000 0.000000 0.533333 0.100000 0.266667 0.300000 0.233333 0.233333 0.266667 0.266667 0.233333 0.266667 0.266667 0.533333 0.300000 0.233333 0.400000 0.400000 0.266667 0.500000 0.266667 0.233333 0.233333 0.466667 0.300000 0.233333 0.266667 0.333333
+0.533333 0.200000 0.233333 0.533333 0.000000 0.533333 0.500000 0.333333 0.533333 0.600000 0.366667 0.633333 0.533333 0.566667 0.433333 0.100000 0.600000 0.600000 0.300000 0.300000 0.500000 0.133333 0.433333 0.466667 0.600000 0.166667 0.600000 0.533333 0.633333 0.300000
+0.166667 0.433333 0.400000 0.100000 0.533333 0.000000 0.266667 0.300000 0.233333 0.233333 0.266667 0.266667 0.233333 0.266667 0.266667 0.533333 0.300000 0.233333 0.400000 0.400000 0.266667 0.500000 0.266667 0.233333 0.233333 0.466667 0.300000 0.233333 0.266667 0.333333
+0.266667 0.400000 0.366667 0.266667 0.500000 0.266667 0.000000 0.266667 0.266667 0.333333 0.233333 0.366667 0.266667 0.300000 0.233333 0.500000 0.333333 0.333333 0.366667 0.366667 0.100000 0.466667 0.233333 0.133333 0.333333 0.433333 0.333333 0.266667 0.366667 0.300000
+0.300000 0.233333 0.200000 0.300000 0.333333 0.300000 0.266667 0.000000 0.300000 0.366667 0.133333 0.400000 0.300000 0.333333 0.200000 0.333333 0.366667 0.366667 0.200000 0.200000 0.266667 0.300000 0.200000 0.233333 0.366667 0.266667 0.366667 0.300000 0.400000 0.133333
+0.233333 0.433333 0.400000 0.233333 0.533333 0.233333 0.266667 0.300000 0.000000 0.300000 0.266667 0.333333 0.100000 0.200000 0.266667 0.533333 0.233333 0.300000 0.400000 0.400000 0.266667 0.500000 0.266667 0.233333 0.300000 0.466667 0.233333 0.166667 0.333333 0.333333
+0.166667 0.500000 0.466667 0.233333 0.600000 0.233333 0.333333 0.366667 0.300000 0.000000 0.333333 0.200000 0.300000 0.333333 0.333333 0.600000 0.366667 0.166667 0.466667 0.466667 0.333333 0.566667 0.333333 0.300000 0.100000 0.533333 0.366667 0.300000 0.200000 0.400000
+0.266667 0.266667 0.233333 0.266667 0.366667 0.266667 0.233333 0.133333 0.266667 0.333333 0.000000 0.366667 0.266667 0.300000 0.166667 0.366667 0.333333 0.333333 0.233333 0.233333 0.233333 0.333333 0.166667 0.200000 0.333333 0.300000 0.333333 0.266667 0.366667 0.166667
+0.200000 0.533333 0.500000 0.266667 0.633333 0.266667 0.366667 0.400000 0.333333 0.200000 0.366667 0.000000 0.333333 0.366667 0.366667 0.633333 0.400000 0.133333 0.500000 0.500000 0.366667 0.600000 0.366667 0.333333 0.200000 0.566667 0.400000 0.333333 0.100000 0.433333
+0.233333 0.433333 0.400000 0.233333 0.533333 0.233333 0.266667 0.300000 0.100000 0.300000 0.266667 0.333333 0.000000 0.200000 0.266667 0.533333 0.233333 0.300000 0.400000 0.400000 0.266667 0.500000 0.266667 0.233333 0.300000 0.466667 0.233333 0.166667 0.333333 0.333333
+0.266667 0.466667 0.433333 0.266667 0.566667 0.266667 0.300000 0.333333 0.200000 0.333333 0.300000 0.366667 0.200000 0.000000 0.300000 0.566667 0.133333 0.333333 0.433333 0.433333 0.300000 0.533333 0.300000 0.266667 0.333333 0.500000 0.133333 0.133333 0.366667 0.366667
+0.266667 0.333333 0.300000 0.266667 0.433333 0.266667 0.233333 0.200000 0.266667 0.333333 0.166667 0.366667 0.266667 0.300000 0.000000 0.433333 0.333333 0.333333 0.300000 0.300000 0.233333 0.400000 0.100000 0.200000 0.333333 0.366667 0.333333 0.266667 0.366667 0.233333
+0.533333 0.200000 0.233333 0.533333 0.100000 0.533333 0.500000 0.333333 0.533333 0.600000 0.366667 0.633333 0.533333 0.566667 0.433333 0.000000 0.600000 0.600000 0.300000 0.300000 0.500000 0.133333 0.433333 0.466667 0.600000 0.166667 0.600000 0.533333 0.633333 0.300000
+0.300000 0.500000 0.466667 0.300000 0.600000 0.300000 0.333333 0.366667 0.233333 0.366667 0.333333 0.400000 0.233333 0.133333 0.333333 0.600000 0.000000 0.366667 0.466667 0.466667 0.333333 0.566667 0.333333 0.300000 0.366667 0.533333 0.100000 0.166667 0.400000 0.400000
+0.166667 0.500000 0.466667 0.233333 0.600000 0.233333 0.333333 0.366667 0.300000 0.166667 0.333333 0.133333 0.300000 0.333333 0.333333 0.600000 0.366667 0.000000 0.466667 0.466667 0.333333 0.566667 0.333333 0.300000 0.166667 0.533333 0.366667 0.300000 0.133333 0.400000
+0.400000 0.200000 0.166667 0.400000 0.300000 0.400000 0.366667 0.200000 0.400000 0.466667 0.233333 0.500000 0.400000 0.433333 0.300000 0.300000 0.466667 0.466667 0.000000 0.100000 0.366667 0.266667 0.300000 0.333333 0.466667 0.233333 0.466667 0.400000 0.500000 0.166667
+0.400000 0.200000 0.166667 0.400000 0.300000 0.400000 0.366667 0.200000 0.400000 0.466667 0.233333 0.500000 0.400000 0.433333 0.300000 0.300000 0.466667 0.466667 0.100000 0.000000 0.366667 0.266667 0.300000 0.333333 0.466667 0.233333 0.466667 0.400000 0.500000 0.166667
+0.266667 0.400000 0.366667 0.266667 0.500000 0.266667 0.100000 0.266667 0.266667 0.333333 0.233333 0.366667 0.266667 0.300000 0.233333 0.500000 0.333333 0.333333 0.366667 0.366667 0.000000 0.466667 0.233333 0.133333 0.333333 0.433333 0.333333 0.266667 0.366667 0.300000
+0.500000 0.166667 0.200000 0.500000 0.133333 0.500000 0.466667 0.300000 0.500000 0.566667 0.333333 0.600000 0.500000 0.533333 0.400000 0.133333 0.566667 0.566667 0.266667 0.266667 0.466667 0.000000 0.400000 0.433333 0.566667 0.133333 0.566667 0.500000 0.600000 0.266667
+0.266667 0.333333 0.300000 0.266667 0.433333 0.266667 0.233333 0.200000 0.266667 0.333333 0.166667 0.366667 0.266667 0.300000 0.100000 0.433333 0.333333 0.333333 0.300000 0.300000 0.233333 0.400000 0.000000 0.200000 0.333333 0.366667 0.333333 0.266667 0.366667 0.233333
+0.233333 0.366667 0.333333 0.233333 0.466667 0.233333 0.133333 0.233333 0.233333 0.300000 0.200000 0.333333 0.233333 0.266667 0.200000 0.466667 0.300000 0.300000 0.333333 0.333333 0.133333 0.433333 0.200000 0.000000 0.300000 0.400000 0.300000 0.233333 0.333333 0.266667
+0.166667 0.500000 0.466667 0.233333 0.600000 0.233333 0.333333 0.366667 0.300000 0.100000 0.333333 0.200000 0.300000 0.333333 0.333333 0.600000 0.366667 0.166667 0.466667 0.466667 0.333333 0.566667 0.333333 0.300000 0.000000 0.533333 0.366667 0.300000 0.200000 0.400000
+0.466667 0.133333 0.166667 0.466667 0.166667 0.466667 0.433333 0.266667 0.466667 0.533333 0.300000 0.566667 0.466667 0.500000 0.366667 0.166667 0.533333 0.533333 0.233333 0.233333 0.433333 0.133333 0.366667 0.400000 0.533333 0.000000 0.533333 0.466667 0.566667 0.233333
+0.300000 0.500000 0.466667 0.300000 0.600000 0.300000 0.333333 0.366667 0.233333 0.366667 0.333333 0.400000 0.233333 0.133333 0.333333 0.600000 0.100000 0.366667 0.466667 0.466667 0.333333 0.566667 0.333333 0.300000 0.366667 0.533333 0.000000 0.166667 0.400000 0.400000
+0.233333 0.433333 0.400000 0.233333 0.533333 0.233333 0.266667 0.300000 0.166667 0.300000 0.266667 0.333333 0.166667 0.133333 0.266667 0.533333 0.166667 0.300000 0.400000 0.400000 0.266667 0.500000 0.266667 0.233333 0.300000 0.466667 0.166667 0.000000 0.333333 0.333333
+0.200000 0.533333 0.500000 0.266667 0.633333 0.266667 0.366667 0.400000 0.333333 0.200000 0.366667 0.100000 0.333333 0.366667 0.366667 0.633333 0.400000 0.133333 0.500000 0.500000 0.366667 0.600000 0.366667 0.333333 0.200000 0.566667 0.400000 0.333333 0.000000 0.433333
+0.333333 0.200000 0.166667 0.333333 0.300000 0.333333 0.300000 0.133333 0.333333 0.400000 0.166667 0.433333 0.333333 0.366667 0.233333 0.300000 0.400000 0.400000 0.166667 0.166667 0.300000 0.266667 0.233333 0.266667 0.400000 0.233333 0.400000 0.333333 0.433333 0.000000
diff --git a/Makefile b/Makefile
new file mode 100644
index 0000000..49dc440
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,77 @@
+#########################
+# Configure project
+#########################
+PROGRAM=mqtc
+INCLUDE=
+
+
+#########################
+# Configure build
+#########################
+
+COMPILER=gcc
+#
+# optimize enable all include
+# level 3 warnings paths see note
+# \ | | /
+CC_FLAGS=-O3 -Wall $(INCLUDE) #-ffast-math
+LD_FLAGS=-lm
+# /
+# math
+#
+#
+# NOTE on -ffast-math
+#
+# First, breaks strict IEEE compliance, e.g. allows re-ordering of
+# instructions to a mathematical equivalent, which may not be IEEE
+# floating-point equivalent.
+#
+# Second, disables setting errno after single-instruction math functions,
+# avoiding a write to a thread-local variable (can produce 100% speedup on
+# certain architectures).
+#
+# Third, assumes finite math only, meaning no checks for NaN (or 0) are
+# made where they would normally be. It is assumed these values will never
+# appear.
+#
+# Fourth, enables reciprocal approximations for division and reciprocal
+# square root.
+#
+# Fifth, disables signed zero (even if the compile target supports it)
+# and rounding math, which enables optimizations e.g. constant folding.
+#
+# Sixth, generates code assuming no hardware interrupts occur in math
+# due to signal()/trap(). If these cannot be disabled on the compile
+# target and consequently occur, they will not be handled.
+#
+
+#########################
+# Configure files
+#########################
+
+SOURCES=main.c \
+ input.c \
+ logs.c \
+ util/list.c \
+ util/math.c \
+ prng/mersenne.c \
+ prng/prng.c \
+ prng/coin.c \
+ prng/dice.c \
+ prng/alias.c \
+ tree/pnode.c \
+ tree/ptree.c \
+ tree/print.c
+
+STATICS=
+OBJECTS=$(SOURCES:.c=.o)
+
+
+#########################
+# Configure rules
+#########################
+all: $(SOURCES)
+ $(COMPILER) $(CC_FLAGS) $(SOURCES) $(STATICS) -o $(PROGRAM) $(LD_FLAGS)
+
+clean:
+ rm -f $(OBJECTS) $(PROGRAM) gmon.out
diff --git a/input.c b/input.c
new file mode 100644
index 0000000..250676d
--- /dev/null
+++ b/input.c
@@ -0,0 +1,159 @@
+#include <stdlib.h>
+#include <stdio.h>
+#include <stdbool.h>
+#include <string.h>
+
+int read_float_line(FILE *input, float **vector, int *count)
+{
+ float data_buffer[4096];
+ char line_buffer[4096];
+
+ char *ptr;
+ float datum;
+ int status;
+ int i = 0;
+ int n = 0;
+
+ if (NULL != fgets(line_buffer, 4096, input)) {
+
+ ptr = line_buffer;
+
+ for (;;) {
+
+ /* Skip any leading whitespace */
+ while ((ptr[i]==' ' || ptr[i]=='\t')) {
+ i++;
+ }
+
+ /*
+ * If we hit a newline, that's
+ * the end of the loop.
+ */
+ if (ptr[i] == '\n') {
+ break;
+ } else {
+ /*
+ * Otherwise, we must have
+ * stopped in front of a
+ * non-whitespace character.
+ *
+ * Let's scan it.
+ */
+ status = sscanf(&(ptr[i]), "%f", &datum);
+
+ if (status == 0 || status == EOF) {
+ printf("Scan error.\n");
+ return 0;
+ } else {
+ /* Add the new item to our data */
+ data_buffer[n] = datum;
+ n++;
+
+ /*
+ * Scan until the end of the
+ * thing we just read.
+ */
+ while ((ptr[i]!=' ' && ptr[i]!='\n')) {
+ i++;
+ }
+ }
+ }
+ }
+ } else {
+ return EOF;
+ }
+
+ /*
+ * If the user gave us a count pointer
+ * to fill in, go ahead and do that.
+ */
+ if (count != NULL) {
+ *count = n;
+ }
+
+ /*
+ * If the user gave us a NULL pointer,
+ * allocate some memory for it and
+ * put the contents of data_buffer at
+ * that memory location.
+ */
+ if (vector != NULL) {
+ if (*vector == NULL) {
+ *vector = calloc(1, n*sizeof(float));
+ }
+ for (i=0; i<n; i++) {
+ (*vector)[i] = data_buffer[i];
+ }
+ }
+
+ return 1;
+}
+
+
+
+/**
+ * read_square_matrix()
+ * --------------------
+ * Read a (possibly un-seekable) input stream into memory.
+ *
+ * @input: Input file stream
+ * @count: Number of rows in the matrix
+ * Return: Square matrix of floating-point values.
+ *
+ * NOTE
+ * The method used for reading and allocating is slightly
+ * complicated, because you can't just read the stream
+ * once to count the number of '\n' characters, and use
+ * that to determine the number of lines, then allocate
+ * them, then read it again. STDIN, for example, cannot be
+ * seeked on, so functions like rewind() will not work.
+ */
+float **read_square_matrix(FILE *input, int *count)
+{
+ float **matrix;
+ float *first_line = NULL;
+ int n = 0;
+ int i = 0;
+
+ /*
+ * Parse the first line into an array of floats.
+ */
+ if (EOF == read_float_line(input, &first_line, &n)) {
+ fprintf(stderr, "Couldn't read first line.\n");
+ return NULL;
+ }
+
+ /*
+ * Now we know how many fields there are, we can
+ * allocate the right number of row pointers
+ * (since the matrix is square).
+ */
+ if (NULL == (matrix=calloc(n, sizeof(float *)))) {
+ fprintf(stderr, "Out of memory.\n");
+ return NULL;
+ }
+
+ /*
+ * Assign the first row (We already scanned it in)
+ */
+ matrix[0] = first_line;
+
+ /*
+ * Scan the remaining lines into the allocated
+ * matrix region.
+ */
+ for (i=1; i<n; i++) {
+ float *line = NULL;
+ if (EOF != read_float_line(input, &line, NULL)) {
+ matrix[i] = line;
+ }
+ }
+
+ /* Assign for the caller to have a count. */
+ if (count != NULL) {
+ *count = n;
+ }
+
+ return matrix;
+}
+
diff --git a/input.h b/input.h
new file mode 100644
index 0000000..5553fba
--- /dev/null
+++ b/input.h
@@ -0,0 +1,6 @@
+#ifndef __MQTC_INPUT
+#define __MQTC_INPUT
+
+float **read_square_matrix(FILE *input, int *count);
+
+#endif
diff --git a/logs.c b/logs.c
new file mode 100644
index 0000000..9cc16ab
--- /dev/null
+++ b/logs.c
@@ -0,0 +1,68 @@
+#include <stdlib.h>
+#include <stdio.h>
+#include <stdarg.h>
+
+FILE *F_cost;
+FILE *F_fitness;
+FILE *F_alias;
+FILE *F_mutate;
+
+void open_logs(void)
+{
+ F_alias = fopen("./log/alias.log", "w+");
+ F_mutate = fopen("./log/mutate.log", "w+");
+ F_fitness = fopen("./log/fitness.log", "w+");
+ F_cost = fopen("./log/cost.log", "w+");
+}
+
+void close_logs(void)
+{
+ fclose(F_cost);
+ fclose(F_alias);
+ fclose(F_mutate);
+ fclose(F_fitness);
+}
+
+void log_cost(const char *fmt, ...)
+{
+ va_list ap;
+
+ if (F_cost != NULL && ferror(F_cost) == 0) {
+ va_start(ap, fmt);
+ vfprintf(F_cost, fmt, ap);
+ va_end(ap);
+ }
+}
+
+void log_fitness(const char *fmt, ...)
+{
+ va_list ap;
+
+ if (F_fitness != NULL && ferror(F_fitness) == 0) {
+ va_start(ap, fmt);
+ vfprintf(F_fitness, fmt, ap);
+ va_end(ap);
+ }
+}
+
+void log_alias(const char *fmt, ...)
+{
+ va_list ap;
+
+ if (F_alias != NULL && ferror(F_alias) == 0) {
+ va_start(ap, fmt);
+ vfprintf(F_alias, fmt, ap);
+ va_end(ap);
+ }
+}
+
+void log_mutate(const char *fmt, ...)
+{
+ va_list ap;
+
+ if (F_mutate != NULL && ferror(F_mutate) == 0) {
+ va_start(ap, fmt);
+ vfprintf(F_mutate, fmt, ap);
+ va_end(ap);
+ }
+}
diff --git a/logs.h b/logs.h
new file mode 100644
index 0000000..ac5f2c4
--- /dev/null
+++ b/logs.h
@@ -0,0 +1,14 @@
+#ifndef __MQTC_LOGS
+#define __MQTC_LOGS
+#include <stdarg.h>
+
+void open_logs(void);
+void close_logs(void);
+
+void log_cost (const char *fmt, ...);
+void log_fitness(const char *fmt, ...);
+void log_alias (const char *fmt, ...);
+void log_mutate (const char *fmt, ...);
+
+
+#endif
diff --git a/main.c b/main.c
new file mode 100644
index 0000000..9b2ccbe
--- /dev/null
+++ b/main.c
@@ -0,0 +1,206 @@
+#include <stdlib.h>
+#include <stdio.h>
+#include <stdbool.h>
+#include <string.h>
+#include <time.h>
+#include <math.h>
+#include "tree/ytree.h"
+#include "input.h"
+
+int DATA_COUNT;
+
+/**
+ * sufficient_k()
+ * ``````````````
+ * @n : Number of leaf nodes in the phylogenetic tree
+ * Return: Minimum k allowing complete mutation.
+ *
+ * NOTE
+ * For the class {T} of phylogenetic trees with n leaf nodes,
+ * this function returns the minimal k such that for an
+ * arbitrary tree t in {T}, there exists a k-mutation which
+ * can transform t into arbitrary t' also in {T}.
+ *
+ * That is, given at least k mutations, any tree in the class
+ * can be transformed into any other tree also in the class.
+ */
+static inline int sufficient_k(int n)
+{
+ return (5 * n) - 16;
+}
+
+
+/**
+ * build_pmf()
+ * -----------
+ * Construct the probability mass function for the alias.
+ *
+ * @limit: Maximum value to sample in the distribution
+ * Return: Array of probability values, suitable for building an alias.
+ *
+ * NOTE
+ * The actual distribution being used is from [1], namely,
+ *
+ * p(k) = 1 / ( (k+2) * (log(k+2)^2) ).
+ *
+ * This is a shifted version of the more general equation
+ *
+ * p(k) = c / ( (k) * (log(k)^2) ),
+ *
+ * with c ~= 2.1.
+ *
+ * Among the family of smooth analytic functions that can be expressed
+ * as a series of fractional powers and logarithms, this equation lies
+ * on the exact edge of the convergence zone.
+ *
+ * SUM(1/k)
+ * SUM(1/(k*log(k)) ----\ divergent
+ * SUM(1/(k*log(k)*log(log(k))) ----/ functions
+ *
+ * SUM(1/k^2)
+ * SUM(1/(k*(log(k)^2))) ----\ convergent
+ * SUM(1/(k*(log(k)*(log(k)^2)))) ----/ functions
+ *
+ * The probabilities of a discrete probability distribution must sum to 1,
+ * so any function used to describe the distribution must converge as k goes
+ * to infinity.
+ *
+ * 1/(k*(log(k)^2)) is thus one of the maximal "fat tail" distributions that
+ * exists. This will concentrate probability on larger values of k, which is
+ * what we want for constructing k-mutations.
+ */
+float *build_pmf(int limit)
+{
+ float *value; /* value array */
+ float k; /* index */
+ float p; /* sum of probabilities */
+
+ value = calloc(limit, sizeof(float));
+
+ p = 0.0;
+
+ for (k=1; k<limit; k+=1.0) {
+ value[(int)k] = 1.0/((k+2.0)*powf(log2f(k+2.0), 2.0));
+ p += value[(int)k];
+ }
+
+ /* Ensure the probabilities sum to 1. */
+ value[0] = 1.0 - p;
+
+ return value;
+}
+
+/**
+ * run_mutations()
+ * ---------------
+ */
+void run_mutations(int gens, FILE *input)
+{
+ #define N_TREES 3
+ struct ytree_t *tree[N_TREES];
+ struct ytree_t *best_tree;
+ struct ytree_t *champion;
+ struct alias_t *alias;
+ float *prob;
+ float **data;
+ float best_cost = 0.0;
+ float init_cost[N_TREES];
+ float this_cost[N_TREES];
+ int i;
+ int j;
+ int m;
+
+ /*
+ * Read the input matrix and build the
+ * phylogenetic tree from it.
+ */
+
+ data = read_square_matrix(input, &DATA_COUNT);
+
+ for (i=0; i<N_TREES; i++) {
+ tree[i] = ytree_create(DATA_COUNT, data);
+ }
+
+ /*
+ * Once we know DATA_COUNT (set in read_square_matrix()),
+ * we can use that as our N to build the alias with a
+ * 'sufficient' number of possibilities to transform any
+ * tree into any other tree.
+ *
+ * That sufficiency is given by f(n) = 5n-16 [Cilibrasi 2011]
+ *
+ * We use the alias to sample from the non-uniform
+ * probability mass function and obtain the k for
+ * which we apply a k-mutation.
+ */
+
+ prob = build_pmf(sufficient_k(DATA_COUNT));
+ alias = alias_create(sufficient_k(DATA_COUNT), prob);
+
+ /*
+ * Initialize the best tree and the best
+ * cost of the tree.
+ */
+
+ for (i=0; i<N_TREES; i++) {
+ init_cost[i] = ytree_cost_scaled(tree[i]);
+ if (init_cost[i] > best_cost || i == 0) {
+ best_cost = init_cost[i];
+ best_tree = tree[i];
+ }
+ }
+
+ champion = ytree_copy(best_tree);
+
+ /*
+ * Hill-climb for the optimal cost by
+ * mutating the existing tree according
+ * to the non-uniform probability given
+ * in the alias.
+ */
+
+ for (i=0; i<gens; i++) {
+ for (j=0; j<N_TREES; j++) {
+ tree[j] = ytree_mutate_mmc2(tree[j], alias, &m);
+
+ this_cost[j] = ytree_cost_scaled(tree[j]);
+
+ if (this_cost[j] > best_cost) {
+ best_cost = this_cost[j];
+ best_tree = tree[j];
+ }
+ }
+
+ if (best_tree != NULL) {
+ ytree_free(champion);
+ champion = ytree_copy(best_tree);
+ best_tree = NULL;
+ }
+
+ if (best_cost == 1.0) {
+ /* Halt */
+ break;
+ }
+ }
+
+ ynode_print(champion->root, "%d");
+ printf("best:%f init:", best_cost);
+ for (i=0; i<N_TREES; i++) {
+ printf("%f ", init_cost[i]);
+ }
+ printf("\n");
+}
+
+
+int main(int argc, char *argv[])
+{
+ if (argc == 2) {
+ run_mutations(atoi(argv[1]), stdin);
+ return 1;
+ } else {
+ printf("Usage: cat <datafile> | %s <# generations>", argv[0]);
+ }
+
+ return 0;
+}
+
diff --git a/prng/alias.c b/prng/alias.c
new file mode 100644
index 0000000..e6acee3
--- /dev/null
+++ b/prng/alias.c
@@ -0,0 +1,199 @@
+#include <stdlib.h>
+#include <math.h>
+#include "../util/list.h"
+#include "dice.h"
+#include "coin.h"
+#include "alias.h"
+
+/******************************************************************************
+ * ALIAS METHOD
+ * ------------
+ * Allows the sampling of a non-uniform discrete probability
+ * distribution.
+ *
+ ******************************************************************************/
+
+struct alias_item {
+ struct list_t node;
+ float pr;
+ int idx;
+};
+
+
+/**
+ * alias_create()
+ * --------------
+ * Create a new alias structure, and run the setup function.
+ *
+ * @prob_count: Number of values in distribution (length of @prob_array)
+ * @prob_array: Discrete probability distribution.
+ * Return : Pointer to alias structure.
+ */
+struct alias_t *alias_create(int prob_count, float *prob_array)
+{
+ struct alias_t *alias;
+
+ alias = calloc(1, sizeof(struct alias_t));
+ alias->n = (float)prob_count;
+ alias->pr = calloc(prob_count, sizeof(float));
+ alias->al = calloc(prob_count, sizeof(float));
+
+ alias_setup(alias, prob_array);
+
+ return alias;
+}
+
+
+/**
+ * alias_destroy()
+ * ---------------
+ * Free the probability and alias tables associated with an alias structure.
+ *
+ * @alias: Pointer to an alias structure.
+ * Return: Nothing
+ */
+void alias_destroy(struct alias_t *alias)
+{
+ if (alias->pr != NULL) {
+ free(alias->pr);
+ }
+ if (alias->al != NULL) {
+ free(alias->al);
+ }
+ if (alias != NULL) {
+ free(alias);
+ }
+ return;
+}
+
+
+/**
+ * alias_setup()
+ * -------------
+ * Setup the probability and alias tables for an alias structure.
+ *
+ * @alias : Pointer to the alias structure.
+ * @prob_array: Discrete probability distribution.
+ * Return : Nothing.
+ */
+void alias_setup(struct alias_t *alias, float *prob_array)
+{
+ struct list_t small; /* Small values worklist */
+ struct list_t large; /* Large values worklist */
+ struct alias_item *sm;
+ struct alias_item *lg;
+ int i;
+
+ list_init(&small);
+ list_init(&large);
+
+ /*
+ * Load items from the array into the
+ * two worklists.
+ */
+
+ for (i=0; i<alias->n; i++) {
+
+ struct alias_item *item;
+
+ item = calloc(1, sizeof(struct alias_item));
+ item->pr = prob_array[i] * (float)alias->n; /* SCALED */
+ item->idx = i;
+
+ if (item->pr<1.0) {
+ list_add(&small, item);
+ } else {
+ list_add(&large, item);
+ }
+ }
+
+ /*
+ * Build the probability and alias tables.
+ *
+ * Depending on the probability distribution,
+ * 'small' or 'large' may empty first, so we
+ * treat them together, while it lasts, and
+ * then separately.
+ */
+
+ for (i=0; !list_empty(&small) && !list_empty(&large); i++) {
+ /*
+ * Small and large not empty
+ */
+ sm = (struct alias_item *)list_head(&small);
+ lg = (struct alias_item *)list_head(&large);
+
+ list_remove_from(&small, sm);
+ list_remove_from(&large, lg);
+
+ alias->pr[sm->idx] = sm->pr;
+ alias->al[sm->idx] = lg->idx;
+
+ lg->pr = (lg->pr + sm->pr) - 1.0;
+
+ if (lg->pr<1.0) {
+ list_add(&small, lg);
+ } else {
+ list_add(&large, lg);
+ }
+
+ free(sm);
+ }
+
+ for (; !list_empty(&large); i++) {
+ /*
+ * Large is not empty
+ */
+ lg = (struct alias_item *)list_head(&large);
+ list_remove_from(&large, lg);
+
+ alias->pr[lg->idx] = 1.0;
+
+ free(lg);
+ }
+
+ for (; !list_empty(&small); i++) {
+ /*
+ * Small is not empty
+ */
+ sm = (struct alias_item *)list_head(&small);
+ list_remove_from(&small, sm);
+
+ alias->pr[sm->idx] = 1.0;
+
+ free(sm);
+ }
+
+ /*
+ * All worklist memory should be free'd
+ */
+ return;
+}
+
+
+/**
+ * alias_sample()
+ * --------------
+ * Sample from a discrete probability distribution.
+ *
+ * @alias: Pointer to alias structure.
+ * Return: Random integer value from (0,1,2,...,@alias->n).
+ */
+int alias_sample(struct alias_t *alias)
+{
+ int value;
+ int i;
+
+ for (i=0; i<alias->n; i++) {
+ value = dice_roll(alias->n);
+
+ if (coin_flip(alias->pr[value], HEADS) == HEADS) {
+ return value;
+ } else {
+ return (int)alias->al[value];
+ }
+ }
+
+ return (alias->n+1); /* impossible */
+}
+
diff --git a/prng/alias.h b/prng/alias.h
new file mode 100644
index 0000000..cc8f727
--- /dev/null
+++ b/prng/alias.h
@@ -0,0 +1,23 @@
+#ifndef __MATH_ALIAS_H
+#define __MATH_ALIAS_H
+
+/******************************************************************************
+ * VOSE-WALKER ALIAS METHOD
+ * ------------------------
+ * Allow the sampling of a non-uniform discrete probability
+ * distribution.
+ *
+ ******************************************************************************/
+
+struct alias_t {
+ int n;
+ float *pr;
+ float *al;
+};
+
+struct alias_t *alias_create (int prob_count, float *prob_array);
+void alias_setup (struct alias_t *alias, float *prob_array);
+void alias_destroy(struct alias_t *alias);
+int alias_sample (struct alias_t *alias);
+
+#endif
diff --git a/prng/coin.c b/prng/coin.c
new file mode 100644
index 0000000..9c90d89
--- /dev/null
+++ b/prng/coin.c
@@ -0,0 +1,58 @@
+#include <stdlib.h>
+#include "prng.h"
+#include "coin.h"
+
+/**
+ * coin_flip()
+ * -----------
+ * Simulate a biased coin flip, applying bias to a given side.
+ *
+ * @bias : Probability of @side appearing
+ * @side : 0 (heads) or 1 (tails).
+ * Return: 0 (heads) or 1 (tails).
+ *
+ * NOTE
+ * To flip a coin with arbitrary-valued sides (rather than boolean),
+ * see the 'coin_flip_value' macro, defined in coin.h.
+ */
+int coin_flip(double bias, int side)
+{
+ if (side == HEADS) {
+ return (prng_uniform_random()<bias) ? HEADS : TAILS;
+ } else if (side == TAILS) {
+ return (prng_uniform_random()<bias) ? TAILS : HEADS;
+ } else {
+ /* Error-ish */
+ return side;
+ }
+}
+
+
+/**
+ * coin_heads()
+ * ------------
+ * Simulate a biased coin flip, applying the bias to heads.
+ *
+ * @bias : Probability of 0 (heads) appearing
+ * Return: 0 (heads) or 1 (tails).
+ */
+int coin_heads(double bias)
+{
+ return (prng_uniform_random()<bias) ? HEADS : TAILS;
+}
+
+
+/**
+ * coin_tails()
+ * ------------
+ * Simulate a biased coin flip, applying the bias to tails.
+ *
+ * @bias : Probability of 1 (tails) appearing
+ * Return: 0 (heads) or 1 (tails).
+ */
+int coin_tails(double bias)
+{
+ return (prng_uniform_random()<bias) ? TAILS : HEADS;
+}
+
+
diff --git a/prng/coin.h b/prng/coin.h
new file mode 100644
index 0000000..1175bf3
--- /dev/null
+++ b/prng/coin.h
@@ -0,0 +1,35 @@
+#ifndef __MATH_COIN_H
+#define __MATH_COIN_H
+
+#define HEADS 0
+#define TAILS 1
+
+int coin_flip (double bias, int side);
+int coin_heads(double bias);
+int coin_tails(double bias);
+
+/**
+ * coin_flip_value()
+ * -----------------
+ * Polymorphic macro function allowing choice between two arbitrary values.
+ *
+ * @bias : Probability of @s1 being chosen.
+ * @s1 : Arbitrary value representing 'HEADS'
+ * @s2 : Arbitrary value representing 'TAILS'
+ * Return: Nothing, evaluates to either @s1 or @s2.
+ */
+#define coin_flip_value(bias, s1, s2) \
+ (coin_flip_heads(bias)) ? s1 : s2
+
+
+/**
+ * coin_fair()
+ * -----------
+ * Flip a fair coin.
+ * Return: 0 (heads) or 1 (tails).
+ */
+#define coin_fair() \
+ coin_heads(0.5)
+
+
+#endif
diff --git a/prng/coin2.c b/prng/coin2.c
new file mode 100644
index 0000000..1696b05
--- /dev/null
+++ b/prng/coin2.c
@@ -0,0 +1,281 @@
+#include <stdlib.h>
+#include <math.h>
+#include "prng.h"
+#include "../util/list.h"
+
+/******************************************************************************
+ * GENERAL M-COIN WITH ARBITRARY DISCRETE BIAS
+ *
+ * An M-coin is an abstract probabilistic object, seen
+ * as a coin with M sides, with each side having some
+ * probability of appearing when the coin is flipped.
+ *
+ * A flip of an M-coin is just a way of talking about
+ * taking a random sample from a discrete probability
+ * distribution.
+ *
+ * Here, the sampling and construction of the probability
+ * distribution is done using the Vose-Walker alias method.
+ *
+ ******************************************************************************/
+
+struct alias_item {
+ struct list_t node;
+ float pr;
+ int idx;
+};
+
+/**
+ * coin_setup()
+ * ------------
+ * Setup the probability and alias tables for a coin structure.
+ *
+ * @coin : Pointer to the coin structure.
+ * @prob_array: Discrete probability distribution.
+ * Return : Nothing.
+ */
+void coin_setup(struct coin_t *coin, int num_sides, float *prob_array)
+{
+ struct list_t small; /* Small values worklist */
+ struct list_t large; /* Large values worklist */
+ struct alias_item *sm;
+ struct alias_item *lg;
+ int i;
+
+ list_init(&small);
+ list_init(&large);
+
+ /*
+ * Allocate resources for the coin struct.
+ */
+
+ coin->n = num_sides;
+
+ /* Don't allocate twice */
+ if (coin->pr != NULL) {
+ free(coin->pr);
+ }
+ if (coin->al != NULL) {
+ free(coin->al);
+ }
+
+ coin->pr = calloc(num_sides, sizeof(float));
+ coin->al = calloc(num_sides, sizeof(float));
+
+ /*
+ * Load items from the array into the
+ * two worklists.
+ */
+
+ for (i=0; i<coin->n; i++) {
+
+ struct alias_item *item;
+
+ item = calloc(1, sizeof(struct alias_item));
+ item->pr = prob_array[i] * (float)alias->n; /* SCALED */
+ item->idx = i;
+
+ if (item->pr<1.0) {
+ list_add(&small, item);
+ } else {
+ list_add(&large, item);
+ }
+ }
+
+ /*
+ * Build the probability and alias tables.
+ *
+ * Depending on the probability distribution,
+ * 'small' or 'large' may empty first, so we
+ * treat them together, while it lasts, and
+ * then separately.
+ */
+
+ for (i=0; !list_empty(&small) && !list_empty(&large); i++) {
+ /*
+ * Small and large not empty
+ */
+ sm = (struct alias_item *)list_head(&small);
+ lg = (struct alias_item *)list_head(&large);
+
+ list_remove_from(&small, sm);
+ list_remove_from(&large, lg);
+
+ coin->pr[sm->idx] = sm->pr;
+ coin->al[sm->idx] = lg->idx;
+
+ lg->pr = (lg->pr + sm->pr) - 1.0;
+
+ if (lg->pr<1.0) {
+ list_add(&small, lg);
+ } else {
+ list_add(&large, lg);
+ }
+
+ free(sm);
+ }
+
+ for (; !list_empty(&large); i++) {
+ /*
+ * Large is not empty
+ */
+ lg = (struct alias_item *)list_head(&large);
+ list_remove_from(&large, lg);
+
+ coin->pr[lg->idx] = 1.0;
+
+ free(lg);
+ }
+
+ for (; !list_empty(&small); i++) {
+ /*
+ * Small is not empty
+ */
+ sm = (struct alias_item *)list_head(&small);
+ list_remove_from(&small, sm);
+
+ coin->pr[sm->idx] = 1.0;
+
+ free(sm);
+ }
+
+ /*
+ * All worklist memory should be free'd
+ */
+ return;
+}
+
+
+/**
+ * coin_flip()
+ * -----------
+ * Sample from a discrete probability distribution.
+ *
+ * @coin : Pointer to alias structure.
+ * Return: Random integer value from (0,1,2,...,@alias->n).
+ */
+int coin_flip(struct coin_t *coin)
+{
+ int value;
+ int i;
+
+ for (i=0; i<coin->n; i++) {
+ value = roll_die(coin->n);
+
+ if (flip_coin(coin->pr[value], HEADS) == HEADS) {
+ return value;
+ } else {
+ return (int)coin->al[value];
+ }
+ }
+
+ return (coin->n+1); /* impossible */
+}
+
+/******************************************************************************
+ * FAIR AND BIASED 2-COIN FUNCTIONS
+ ******************************************************************************/
+
+/**
+ * flip_coin()
+ * -----------
+ * Simulate a biased coin flip, applying bias to a given side.
+ *
+ * @bias : Probability of @side appearing
+ * @side : 0 (heads) or 1 (tails).
+ * Return: 0 (heads) or 1 (tails).
+ *
+ * NOTE
+ * To flip a coin with arbitrary-valued sides (rather than boolean),
+ * see the 'coin_flip_value' macro, defined in coin.h.
+ */
+int flip_coin(double bias, int side)
+{
+ if (side == HEADS) {
+ return (prng_uniform_random()<bias) ? HEADS : TAILS;
+ } else if (side == TAILS) {
+ return (prng_uniform_random()<bias) ? TAILS : HEADS;
+ } else {
+ /* Error-ish */
+ return side;
+ }
+}
+
+
+/**
+ * flip_heads()
+ * ------------
+ * Simulate a biased coin flip, applying the bias to heads.
+ *
+ * @bias : Probability of 0 (heads) appearing
+ * Return: 0 (heads) or 1 (tails).
+ */
+int flip_heads(double bias)
+{
+ return (prng_uniform_random()<bias) ? HEADS : TAILS;
+}
+
+
+/**
+ * flip_tails()
+ * ------------
+ * Simulate a biased coin flip, applying the bias to tails.
+ *
+ * @bias : Probability of 1 (tails) appearing
+ * Return: 0 (heads) or 1 (tails).
+ */
+int flip_tails(double bias)
+{
+ return (prng_uniform_random()<bias) ? TAILS : HEADS;
+}
+
+
+/**
+ * flip_fair()
+ * -----------
+ * Simulate a fair flip of a 2-coin. Should feel familiar.
+ * Return: 0 (HEADS) or 1 (TAILS)
+ */
+int flip_fair(void)
+{
+ return (prng_uniform_random()<0.5) ? HEADS : TAILS;
+}
+
+/******************************************************************************
+ * FAIR M-COIN FUNCTIONS
+ ******************************************************************************/
+
+/**
+ * roll_die()
+ * ----------
+ * Simulate a fair roll of an n-sided die.
+ *
+ * @sides: Number of sides (of the die)
+ * Return: integer between 0 and (n-1).
+ */
+int roll_die(int sides)
+{
+ return (int)(floor(prng_uniform_random()*(double)sides))%sides;
+}
+
+
+/**
+ * roll_dice()
+ * -----------
+ * Simulate a fair roll of an n-sided die.
+ *
+ * @ndice: Number of dice
+ * @sides: Number of sides (of the die)
+ * Return: Sum of rolls
+ */
+int roll_dice(int ndice, int sides)
+{
+ int sum = 0;
+
+ while (ndice-->0) {
+ sum += roll_die(sides);
+ }
+
+ return sum;
+}
+
diff --git a/prng/coin2.h b/prng/coin2.h
new file mode 100644
index 0000000..ccc2750
--- /dev/null
+++ b/prng/coin2.h
@@ -0,0 +1,34 @@
+#ifndef __MATH_COIN_H
+#define __MATH_COIN_H
+
+/* Useful macros */
+#define HEADS 0
+#define TAILS 1
+
+/******************************************************************************
+ * GENERAL M-COIN
+ ******************************************************************************/
+struct coin_t {
+ int n; /* Number of sides */
+ float *pr; /* Length n probability array pr[i]=prob of side i */
+ float *al; /* Alias table, used in Vose-Walker Alias method */
+};
+
+void coin_setup(struct coin_t *coin, int num_sides, float *prob_array);
+int coin_flip (struct coin_t *coin);
+
+/******************************************************************************
+ * FAIR AND BIASED 2-COIN
+ ******************************************************************************/
+int flip_coin (double bias, int side);
+int flip_heads(double bias);
+int flip_tails(double bias);
+int flip_fair (void);
+
+/******************************************************************************
+ * FAIR M-COIN
+ ******************************************************************************/
+int roll_die (int num_sides);
+int roll_dice (int num_dice, int num_sides);
+
+#endif
diff --git a/prng/dice.c b/prng/dice.c
new file mode 100644
index 0000000..bec9726
--- /dev/null
+++ b/prng/dice.c
@@ -0,0 +1,40 @@
+#include <stdlib.h>
+#include <math.h>
+#include "prng.h"
+#include "dice.h"
+#include "alias.h"
+
+/**
+ * dice_roll()
+ * -----------
+ * Simulate a fair roll of an n-sided die.
+ *
+ * @sides: Number of sides (of the die)
+ * Return: integer between 0 and (n-1).
+ */
+int dice_roll(int sides)
+{
+ return (int)(floor(prng_uniform_random()*(double)sides))%sides;
+}
+
+
+/**
+ * dice_roll_multi()
+ * -----------------
+ * Simulate a fair roll of an n-sided die.
+ *
+ * @ndice: Number of dice
+ * @sides: Number of sides (of the die)
+ * Return: Sum of rolls
+ */
+int dice_roll_multi(int ndice, int sides)
+{
+ int sum = 0;
+
+ while (ndice-->0) {
+ sum += dice_roll(sides);
+ }
+
+ return sum;
+}
+
diff --git a/prng/dice.h b/prng/dice.h
new file mode 100644
index 0000000..b6019b6
--- /dev/null
+++ b/prng/dice.h
@@ -0,0 +1,8 @@
+#ifndef __MATH_DICE_H
+#define __MATH_DICE_H
+
+int dice_roll (int sides);
+int dice_roll_multi(int ndice, int sides);
+
+#endif
+
diff --git a/prng/mersenne.c b/prng/mersenne.c
new file mode 100644
index 0000000..c0e140c
--- /dev/null
+++ b/prng/mersenne.c
@@ -0,0 +1,168 @@
+/*
+ * mersenne.c -- "Mersenne twister" pseudorandom number generator
+ *
+ * Copyright (C) 2012 Jason Linehan
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+/******************************************************************************
+ * Implements a pseudorandom number generator using the Mersenne twister
+ * algorithm.
+ *
+ * Once initialized, the function mt_random() will return a random unsigned
+ * 32-bit integer, which may take any of the values in the interval [0,2^32)
+ *
+ * Profiling has indicated this algorithm to be about twice as fast as the
+ * library rand() function.
+ *
+ * The original code on which this implementation is based was written by
+ * Michael Brundage and has been placed in the public domain.
+ ******************************************************************************/
+#include <stdlib.h>
+#include <stdint.h>
+#include <time.h>
+#include "mersenne.h"
+
+#define MT_LEN 624
+#define MT_MAX 2147483648 /* Maximum value of mt output */
+
+#define MT_IA 397
+#define MT_IB (MT_LEN - MT_IA)
+
+#define UPPER_MASK 0x80000000UL /* Most significant w-r bits */
+#define LOWER_MASK 0x7FFFFFFFUL /* Least significant r bits */
+#define MATRIX_A 0x9908B0DFUL /* Constant vector A */
+
+#define TWIST(b,i,j) ((b)[i] & UPPER_MASK) | ((b)[j] & LOWER_MASK)
+
+ /* Magic is (x * MATRIX_A) for x={0,1} */
+#define MAGIC(s) ((((s)&0x1UL) == 0) ? 0x0UL : MATRIX_A)
+
+
+
+/**
+ * mt_initialize()
+ * ---------------
+ * Primes the generator with pseudo-random values.
+
+ * @mt : The mersenne twister struct
+ * Return: nothing
+ *
+ * NOTE
+ * The twister uses its own values to generate subsequent
+ * value sets. However, the first time it runs, the value
+ * array is empty, so we must be initialize it with some
+ * pseudo-random values to prime the pump.
+ */
+void mt_initialize(struct mt_t *mt)
+{
+ int i;
+
+ /* Set the pseudo-random seed value */
+ srand(time(NULL));
+
+ for (i=0; i<MT_LEN; i++) {
+ mt->value[i] = rand(); /* feed me */
+ }
+
+ /* never again */
+ mt->is_initialized = 1;
+}
+
+
+/**
+ * mt_random()
+ * -----------
+ * Generate the pseudo-random value
+ *
+ * @mt : Mersenne twister struct
+ * Return: A pseudorandom unsigned long integer.
+ */
+unsigned long mt_random_int32(struct mt_t *mt)
+{
+ unsigned long *value;
+ unsigned long s;
+ int i;
+
+ value = mt->value;
+
+ /*
+ * If we have run out of values, we memoize the
+ * generator by generating the next MT_LEN
+ * values all at once.
+ */
+ if (mt->index == MT_LEN || !mt->is_initialized) {
+ if (!(mt->is_initialized)) {
+ /*
+ * We must prime the generator the
+ * first time it runs.
+ */
+ mt_initialize(mt);
+ }
+
+ for (i=0; i<MT_IB; i++) {
+ s = TWIST(value, i, i+1);
+ value[i] = value[i+MT_IA] ^ (s >> 1) ^ MAGIC(s);
+ }
+
+ for (; i<MT_LEN-1; i++) {
+ s = TWIST(value, i, i+1);
+ value[i] = value[i-MT_IB] ^ (s >> 1) ^ MAGIC(s);
+ }
+
+ s = TWIST(value, MT_LEN-1, 0);
+
+ value[MT_LEN-1] = value[MT_IA-1] ^ (s >> 1) ^ MAGIC(s);
+
+ /* Reset the index */
+ mt->index = 0;
+ }
+
+ /* Get the next value and increment the generator index. */
+ s = mt->value[mt->index++];
+
+ /* Temper the value before outputting. */
+ s ^= (s >> 11);
+ s ^= (s << 7) & 0x9d2c5680UL;
+ s ^= (s << 15) & 0xefc60000UL;
+ s ^= (s >> 18);
+
+ return s;
+}
+
+
+/* Generates random number on real interval [0,1] */
+double mt_random_real_0(struct mt_t *mt)
+{
+ /* Divide by 2^32 - 1 */
+ return mt_random_int32(mt) * (1.0/4294967295.0);
+}
+
+/* Generates random number on real interval [0,1) */
+double mt_random_real_1(struct mt_t *mt)
+{
+ /* Divide by 2^32 */
+ return mt_random_int32(mt) * (1.0/4294967296.0);
+}
+
+/* Generates random number on real interval (0,1) */
+double mt_random_real_2(struct mt_t *mt)
+{
+ /* Divide by 2^32 */
+ return (((double)mt_random_int32(mt)) + 0.5) * (1.0/4294967296.0);
+}
+
+
diff --git a/prng/mersenne.h b/prng/mersenne.h
new file mode 100644
index 0000000..ec1cee1
--- /dev/null
+++ b/prng/mersenne.h
@@ -0,0 +1,23 @@
+#ifndef __MERSENNE_H
+#define __MERSENNE_H
+
+#include <math.h>
+
+#define MT_LEN 624
+
+struct mt_t {
+ /* The index of the generator */
+ int index;
+ /* If the generator has been initialized the first time. */
+ int is_initialized;
+ /* Stores the state of the generator */
+ unsigned long value[MT_LEN];
+};
+
+
+unsigned long mt_random_int32 (struct mt_t *mt);
+double mt_random_real_0(struct mt_t *mt); /* [0,1] */
+double mt_random_real_1(struct mt_t *mt); /* [0,1) */
+double mt_random_real_2(struct mt_t *mt); /* (0,1) */
+
+#endif
diff --git a/prng/prng.c b/prng/prng.c
new file mode 100644
index 0000000..a40599b
--- /dev/null
+++ b/prng/prng.c
@@ -0,0 +1,46 @@
+#include <stdlib.h>
+#include <math.h>
+#include "mersenne.h"
+
+/*
+ * Generator structure to generate random
+ * values via the Mersenne Twister
+ * algorithm.
+ */
+static struct mt_t Generator = {0};
+
+
+/**
+ * prng_uniform_random()
+ * `````````````````````
+ * Generate a uniform random value on [0,1]
+ * Return: random real
+ */
+double prng_uniform_random(void)
+{
+ return mt_random_real_0(&Generator);
+}
+
+/**
+ * prng_uniform_random_open_right()
+ * ````````````````````````````````
+ * Generate a uniform random value on [0,1)
+ * Return: random real
+ */
+double prng_uniform_random_open_right(void)
+{
+ return mt_random_real_1(&Generator);
+}
+
+/**
+ * prng_uniform_random_open()
+ * ``````````````````````````
+ * Generate a uniform random value on (0,1)
+ * Return: random real
+ */
+double prng_uniform_random_open(void)
+{
+ return mt_random_real_2(&Generator);
+}
+
+
diff --git a/prng/prng.h b/prng/prng.h
new file mode 100644
index 0000000..256b8cf
--- /dev/null
+++ b/prng/prng.h
@@ -0,0 +1,10 @@
+#ifndef __JDL_PRNG
+#define __JDL_PRNG
+
+#include "mersenne.h"
+
+double prng_uniform_random(void);
+double prng_uniform_random_open_right(void);
+double prng_uniform_random_open(void);
+
+#endif
diff --git a/prng/tbl/5tbl.c b/prng/tbl/5tbl.c
new file mode 100644
index 0000000..6b775e8
--- /dev/null
+++ b/prng/tbl/5tbl.c
@@ -0,0 +1,183 @@
+/*
+This main asks for the Poisson parameter lambda, then generates 10^8
+random Poisson-lambda variates using the compacted 5-table method.
+It then does a goodness-of-fit test on the output.
+It then asks for the binomial parameters n,p, generates 10^8
+binomial(n,p) variates with the compacted 5-table method and tests
+the output. Generation speed is 10-15 nanoseconds for each variate,
+some 70-100 million per second.
+Output is displayed on the screen and also sent to the file tests.out.
+*/
+
+
+
+#include <stdio.h>
+#include <math.h>
+#include <stdlib.h>
+
+# define dg(m,k) ((m>>(30-6*k))&63) /* gets kth digit of m (base 64) */
+
+static int *P; /* Probabilities as an array of 30-bit integers*/
+static int size,t1,t2,t3,t4,offset,last; /* size of P[], limits for table lookups */
+static unsigned long jxr=182736531; /* Xorshift RNG */
+static short int *AA,*BB,*CC,*DD,*EE; /* Tables for condensed table-lookup */
+static FILE *fp;
+
+void get5tbls() /* Creates the 5 tables after array P[size] has been created */
+{ int i,j,k,m,na=0,nb=0,nc=0,nd=0,ne=0;
+ /* get table sizes, malloc */
+for(i=0;i<size;i++){m=P[i];na+=dg(m,1);nb+=dg(m,2);nc+=dg(m,3);nd+=dg(m,4);ne+=dg(m,5);}
+AA=malloc(na*sizeof(int));BB=malloc(nb*sizeof(int));
+CC=malloc(nc*sizeof(int));DD=malloc(nd*sizeof(int));EE=malloc(ne*sizeof(int));
+printf(" Table sizes:%d,%d,%d,%d,%d,total=%d\n",na,nb,nc,nd,ne,na+nb+nc+nd+ne);
+fprintf(fp," Table sizes:%d,%d,%d,%d,%d,total=%d\n",na,nb,nc,nd,ne,na+nb+nc+nd+ne);
+t1=na<<24; t2=t1+(nb<<18); t3=t2+(nc<<12); t4=t3+(nd<<6);
+na=nb=nc=nd=ne=0;
+ /* Fill tables AA,BB,CC,DD,EE */
+for(i=0;i<size;i++){m=P[i]; k=i+offset;
+ for(j=0;j<dg(m,1);j++) AA[na+j]=k; na+=dg(m,1);
+ for(j=0;j<dg(m,2);j++) BB[nb+j]=k; nb+=dg(m,2);
+ for(j=0;j<dg(m,3);j++) CC[nc+j]=k; nc+=dg(m,3);
+ for(j=0;j<dg(m,4);j++) DD[nd+j]=k; nd+=dg(m,4);
+ for(j=0;j<dg(m,5);j++) EE[ne+j]=k; ne+=dg(m,5);
+ }
+}// end get5tbls
+
+void PoissP(double lam) /* Creates Poisson Probabilites */
+{int i,j=-1,nlam; /* P's are 30-bit integers, assumed denominator 2^30 */
+double p=1,c,t=1.;
+ /* generate P's from 0 if lam<21.4 */
+if(lam<21.4){p=t=exp(-lam); for(i=1;t*2147483648.>1;i++) t*=(lam/i);
+ size=i-1; last=i-2;
+ /* Given size, malloc and fill P array, (30-bit integers) */
+ P=malloc(size*sizeof(int)); P[0]=exp(-lam)*(1<<30)+.5;
+ for(i=1;i<size;i++) {p*=(lam/i); P[i]=p*(1<<30)+.5; }
+ }
+ /* If lam>=21.4, generate from largest P up,then largest down */
+if(lam>=21.4)
+{nlam=lam; /*first find size */
+c=lam*exp(-lam/nlam);
+for(i=1;i<=nlam;i++) t*=(c/i);
+p=t;
+for(i=nlam+1;t*2147483648.>1;i++) t*=(lam/i);
+last=i-2;
+t=p; j=-1;
+for(i=nlam-1;i>=0;i--){t*=((i+1)/lam); if(t*2147483648.<1){j=i;break;} }
+offset=j+1; size=last-offset+1;
+ /* malloc and fill P array as 30-bit integers */
+P=malloc(size*sizeof(int));
+t=p; P[nlam-offset]=p*(1<<30)+0.5;
+for(i=nlam+1;i<=last;i++){t*=(lam/i); P[i-offset]=t*(64<<24)+.5;}
+t=p;
+for(i=nlam-1;i>=offset;i--){t*=((i+1)/lam);P[i-offset]=t*(1<<30)+0.5;}
+} // end lam>=21.4
+} //end PoissP
+
+void BinomP(int n, double p) /*Creates Binomial Probabilities */
+ /* Note: if n large and p near 1, generate j=Binom(n,1-p), return n-j*/
+{double h,t,p0;
+ int i,kb=0,ke,k=0,s;
+/* first find size of P array */
+p0=t=exp(n*log(1.-p));
+h=p/(1.-p);
+ke=n;
+if(t*2147483648.>1) {k=1;kb=0;}
+for(i=1;i<=n;i++)
+ { t*=(n+1-i)*h/i;
+ if(k==0 && t*2147483648.>1) {k=1;kb=i;}
+ if(k==1 && t*2147483648.<1) {k=2;ke=i-1;}
+ }
+ size=ke-kb+1; offset=kb;
+ /* Then malloc and assign P values as 30-bit integers */
+P=malloc(size*sizeof(int));
+t=p0; for(i=1;i<=kb;i++) t*=(n+1-i)*h/i;
+ s=t*1073741824.+.5; P[0]=s;
+for(i=kb+1;i<=ke;i++) {t*=(n+1-i)*h/i;
+ P[i-kb]=t*1073741824.+.5;
+ s+=P[i-kb];}
+i=n*p-offset; P[i]+=(s-1073741824);
+} //end BinomP
+
+
+
+/* Discrete random variable generating function */
+
+int Dran() /* Uses 5 compact tables */
+{unsigned long j;
+ jxr^=jxr<<13; jxr^=jxr>>17; jxr^=jxr<<5; j=(jxr>>2);
+if(j<t1) return AA[j>>24];
+if(j<t2) return BB[(j-t1)>>18];
+if(j<t3) return CC[(j-t2)>>12];
+if(j<t4) return DD[(j-t3)>>6];
+return EE[j-t4];
+}
+
+
+double Phi(double x)
+ {long double s=x,t=0,b=x,q=x*x,i=1;
+ while(s!=t) s=(t=s)+(b*=q/(i+=2));
+ return .5+s*exp(-.5*q-.91893853320467274178L);
+ }
+
+double chisq(double z,int n)
+{double s=0.,t=1.,q,h;
+ int i;
+ if(z<=0.) return (0.);
+ if(n>3000) return Phi((exp(log(z/n)/3.)-1.+2./(9*n))/sqrt(2./(9*n)));
+ h=.5*z;
+ if(n&1){q=sqrt(z); t=2*exp(-.5*z-.918938533204673)/q;
+ for(i=1;i<=(n-2);i+=2){ t=t*z/i; s+=t;}
+ return(2*Phi(q)-1-s);
+ }
+ for(i=1;i<n/2;i++) { t=t*h/i; s+=t;}
+ return (1.-(1+s)*exp(-h));
+}
+
+
+void Dtest(n)
+/* requires static 'size', static int array P[size] */
+/* generates n Dran()'s, tests output */
+{ double x=0,y,s=0,*E;
+ int kb=0,ke=1000,i,j=0,*M;
+ E=malloc(size*sizeof(double));
+ M=malloc(size*sizeof(int));
+ for(i=0;i<size;i++) {E[i]=(n+0.)*P[i]/1073741824.;M[i]=0;}
+ s=0; for(i=0;i<size;i++) {s+=E[i]; if(s>10){kb=i;E[kb]=s;break;} }
+ s=0; for(i=size-1;i>0;i--) {s+=E[i]; if(s>10){ke=i;E[ke]=s;break;} }
+ j=0; for(i=0;i<=kb;i++) j+=M[i]; M[kb]=j;
+ j=0; for(i=ke;i<size;i++) j+=M[i]; M[ke]=j;
+ for(i=0;i<size;i++) M[i]=0; s=0; x=0;
+ for(i=0;i<n;i++) {j=Dran(); if(j<kb+offset) j=kb+offset;
+ if(j>ke+offset) j=ke+offset;
+ M[j-offset]++;}
+printf("\n D Observed Expected (O-E)^2/E sum\n");
+fprintf(fp,"\n D Observed Expected (O-E)^2/E sum\n");
+ for(i=kb;i<=ke;i++){ y=M[i]-E[i]; y=y*y/E[i]; s+=y;
+
+ printf("%4d %10d %12.2f %7.2f %7.2f\n",i+offset,M[i],E[i],y,s);
+ fprintf(fp,"%4d %10d %12.2f %7.2f %7.2f\n",i+offset,M[i],E[i],y,s);
+ }
+ printf(" chisquare for %d d.f.=%7.2f, p=%7.5f\n",ke-kb,s,chisq(s,ke-kb));
+ fprintf(fp," chisquare for %d d.f.=%7.2f, p=%7.5f\n",ke-kb,s,chisq(s,ke-kb));
+}
+
+int main(){
+int j=0,n,nsmpls=100000000;
+double lam,p;
+fp=fopen("tests.out","w");
+printf(" Enter lambda:\n");
+scanf("%lf",&lam);
+PoissP(lam); get5tbls();
+//printf("start"); for(n=0;n<1000000000;n++) j+=Dran(); printf("END\n"); //15 nanos
+Dtest(nsmpls);
+fprintf(fp," Above results for sample of %d from Poisson, lambda=%3.2f\n\n",nsmpls,lam);
+free(P);
+printf("\n Enter n and p for Binomial:\n");
+scanf("%d %lf",&n,&p);
+BinomP(n,p);get5tbls();
+Dtest(nsmpls);
+fprintf(fp," Above result for sample of %d from Binomial(%d,%3.3f)\n\n",nsmpls,n,p);
+free(P);
+printf(" Test results sent to file tests.out\n");
+return 0;
+}
diff --git a/prng/tbl/5tbl.c.zip b/prng/tbl/5tbl.c.zip
new file mode 100644
index 0000000..64ae760
--- /dev/null
+++ b/prng/tbl/5tbl.c.zip
Binary files differ
diff --git a/prng/tbl/__MACOSX/._5tbl.c b/prng/tbl/__MACOSX/._5tbl.c
new file mode 100644
index 0000000..e4ee46c
--- /dev/null
+++ b/prng/tbl/__MACOSX/._5tbl.c
Binary files differ
diff --git a/tree/.node_cost.c.swp b/tree/.node_cost.c.swp
new file mode 100644
index 0000000..72c1f20
--- /dev/null
+++ b/tree/.node_cost.c.swp
Binary files differ
diff --git a/tree/node_alloc.c b/tree/node_alloc.c
new file mode 100644
index 0000000..a13ef9f
--- /dev/null
+++ b/tree/node_alloc.c
@@ -0,0 +1,87 @@
+#include "ytree.h"
+
+static int ID = 0;
+
+/******************************************************************************
+ * CREATE/COPY/DELETE
+ ******************************************************************************/
+
+/**
+ * ynode_create()
+ * --------------
+ * Allocate a node structure and set its value.
+ *
+ * @value: Value to create a node for.
+ * Return: Pointer to created node, or NULL;
+ */
+struct ynode_t *ynode_create(ynode_value_t value, ynode_value_t key)
+{
+ struct ynode_t *n;
+
+ n = calloc(1, sizeof(struct ynode_t));
+
+ n->value = value;
+ n->id = ID++;
+ n->key = key;
+
+ return n;
+}
+
+
+/**
+ * ynode_create_root()
+ * -------------------
+ * Allocate a node structure without a value. (usually the root)
+ *
+ * Return: Pointer to created node, or NULL;
+ */
+struct ynode_t *ynode_create_root(void)
+{
+ struct ynode_t *n;
+
+ n = ynode_create((ynode_value_t)0, 0);
+ n->P = NULL;
+
+ return n;
+}
+
+
+/**
+ * ynode_copy()
+ * ------------
+ * Copy a tree node and its attributes, but not its relations.
+ *
+ * @n : Pointer to tree node (will be copied).
+ * Return: Pointer to tree node (note that relations will be NULL).
+ */
+struct ynode_t *ynode_copy(struct ynode_t *n)
+{
+ struct ynode_t *copy;
+
+ copy = ynode_create(n->value, n->key);
+
+ copy->id = n->id;
+
+ copy->L = NULL;
+ copy->R = NULL;
+ copy->P = NULL;
+
+ return copy;
+}
+
+
+/**
+ * ynode_destroy()
+ * ---------------
+ * Free memory associated with a node.
+ *
+ * @n : Pointer to tree node.
+ * Return: nothing
+ */
+void ynode_destroy(struct ynode_t *n)
+{
+ if (n != NULL) {
+ free(n);
+ }
+ return;
+}
diff --git a/tree/node_check.c b/tree/node_check.c
new file mode 100644
index 0000000..7ad29ad
--- /dev/null
+++ b/tree/node_check.c
@@ -0,0 +1,204 @@
+#include "ytree.h"
+
+/******************************************************************************
+ * NODE PROPERTIES
+ ******************************************************************************/
+
+/**
+ * ynode_is_subtree_of()
+ * ---------------------
+ * Determine whether a node is a subtree of another.
+ *
+ * @a : Node to be checked
+ * @b : Node to be checked
+ * Return: 1 (TRUE), or 0 (FALSE)
+ */
+int ynode_is_subtree_of(struct ynode_t *a, struct ynode_t *b)
+{
+ if (a->P == NULL) {
+ /* The root is a subtree of nothing*/
+ return 0;
+ }
+ if (b->P == NULL) {
+ /* Everything is a subtree of the root */
+ return 1;
+ }
+
+ while (a->P != NULL) {
+ if (a->id == b->id) {
+ return 1;
+ } else {
+ a = a->P;
+ }
+ }
+ return 0;
+}
+
+
+/**
+ * ynode_is_disjoint()
+ * -------------------
+ * Determine whether two nodes are disjoint.
+ *
+ * @a : Node to be checked
+ * @b : Node to be checked
+ * Return: 1 (TRUE), or 0 (FALSE)
+ */
+int ynode_is_disjoint(struct ynode_t *a, struct ynode_t *b)
+{
+ return (!ynode_is_subtree_of(a, b) && !ynode_is_subtree_of(b, a));
+}
+
+
+/**
+ * ynode_is_sibling()
+ * ------------------
+ * Determine whether two nodes are siblings.
+ *
+ * @a : Node to be checked
+ * @b : Node to be checked
+ * Return: 1 (TRUE), or 0 (FALSE)
+ */
+int ynode_is_sibling(struct ynode_t *a, struct ynode_t *b)
+{
+ return (a && b && a->P && b->P && a->P->id == b->P->id);
+}
+
+
+/**
+ * ynode_is_equal()
+ * ----------------
+ * Determine whether two nodes are equal.
+ *
+ * @a : Node to be checked
+ * @b : Node to be checked
+ * Return: 1 (TRUE), or 0 (FALSE)
+ */
+int ynode_is_equal(struct ynode_t *a, struct ynode_t *b)
+{
+ return (a && b && a->id == b->id);
+}
+
+
+/**
+ * ynode_is_left_child()
+ * ---------------------
+ * Determine whether a node is a left-child.
+ *
+ * @a : Node to be checked
+ * Return: 1 (TRUE), or 0 (FALSE)
+ */
+int ynode_is_left_child(struct ynode_t *a)
+{
+ return (a && a->P && a->P->L && a->P->L->id == a->id);
+}
+
+
+/**
+ * ynode_is_right_child()
+ * ----------------------
+ * Determine whether a node is a right-child.
+ *
+ * @a : Node to be checked
+ * Return: 1 (TRUE), or 0 (FALSE)
+ */
+int ynode_is_right_child(struct ynode_t *a)
+{
+ return (a && a->P && a->P->R && a->P->R->id == a->id);
+}
+
+
+/**
+ * ynode_is_full()
+ * ---------------
+ * Determine whether an (internal) node is full.
+ *
+ * @a : Node to be checked
+ * Return: 1 (TRUE), or 0 (FALSE)
+ */
+int ynode_is_full(struct ynode_t *a)
+{
+ return (a && a->L != NULL && a->R != NULL);
+}
+
+
+/**
+ * ynode_is_internal()
+ * -------------------
+ * Determine whether a node is internal.
+ *
+ * @a : Node to be checked
+ * Return: 1 (TRUE), or 0 (FALSE)
+ */
+int ynode_is_internal(struct ynode_t *a)
+{
+ /*return (a && (a->L != NULL || a->R != NULL));*/
+ return (a && a->value == YTREE_INTERNAL_NODE_LABEL);
+}
+
+
+/**
+ * ynode_is_leaf()
+ * ---------------
+ * Determine whether a node is a leaf.
+ *
+ * @a : Node to be checked
+ * Return: 1 (TRUE), or 0 (FALSE)
+ */
+int ynode_is_leaf(struct ynode_t *a)
+{
+ return (a && a->value != YTREE_INTERNAL_NODE_LABEL && a->L == NULL && a->R == NULL);
+}
+
+
+/**
+ * ynode_is_root()
+ * ---------------
+ * Determine whether a node is the root.
+ *
+ * @a : Node to be checked
+ * Return: 1 (TRUE), or 0 (FALSE)
+ */
+int ynode_is_root(struct ynode_t *a)
+{
+ return (a && a->P == NULL);
+}
+
+
+/**
+ * ynode_is_ternary()
+ * ------------------
+ * Verify the shape properties of a ternary tree.
+ *
+ * @n : Node (pseudo-root) to check from.
+ * Return: 1 (TRUE) or 0 (FALSE).
+ */
+int ynode_is_ternary(struct ynode_t *n)
+{
+ if (n == NULL) {
+ return 0;
+ }
+
+ int lc = ynode_count_leaves(n);
+ int ic = ynode_count_internal(n);
+
+ /*
+ * For n items stored in the tree, there
+ * should be n leaf nodes, and (n-2)
+ * internal nodes.
+ *
+ * Unless of course there are <2 items
+ * in the tree, in which case no new
+ * nodes will be created, excepting the
+ * root.
+ */
+ if ((lc<=2 && ic==0) || ic==(lc-2)) {
+ return 1;
+ } else {
+ #ifdef YTREE_PRINT_DEBUG
+ fprintf(stderr, "leaves:%d internal:%d\n", lc, ic);
+ #endif
+ return 0;
+ }
+}
+
diff --git a/tree/node_cost.c b/tree/node_cost.c
new file mode 100644
index 0000000..4004031
--- /dev/null
+++ b/tree/node_cost.c
@@ -0,0 +1,257 @@
+#include "ytree.h"
+
+/******************************************************************************
+ * NODE COST
+ ******************************************************************************/
+
+/**
+ * ynode_get_cost()
+ * ----------------
+ * Determine the cost value of a particular node.
+ *
+ * @n : Pointer to node.
+ * @distance: Distance matrix from which to compute the cost.
+ * Return : Cost value of @n.
+ */
+float ynode_get_cost(struct ynode_t *n, float **distance)
+{
+ struct ynode_t *root;
+
+ ynode_value_t *value_L;
+ ynode_value_t *value_R;
+ ynode_value_t *value_P;
+
+ int count_L;
+ int count_R;
+ int count_P;
+
+ int combi_L;
+ int combi_R;
+ int combi_P;
+
+ float distance_LR = 0.0;
+ float distance_PL = 0.0;
+ float distance_PR = 0.0;
+
+ float cost_LR;
+ float cost_PL;
+ float cost_PR;
+
+ float cost = 0.0;
+
+ int i;
+ int j;
+
+ if (n == NULL) {
+ return -1.0;
+ }
+
+ root = ynode_get_root(n);
+
+ count_L = ynode_count_leaves(n->L);
+ count_R = ynode_count_leaves(n->R);
+ count_P = ynode_count_leaves_not(root, n);
+
+ combi_L = binomial(count_L, 2);
+ combi_R = binomial(count_R, 2);
+ combi_P = binomial(count_P, 2);
+
+ value_L = ynode_get_values(n->L);
+ value_R = ynode_get_values(n->R);
+ value_P = ynode_get_values_not(root, n);
+
+ for (i=0; i<count_L; i++) {
+ for (j=0; j<count_R; j++) {
+ if (value_L[i] == value_R[j]) {
+ printf("CHAOS REIGNS\n");
+ abort();
+ }
+ if (value_L[i]<=DATA_COUNT && value_R[j]<=DATA_COUNT) {
+ distance_LR += distance[value_L[i]][value_R[j]];
+ }
+ }
+ }
+
+ for (i=0; i<count_P; i++) {
+ for (j=0; j<count_L; j++) {
+ if (value_P[i] == value_L[j]) {
+ printf("CHAOS REIGNS\n");
+ abort();
+ }
+ if (value_P[i]<=DATA_COUNT && value_L[j]<=DATA_COUNT) {
+ distance_PL += distance[value_P[i]][value_L[j]];
+ }
+ }
+ }
+
+ for (i=0; i<count_P; i++) {
+ for (j=0; j<count_R; j++) {
+ if (value_P[i] == value_R[j]) {
+ printf("CHAOS REIGNS\n");
+ abort();
+ }
+ if (value_P[i]<=DATA_COUNT && value_R[j]<=DATA_COUNT) {
+ distance_PR += distance[value_P[i]][value_R[j]];
+ }
+ }
+ }
+
+ cost_LR = (float)combi_P * distance_LR;
+ cost_PL = (float)combi_R * distance_PL;
+ cost_PR = (float)combi_L * distance_PR;
+
+ cost = cost_LR + cost_PL + cost_PR;
+
+ /*printf("COUNT l:%d r:%d p:%d\n", count_L, count_R, count_P);*/
+ /*printf("BINOM l:%d r:%d p:%d\n", combi_L, combi_R, combi_P);*/
+ /*printf("DIST l:%f r:%f p:%f\n", distance_PR, distance_PL, distance_LR);*/
+ /*printf("COST l:%f r:%f p:%f\n", cost_PR, cost_PL, cost_LR);*/
+ /*printf("TOTAL %f\n", cost);*/
+
+ if (value_L != NULL) {
+ free(value_L);
+ }
+ if (value_R != NULL) {
+ free(value_R);
+ }
+ if (value_P != NULL) {
+ free(value_P);
+ }
+
+ return cost;
+}
+
+
+/**
+ * ynode_get_cost_max()
+ * --------------------
+ * Determine the maximum cost value of a particular node.
+ *
+ * @a : Pointer to node.
+ * @d : Distance matrix from which to compute the cost.
+ * @n : Number of items, i.e. @d is an @nx@n matrix.
+ * Return: Maximum cost value of @n.
+ */
+float ynode_get_cost_max(struct ynode_t *a, float **d, int n)
+{
+ int i;
+ int j;
+ int k;
+ int l;
+
+ int count;
+ int combo;
+ float cost;
+
+ float ijkl; /* ij|kl quartet */
+ float ikjl; /* ik|jl quartet */
+ float iljk; /* il|jk quartet */
+
+ cost = 0.0;
+ count = 0;
+ combo = binomial(n, 4);
+
+ for (i=0; i<n; i++) {
+ for (j=(i+1); j<n; j++) {
+ for (k=(j+1); k<n; k++) {
+ for (l=(k+1); l<n; l++) {
+
+ ijkl = d[i][j] + d[k][l];
+ ikjl = d[i][k] + d[j][l];
+ iljk = d[i][l] + d[j][k];
+
+ /*printf("max(%f, %f, %f) := %f\n", ijkl, ikjl, iljk, max_float(3, ijkl, ikjl, iljk));*/
+
+ cost += max_float(3, ijkl, ikjl, iljk);
+
+ count++;
+ }
+ }
+ }
+ }
+
+ if (count != combo) {
+ fprintf(stderr, "Combination error.\n");
+ return -1.0;
+ }
+
+ /*printf("max cost:%f\n", cost);*/
+
+ return cost;
+}
+
+
+/**
+ * ynode_get_cost_min()
+ * --------------------
+ * Determine the minimum cost value of a particular node.
+ *
+ * @a : Pointer to node.
+ * @d : Distance matrix from which to compute the cost.
+ * @n : Number of items, i.e. @d is an @nx@n matrix.
+ * Return: Minimum cost value of @n.
+ */
+float ynode_get_cost_min(struct ynode_t *a, float **d, int n)
+{
+ int i;
+ int j;
+ int k;
+ int l;
+
+ int count;
+ int combo;
+ float cost;
+
+ float ijkl; /* ij|kl quartet */
+ float ikjl; /* ik|jl quartet */
+ float iljk; /* il|jk quartet */
+
+ cost = 0.0;
+ count = 0;
+ combo = binomial(n, 4);
+
+ for (i=0; i<n; i++) {
+ for (j=(i+1); j<n; j++) {
+ for (k=(j+1); k<n; k++) {
+ for (l=(k+1); l<n; l++) {
+
+ ijkl = d[i][j] + d[k][l];
+ ikjl = d[i][k] + d[j][l];
+ iljk = d[i][l] + d[j][k];
+
+ /*printf("min(%f, %f, %f) := %f\n", ijkl, ikjl, iljk, min_float(3, ijkl, ikjl, iljk));*/
+
+ cost += min_float(3, ijkl, ikjl, iljk);
+
+ count++;
+ }
+ }
+ }
+ }
+
+ if (count != combo) {
+ fprintf(stderr, "Combination error.\n");
+ return -1.0;
+ }
+
+ /*printf("min cost:%f\n", cost);*/
+
+ return cost;
+}
+
+
+/**
+ * ynode_get_cost_scaled()
+ * -----------------------
+ * Scale the measured cost using the max/min cost.
+ *
+ * @c : Measured (unscaled) cost
+ * @M : Maximum possible cost
+ * @m : Minimum possible cost
+ * Return: Normalized value on interval (0,1).
+ */
+float ynode_get_cost_scaled(float c, float M, float m)
+{
+ return (M - c) / (M - m);
+}
+
diff --git a/tree/node_count.c b/tree/node_count.c
new file mode 100644
index 0000000..8d8b8ab
--- /dev/null
+++ b/tree/node_count.c
@@ -0,0 +1,84 @@
+#include "ytree.h"
+
+/******************************************************************************
+ * COUNTS
+ ******************************************************************************/
+
+int Count = 0;
+struct ynode_t *Exclude;
+
+void __impl__ynode_count_leaves(struct ynode_t *n, int i)
+{
+ if (ynode_is_leaf(n) && n->value != YTREE_INTERNAL_NODE_LABEL) {
+ Count++;
+ }
+}
+
+void __impl__ynode_count_leaves_not(struct ynode_t *n, int i)
+{
+ if (ynode_is_leaf(n) && n->value != YTREE_INTERNAL_NODE_LABEL) {
+ if (!ynode_is_subtree_of(n, Exclude)) {
+ Count++;
+ }
+ }
+}
+
+void __impl__ynode_count_internal(struct ynode_t *n, int i)
+{
+ if (ynode_is_internal(n)) {
+ Count++;
+ }
+}
+
+
+/**
+ * ynode_count_leaves()
+ * --------------------
+ * Count the number of leaf nodes in a tree.
+ *
+ * @n : Node from which to count (pseudo-root).
+ * Return: number of leaf nodes
+ */
+int ynode_count_leaves(struct ynode_t *n)
+{
+ Count = 0;
+ ynode_traverse_inorder(n, __impl__ynode_count_leaves);
+
+ return Count;
+}
+
+
+/**
+ * ynode_count_leaves_not()
+ * ------------------------
+ * Count the number of leaf nodes in a tree, excluding some subtree.
+ *
+ * @n : Node from which to count (pseudo-root).
+ * @x : Node whose subtree will be excluded.
+ * Return: number of leaf nodes
+ */
+int ynode_count_leaves_not(struct ynode_t *n, struct ynode_t *x)
+{
+ Count = 0;
+ Exclude = x;
+ ynode_traverse_inorder(n, __impl__ynode_count_leaves_not);
+
+ return Count;
+}
+
+
+/**
+ * ynode_count_internal()
+ * ----------------------
+ * Count the number of internal nodes in a tree.
+ *
+ * @n : Node from which to count (pseudo-root).
+ * Return: number of internal nodes
+ */
+int ynode_count_internal(struct ynode_t *n)
+{
+ Count = 0;
+ ynode_traverse_inorder(n, __impl__ynode_count_internal);
+
+ return Count;
+}
diff --git a/tree/node_get.c b/tree/node_get.c
new file mode 100644
index 0000000..0d3db20
--- /dev/null
+++ b/tree/node_get.c
@@ -0,0 +1,265 @@
+#include "ytree.h"
+
+/******************************************************************************
+ * RANDOM NODES
+ ******************************************************************************/
+
+struct ynode_t *Random_node;
+
+void __impl__ynode_get_random(struct ynode_t *n, int i)
+{
+ /*
+ * The i-th node examined will be the
+ * chosen node with probability 1/i.
+ */
+ if (coin_flip(1/i, HEADS) == HEADS) {
+ Random_node = n;
+ }
+}
+
+void __impl__ynode_get_random_internal(struct ynode_t *n, int i)
+{
+ /*
+ * The i-th node examined will be the
+ * chosen node with probability 1/i.
+ */
+ if (ynode_is_internal(n)) {
+ if (coin_flip(1/i, HEADS) == HEADS) {
+ Random_node = n;
+ }
+ }
+}
+
+/**
+ * ynode_get_random()
+ * ------------------
+ * Select a random node from the tree.
+ *
+ * @n : Pointer to root or pseudo-root node
+ * Return: Pointer to some node under @n.
+ */
+struct ynode_t *ynode_get_random(struct ynode_t *n)
+{
+ Random_node = NULL;
+
+ if (n != NULL) {
+ ynode_traverse_inorder(n, __impl__ynode_get_random);
+ }
+
+ return Random_node;
+}
+
+/**
+ * ynode_get_random_internal()
+ * ---------------------------
+ * Select a random node from the tree.
+ *
+ * @n : Pointer to root or pseudo-root node
+ * Return: Pointer to some node under @n.
+ */
+struct ynode_t *ynode_get_random_internal(struct ynode_t *n)
+{
+ Random_node = NULL;
+
+ if (n != NULL) {
+ ynode_traverse_inorder(n, __impl__ynode_get_random_internal);
+ }
+
+ return Random_node;
+}
+
+/**
+ * ynode_get_random_leaf()
+ * -----------------------
+ * Select a random leaf node from the tree.
+ *
+ * @n : Pointer to root or pseudo-root node
+ * Return: Pointer to some leaf node under @n.
+ */
+struct ynode_t *ynode_get_random_leaf(struct ynode_t *n)
+{
+ if (n == NULL) {
+ return NULL;
+ }
+
+ if (ynode_is_leaf(n)) {
+ return n;
+ }
+
+ if (coin_fair()) {
+ return ynode_get_random_leaf(n->L);
+ } else {
+ return ynode_get_random_leaf(n->R);
+ }
+}
+
+
+
+/******************************************************************************
+ * GETTERS
+ ******************************************************************************/
+
+/**
+ * ynode_get_sibling()
+ * -------------------
+ * Select the sibling of a given node.
+ *
+ * @n : Pointer to a node
+ * Return: Pointer to the sibling of @n, or NULL.
+ */
+struct ynode_t *ynode_get_sibling(struct ynode_t *n)
+{
+ if (n != NULL) {
+ if (n->P == NULL) {
+ /* Has no siblings, since it has no parents */
+ return NULL;
+ }
+ if (n->P->L->id == n->id) {
+ return n->P->R;
+ }
+ if (n->P->R->id == n->id) {
+ return n->P->L;
+ }
+ }
+ return NULL;
+}
+
+
+/**
+ * ynode_get_root()
+ * ----------------
+ * Select the root of a given node.
+ *
+ * @n : Pointer to a node
+ * Return: Pointer to the sibling of @n, or NULL.
+ */
+struct ynode_t *ynode_get_root(struct ynode_t *n)
+{
+ if (n == NULL) {
+ return NULL;
+ } else {
+ while (n->P != NULL) {
+ n = n->P;
+ }
+ }
+ return n;
+}
+
+
+/**
+ * ynode_get_path()
+ * ----------------
+ * Create a path string for a given node.
+ *
+ * @n : Pointer to tree node.
+ * Return: {L,R}+ string representing path (position) in tree.
+ */
+char *ynode_get_path(struct ynode_t *a)
+{
+ struct ynode_t *ptr;
+ char *path;
+ int length = 0;
+ int i;
+
+ if (a == NULL) {
+ fprintf(stderr, "Cannot generate path for empty node\n");
+ return NULL;
+ }
+
+ ptr = a;
+
+ /* Distance to root will be the length of the path. */
+ while (ptr->P != NULL) {
+ length++;
+ ptr = ptr->P;
+ }
+
+ path = calloc(length+1, sizeof(char));
+ path[length] = 0;
+
+ ptr = a;
+ i = length;
+
+ while (ptr->P != NULL) {
+ if (i > 0) {
+ i--;
+ }
+ if (ynode_is_left_child(ptr)) {
+ path[i] = 'L';
+ } else {
+ path[i] = 'R';
+ }
+ ptr = ptr->P;
+ }
+
+ /*printf("%s (%d/%d)\n", path, strlen(path), length);*/
+
+ return path;
+}
+
+
+/******************************************************************************
+ * GET VALUES
+ ******************************************************************************/
+
+int Value_array_length;
+ynode_value_t *Value_array;
+struct ynode_t *Excluding;
+
+void __impl__ynode_get_values(struct ynode_t *n, int i)
+{
+ if (ynode_is_leaf(n)) {
+ Value_array[Value_array_length++] = n->value;
+ }
+ return;
+}
+
+void __impl__ynode_get_values_not(struct ynode_t *n, int i)
+{
+ if (ynode_is_leaf(n) && !ynode_is_subtree_of(n, Excluding)) {
+ Value_array[Value_array_length++] = n->value;
+ }
+ return;
+}
+
+
+/**
+ * ynode_get_values()
+ * ------------------
+ * Retreive an array of leaf values below the given node.
+ *
+ * @n : Pointer to a node in the tree.
+ * Return: Array of leaf values.
+ */
+ynode_value_t *ynode_get_values(struct ynode_t *n)
+{
+ Value_array_length = 0;
+ Value_array = calloc(DATA_COUNT, sizeof(ynode_value_t));
+
+ ynode_traverse_inorder(n, __impl__ynode_get_values);
+
+ /* Caller is responsible for freeing this. */
+ return Value_array;
+}
+
+
+/**
+ * ynode_get_values_not()
+ * ----------------------
+ * Retreive an array of leaf values below the given node, excluding subtree.
+ *
+ * @n : Pointer to a node in the tree.
+ * @x : Pointer to a node in the tree (subtree to exclude).
+ * Return: Array of leaf values.
+ */
+ynode_value_t *ynode_get_values_not(struct ynode_t *n, struct ynode_t *x)
+{
+ Value_array_length = 0;
+ Value_array = calloc(DATA_COUNT, sizeof(ynode_value_t));
+ Excluding = x;
+
+ ynode_traverse_inorder(n, __impl__ynode_get_values_not);
+
+ /* Caller is responsible for freeing this. */
+ return Value_array;
+}
diff --git a/tree/node_insert.c b/tree/node_insert.c
new file mode 100644
index 0000000..f65bcee
--- /dev/null
+++ b/tree/node_insert.c
@@ -0,0 +1,272 @@
+#include "ytree.h"
+
+/******************************************************************************
+ * NODE INSERTION HELPERS
+ ******************************************************************************/
+
+
+/**
+ * ynode_add_left()
+ * ----------------
+ * Insert value as left child of a node.
+ *
+ * @n : (Internal) node to insert at
+ * @value: Value to add at @n.
+ * Return: Pointer to the created node, or NULL.
+ *
+ * . .
+ * \ INSERT(1) / \
+ * 2 1 2
+ */
+struct ynode_t *ynode_add_left(struct ynode_t *n, ynode_value_t value, ynode_value_t key)
+{
+ n->L = ynode_create(value, key);
+ n->L->P = n;
+
+ return n->L;
+}
+
+
+/**
+ * ynode_add_right()
+ * -----------------
+ * Insert value as right child of a node.
+ *
+ * @n : (Internal) node to insert at
+ * @value: Value to add at @n.
+ * Return: Pointer to the created node, or NULL.
+ *
+ * . .
+ * / INSERT(2) / \
+ * 1 1 2
+ */
+struct ynode_t *ynode_add_right(struct ynode_t *n, ynode_value_t value, ynode_value_t key)
+{
+ n->R = ynode_create(value, key);
+ n->R->P = n;
+
+ return n->R;
+}
+
+
+/**
+ * ynode_add_left_level()
+ * ----------------------
+ * Insert value at leaf, and convert leaf to internal node.
+ *
+ * @n : (Leaf) node to insert [will be converted to internal node]
+ * @value: Value to add at @n.
+ * Return: Pointer to the created node, or NULL.
+ *
+ * . .
+ * / INSERT(1) /
+ * 2 .
+ * / \
+ * 1 2
+ */
+struct ynode_t *ynode_add_left_level(struct ynode_t *n, ynode_value_t value, ynode_value_t key)
+{
+ n->L = ynode_create(value, key);
+ n->R = ynode_create(n->value, n->key);
+ n->value = '.';
+
+ n->R->P = n;
+ n->L->P = n;
+
+ return n->L;
+}
+
+
+/**
+ * ynode_add_right_level()
+ * -----------------------
+ * Insert value at leaf, and convert leaf to internal node.
+ *
+ * @n : (Leaf) node to insert [will be converted to internal node]
+ * @value: Value to add at @n.
+ * Return: Pointer to the created node, or NULL.
+ *
+ * . .
+ * / INSERT(2) /
+ * 1 .
+ * / \
+ * 1 2
+ */
+struct ynode_t *ynode_add_right_level(struct ynode_t *n, ynode_value_t value, ynode_value_t key)
+{
+ n->R = ynode_create(value, key);
+ n->L = ynode_create(n->value, n->key);
+ n->value = '.';
+
+ n->R->P = n;
+ n->L->P = n;
+
+ return n->R;
+}
+
+
+/**
+ * ynode_add_before()
+ * ------------------
+ * Insert a node between the given node and it's parent.
+ *
+ * @n : Node before which the new node will be inserted.
+ * @value: Value to place at the new node.
+ * Return: Pointer to the created node, or NULL.
+ */
+struct ynode_t *ynode_add_before(struct ynode_t *a, ynode_value_t value, ynode_value_t key)
+{
+ struct ynode_t *par;
+ struct ynode_t *new;
+
+ /* Store the parent */
+ par = a->P;
+
+ /* Insert the new node between a and parent(a). */
+ new = ynode_create(value, key);
+
+ a->P = new;
+ new->P = par;
+
+ if (par->L && (par->L->id == a->id)) {
+ /* A was the left child */
+ new->L = a;
+ par->L = new;
+ } else {
+ /* A was the right child */
+ new->R = a;
+ par->R = new;
+ }
+
+ return new;
+}
+
+
+/******************************************************************************
+ * NODE INSERT
+ ******************************************************************************/
+
+/**
+ * ynode_insert()
+ * --------------
+ * Insert a value into the tree
+ *
+ * @n : Node (possibly root) under which to insert the value
+ * @value: Value to insert
+ * Return: Pointer to the created node, or NULL.
+ *
+ * NOTE
+ * If inserted at an internal node, the value will sink down
+ * the tree until an available leaf position can be found.
+ *
+ * In sinking the inserted value down to a suitable leaf, if
+ * we use the traditional practice of comparing @value with
+ * the value stored at each node in order to choose a path,
+ * there will be a problem.
+ *
+ * Internal nodes are labeled with INTERNAL_NODE_LABEL, a
+ * special value chosen to be disjoint from the input set,
+ * so that internal nodes can be distinguished from value
+ * nodes.
+ *
+ * This means that a all input values will always be sent in
+ * the same direction out of any predicate with an internal
+ * node, so the tree will always branch the same way.
+ *
+ * The solution is to randomly choose a path.
+ */
+struct ynode_t *ynode_insert(struct ynode_t *n, ynode_value_t value, ynode_value_t key)
+{
+ if (n == NULL) {
+ return NULL;
+ }
+
+ if (ynode_is_internal(n) || ynode_is_root(n)) {
+ if (ynode_is_full(n)) {
+ /*
+ * This is a full internal node, so we
+ * have to sink down the tree to find
+ * a leaf node.
+ */
+ /* Random sink. See NOTE */
+ if (coin_fair()) {
+ return ynode_insert(n->L, value, key);
+ } else {
+ return ynode_insert(n->R, value, key);
+ }
+ } else {
+ if (n->L == NULL && n->R == NULL) {
+ /*
+ * Both leaves are open. Choose a random
+ * leaf to fill.
+ *
+ * When the tree is empty, the root node
+ * will appear like this.
+ */
+ if (coin_fair()) {
+ return ynode_add_left(n, value, key);
+ } else {
+ return ynode_add_right(n, value, key);
+ }
+ } else if (n->L == NULL && n->R != NULL) {
+ /* Left leaf open */
+ return ynode_add_left(n, value, key);
+ } else if (n->L != NULL && n->R == NULL) {
+ /* Right leaf open */
+ return ynode_add_right(n, value, key);
+ }
+ }
+ } else {
+ /*
+ * This is a leaf node, which now must become
+ * an internal node.
+ */
+ /*if (value < n->value) {*/
+ if (coin_fair()) {
+ return ynode_add_left_level(n, value, key);
+ } else {
+ return ynode_add_right_level(n, value, key);
+ }
+ }
+
+ return NULL;
+}
+
+
+
+/**
+ * ynode_insert_on_path()
+ * ----------------------
+ * Insert a node according to a path string.
+ *
+ * @r : Pointer to root node of tree to insert into.
+ * @a : Pointer to tree node (will be inserted)
+ * @path : Path string {L,R}+ for @a to follow.
+ * Return: Nothing.
+ *
+ * NOTE
+ * For creating a path from a node, see ynode_get_path().
+ */
+void ynode_insert_on_path(struct ynode_t *r, struct ynode_t *a, char *path)
+{
+ int l = strlen(path);
+ int i;
+
+ struct ynode_t *node = r;
+
+ for (i=0; i<l-1; i++) {
+ if (path[i] == 'L') {
+ node = node->L;
+ } else {
+ node = node->R;
+ }
+ }
+
+ if (path[l-1] == 'L') {
+ node->L = a;
+ } else {
+ node->R = a;
+ }
+
+ a->P = node;
+}
diff --git a/tree/node_mutate.c b/tree/node_mutate.c
new file mode 100644
index 0000000..01aad8a
--- /dev/null
+++ b/tree/node_mutate.c
@@ -0,0 +1,322 @@
+#include "ytree.h"
+
+/******************************************************************************
+ * ynode MUTATION HELPERS
+ ******************************************************************************/
+
+/**
+ * ynode_contract()
+ * ----------------
+ * Contract a relation, destroying the child.
+ *
+ * @n : Parent node (child will be contracted)
+ * Return: Child node (can be destroyed)
+ *
+ * NOTE
+ * This operation takes place between a parent node
+ * and a child node. It destroys an un-necessary internal
+ * linkage.
+ *
+ * x x
+ * | / \
+ * y A B
+ * / \
+ * A B
+ */
+struct ynode_t *ynode_contract(struct ynode_t *n)
+{
+ struct ynode_t *to_remove;
+
+ if (n != NULL) {
+ if (n->L != NULL && n->R != NULL) {
+ printf("You cannot contract a full node and preserve the binary property.\n");
+ return NULL;
+ }
+ if (n->L == NULL && n->R == NULL) {
+ printf("There are no children to be contracted.\n");
+ return NULL;
+ }
+
+ /* Exactly one will be true */
+ if (n->L != NULL) {
+ to_remove = n->L;
+ n->L = n->L->L;
+ n->R = n->L->R;
+ } else {
+ to_remove = n->R;
+ n->L = n->R->L;
+ n->R = n->R->R;
+ }
+
+ if (n->L != NULL) {
+ n->L->P = n;
+ }
+ if (n->R != NULL) {
+ n->R->P = n;
+ }
+
+ return to_remove;
+ }
+ return NULL;
+}
+
+/**
+ * ynode_promote()
+ * ---------------
+ * Promote a child node, overwriting its parent.
+ *
+ * @n : Child node (parent will be destroyed)
+ * Return: Parent node (can be destroyed)
+ *
+ * NOTE
+ * This operation takes place between a parent node
+ * and a child node. It destroys an un-necessary internal
+ * linkage.
+ *
+ * x y
+ * | / \
+ * y A B
+ * / \
+ * A B
+ */
+struct ynode_t *ynode_promote(struct ynode_t *n)
+{
+ struct ynode_t *to_remove;
+
+ if (n != NULL) {
+ if (n->P == NULL) {
+ printf("You cannot promote the root.\n");
+ return NULL;
+ }
+ if (n->P->L != NULL && n->P->R != NULL) {
+ printf("You cannot promote to a full node and preserve the binary property.\n");
+ return NULL;
+ }
+ if (n->P->L == NULL && n->P->R == NULL) {
+ printf("Something has gone horribly wrong.\n");
+ return NULL;
+ }
+ if (n->P->P == NULL) {
+ printf("Cannot promote to the root. You must contract the root instead.\n");
+ return NULL;
+ }
+
+ to_remove = n->P;
+
+ if (ynode_is_left_child(n->P)) {
+ n->P->P->L = n;
+ } else {
+ n->P->P->R = n;
+ }
+
+ n->P = n->P->P;
+
+ return to_remove;
+ }
+ return NULL;
+}
+
+/******************************************************************************
+ * ynode MUTATION
+ ******************************************************************************/
+
+/**
+ * ynode_LEAF_INTERCHANGE()
+ * ------------------------
+ * Given two leaf nodes, exchange their position in the tree.
+ *
+ * @a : Pointer to (leaf) node
+ * @b : Pointer to (leaf) node
+ * Return: Nothing.
+ */
+void ynode_LEAF_INTERCHANGE(struct ynode_t *a, struct ynode_t *b)
+{
+ struct ynode_t *a_par;
+ struct ynode_t *b_par;
+
+ if (a != NULL && b != NULL) {
+
+ if (a->id == b->id) {
+ /* Swap will do nothing */
+ return;
+ }
+
+ if (a->P != NULL && b->P != NULL) {
+ a_par = a->P;
+ b_par = b->P;
+
+ if (ynode_is_left_child(a)) {
+ a_par->L = b;
+ } else {
+ a_par->R = b;
+ }
+
+ if (ynode_is_left_child(b)) {
+ b_par->L = a;
+ } else {
+ b_par->R = a;
+ }
+
+ a->P = b->P;
+ b->P = a_par;
+ }
+ }
+}
+
+
+/**
+ * ynode_SUBTREE_INTERCHANGE()
+ * ---------------------------
+ * Given two nodes, possibly leaves, exchange the subtrees they root.
+ *
+ * @a : Pointer to (leaf/internal) node
+ * @b : Pointer to (leaf/internal) node
+ * Return: Nothing.
+ */
+void ynode_SUBTREE_INTERCHANGE(struct ynode_t *a, struct ynode_t *b)
+{
+ struct ynode_t *a_parent;
+ struct ynode_t *b_parent;
+
+ if (a != NULL && b != NULL) {
+ if (ynode_is_equal(a, b)) {
+ /*printf("Identical subtrees. Interchange declined\n");*/
+ return;
+ }
+ /* The root node cannot participate */
+ if (ynode_is_root(a) || ynode_is_root(b)) {
+ /*printf("Cannot interchange root. Interchange declined\n");*/
+ return;
+ } else {
+
+ if (ynode_is_sibling(a, b)) {
+ /*printf("Sibling nodes.\n");*/
+ if (ynode_is_left_child(a)) {
+ a->P->L = b;
+ a->P->R = a;
+ } else {
+ a->P->L = a;
+ a->P->R = b;
+ }
+ return;
+ }
+
+ if (!ynode_is_disjoint(a, b)) {
+ /*printf("Non-disjoint subtrees. Interchange declined\n");*/
+ return;
+ }
+
+ a_parent = a->P;
+ b_parent = b->P;
+
+ if (ynode_is_left_child(a)) {
+ a_parent->L = b;
+ } else {
+ a_parent->R = b;
+ }
+
+ if (ynode_is_left_child(b)) {
+ b_parent->L = a;
+ } else {
+ b_parent->R = a;
+ }
+
+ a->P = b_parent;
+ b->P = a_parent;
+ }
+ }
+}
+
+
+/**
+ * ynode_SUBTREE_TRANSFER()
+ * ------------------------
+ * Given two nodes, possibly leaves, graft the subtree of one at another.
+ *
+ * @a : Pointer to (leaf/internal) node
+ * @b : Pointer to (leaf/internal) node
+ * Return: Nothing.
+ */
+void ynode_SUBTREE_TRANSFER(struct ynode_t *a, struct ynode_t *b)
+{
+ struct ynode_t *par;
+ struct ynode_t *del;
+ struct ynode_t *sib;
+ struct ynode_t *new;
+
+ if (a != NULL && b != NULL) {
+ if (ynode_is_equal(a, b)) {
+ /*printf("Identical subtrees. Transfer declined\n");*/
+ return;
+ }
+
+ /* The root node cannot participate */
+ if (ynode_is_root(a) || ynode_is_root(b)) {
+ /*printf("Cannot transfer root. Transfer declined\n");*/
+ return;
+ } else {
+
+ if (ynode_is_sibling(a, b)) {
+ /*printf("Sibling nodes: %d and %d.\n", a->id, b->id);*/
+ if (ynode_is_left_child(a)) {
+ a->P->L = b;
+ a->P->R = a;
+ } else {
+ a->P->L = a;
+ a->P->R = b;
+ }
+ return;
+ }
+
+ if (!ynode_is_disjoint(a, b)) {
+ /*printf("Non-disjoint subtrees. Interchange declined\n");*/
+ return;
+ }
+
+ new = ynode_add_before(b, '.', '.');
+ sib = ynode_get_sibling(a);
+ par = a->P;
+ a->P = new;
+
+ if (new->L == NULL) {
+ new->L = a;
+ } else {
+ new->R = a;
+ }
+
+ /* Detach a from the parent */
+ if (par->L && par->L->id == a->id) {
+ par->L = NULL;
+ sib = par->R;
+ } else {
+ par->R = NULL;
+ sib = par->L;
+ }
+
+ if (ynode_is_root(par)) {
+ /*
+ * Preserve the root, but prevent an extra
+ * internal node sneaking in when the root
+ * is the parent of a, by copying the
+ * remaining sibling to the root, and then
+ * removing it.
+ */
+ if (ynode_is_internal(sib)) {
+ par->L = sib->L;
+ par->R = sib->R;
+ par->L->P = par;
+ par->R->P = par;
+ ynode_destroy(sib);
+ }
+ } else {
+ /*
+ * Promote the sibling normally, if the
+ * parent node is not the root node.
+ */
+ del = ynode_promote(sib);
+ ynode_destroy(del);
+ }
+
+ }
+ }
+}
diff --git a/tree/node_print.c b/tree/node_print.c
new file mode 100644
index 0000000..f6121b7
--- /dev/null
+++ b/tree/node_print.c
@@ -0,0 +1,407 @@
+
+/*
+ * Based on an original from:
+ * http://www.openasthra.com/c-tidbits/printing-binary-trees-in-ascii/
+ *
+ * which appears to be down, although still in the Internet Archive.
+ */
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <limits.h>
+#include "ytree.h"
+
+
+struct asciinode_t {
+ struct asciinode_t *L;
+ struct asciinode_t *R;
+
+ /*
+ * Length of the edge to draw from
+ * this node to its children
+ */
+ int edge_length;
+
+ /* Height of the node (its rank) */
+ int height;
+
+ /* Length of the label string */
+ int label_length;
+
+ /*
+ * -1: I am the left child
+ * 0: I am the root
+ * 1: I am the right child
+ */
+ int parent_dir;
+
+ /*
+ * Label value. Max supported is uint32_t in
+ * decimal, 10 digits maximum.
+ */
+ char label[11];
+};
+
+
+#define MAX_HEIGHT 10000
+#define MY_INFINITY (1<<20)
+#define MIN(x, y) ((x) <= (y)) ? x : y
+#define MAX(x, y) ((x) <= (y)) ? y : x
+
+
+int lprofile[MAX_HEIGHT];
+int rprofile[MAX_HEIGHT];
+
+
+/* Adjust the gap between the left and right nodes */
+int gap = 3;
+
+/*
+ * X coordinate of the next character printed; used
+ * for printing the next node in the same rank.
+ */
+int print_next;
+
+
+
+
+
+/******************************************************************************
+ * CONVERT a normal binary tree to the ASCIItree
+ ******************************************************************************/
+
+/**
+ * build_ascii_tree_recursive()
+ * ----------------------------
+ * Recursively build out the ASCII tree.
+ *
+ * @n : Native binary tree node.
+ * @fmt : Format specifier for printing each value.
+ * Return: ASCII tree
+ */
+struct asciinode_t *build_ascii_tree_recursive(struct ynode_t *n, const char *fmt)
+{
+ struct asciinode_t *node;
+
+ if (n == NULL) {
+ return NULL;
+ }
+
+ node = malloc(sizeof(struct asciinode_t));
+ node->L = build_ascii_tree_recursive(n->L, fmt);
+ node->R = build_ascii_tree_recursive(n->R, fmt);
+
+ if (node->L != NULL) {
+ node->L->parent_dir = -1;
+ }
+
+ if (node->R != NULL) {
+ node->R->parent_dir = 1;
+ }
+
+ if (n->value == '.') {
+ sprintf(node->label, ".");
+ } else {
+ sprintf(node->label, fmt, n->key);
+ }
+
+ node->label_length = strlen(node->label);
+
+ return node;
+}
+
+
+/**
+ * build_ascii_tree()
+ * ------------------
+ * Copy the tree into the ASCII tree structure.
+ *
+ * @n : Native binary tree node.
+ * @fmt : Format specifier for printing each value.
+ * Return: ASCII tree
+ */
+struct asciinode_t *build_ascii_tree(struct ynode_t *n, const char *fmt)
+{
+ struct asciinode_t *node;
+
+ if (n == NULL) {
+ return NULL;
+ }
+ node = build_ascii_tree_recursive(n, fmt);
+ node->parent_dir = 0;
+ return node;
+}
+
+
+/**
+ * free_ascii_tree()
+ * -----------------
+ * Free the nodes of the ASCII tree.
+ *
+ * @node : ASCII tree
+ * Return: nothing
+ */
+void free_ascii_tree(struct asciinode_t *node)
+{
+ if (node == NULL) {
+ return;
+ }
+ free_ascii_tree(node->L);
+ free_ascii_tree(node->R);
+ free(node);
+}
+
+
+
+
+
+/******************************************************************************
+ * Utilities used internally by the ASCIItree
+ ******************************************************************************/
+
+
+/**
+ * print_level()
+ * -------------
+ * Print the given level of the ASCII tree
+ *
+ * @node : ASCII tree
+ * @x : x-coordinate of the node
+ * @level: level of the tree
+ * Return: nothing
+ */
+void print_level(struct asciinode_t *node, int x, int level)
+{
+ int i;
+ int is_left;
+
+ if (node == NULL) {
+ return;
+ }
+
+ is_left = (node->parent_dir == -1);
+
+ if (level == 0) {
+ for (i=0; i<(x-print_next-((node->label_length-is_left)/2)); i++) {
+ printf(" ");
+ }
+ print_next += i;
+ printf("%s", node->label);
+ print_next += node->label_length;
+ } else if (node->edge_length >= level) {
+ if (node->L != NULL) {
+ for (i=0; i<(x-print_next-(level)); i++) {
+ printf(" ");
+ }
+ print_next += i;
+ printf("/");
+ print_next += 1;
+ }
+ if (node->R != NULL) {
+ for (i=0; i<(x-print_next+(level)); i++) {
+ printf(" ");
+ }
+ print_next += i;
+ printf("\\");
+ print_next += 1;
+ }
+ } else {
+ print_level(node->L, x-node->edge_length-1, level-node->edge_length-1);
+ print_level(node->R, x+node->edge_length+1, level-node->edge_length-1);
+ }
+}
+
+
+
+/**
+ * compute_lprofile()
+ * ------------------
+ * Construct the lprofile array for the tree.
+ *
+ * NOTE
+ * Assumes the center of the root label is at (x,y).
+ * Assumes that edge_length fields have been computed.
+ *
+ * @node : ASCII tree
+ * @x : x-coordinate of the root node label
+ * @y : x-coordinate of the root node label
+ * Return: nothing
+ */
+void compute_lprofile(struct asciinode_t *node, int x, int y)
+{
+ int is_left;
+ int i;
+
+ if (node == NULL) {
+ return;
+ }
+
+ is_left = (node->parent_dir == -1);
+
+ lprofile[y] = MIN(lprofile[y], x-((node->label_length-is_left)/2));
+
+ if (node->L != NULL) {
+ for (i=1; (i<=node->edge_length) && ((y+i)<MAX_HEIGHT); i++) {
+ lprofile[y+i] = MIN(lprofile[y+i], x-i);
+ }
+ }
+ compute_lprofile(node->L, x-node->edge_length-1, y+node->edge_length+1);
+ compute_lprofile(node->R, x+node->edge_length+1, y+node->edge_length+1);
+}
+
+
+/**
+ * compute_rprofile()
+ * ------------------
+ * Construct the rprofile array for the tree.
+ *
+ * NOTE
+ * Assumes the center of the root label is at (x,y).
+ * Assumes that edge_length fields have been computed.
+ *
+ * @node : ASCII tree
+ * @x : x-coordinate of the root node label
+ * @y : x-coordinate of the root node label
+ * Return: nothing
+ */
+void compute_rprofile(struct asciinode_t *node, int x, int y)
+{
+ int not_left;
+ int i;
+
+ if (node == NULL) {
+ return;
+ }
+
+ not_left = (node->parent_dir != -1);
+
+ rprofile[y] = MAX(rprofile[y], x+((node->label_length-not_left)/2));
+
+ if (node->R != NULL) {
+ for (i=1; (i<=node->edge_length) && ((y+i)<MAX_HEIGHT); i++) {
+ rprofile[y+i] = MAX(rprofile[y+i], x+i);
+ }
+ }
+
+ compute_rprofile(node->L, x-node->edge_length-1, y+node->edge_length+1);
+ compute_rprofile(node->R, x+node->edge_length+1, y+node->edge_length+1);
+}
+
+
+
+/**
+ * compute_edge_lengths()
+ * ----------------------
+ * Construct the edge lengths for the tree.
+ *
+ * @node : ASCII tree
+ * Return: nothing
+ */
+void compute_edge_lengths(struct asciinode_t *node)
+{
+ int h;
+ int hmin;
+ int i;
+ int delta;
+
+ if (node == NULL) {
+ return;
+ }
+
+ compute_edge_lengths(node->L);
+ compute_edge_lengths(node->R);
+
+ /* First, fill in the edge_length of node */
+ if (node->R == NULL && node->L == NULL) {
+ node->edge_length = 0;
+ } else {
+ if (node->L != NULL) {
+ for (i=0; i<node->L->height && i < MAX_HEIGHT; i++) {
+ rprofile[i] = -MY_INFINITY;
+ }
+ compute_rprofile(node->L, 0, 0);
+ hmin = node->L->height;
+ } else {
+ hmin = 0;
+ }
+ if (node->R != NULL) {
+ for (i=0; i<node->R->height && i < MAX_HEIGHT; i++) {
+ lprofile[i] = MY_INFINITY;
+ }
+ compute_lprofile(node->R, 0, 0);
+ hmin = MIN(node->R->height, hmin);
+ } else {
+ hmin = 0;
+ }
+ delta = 4;
+ for (i=0; i<hmin; i++) {
+ delta = MAX(delta, gap + 1 + rprofile[i] - lprofile[i]);
+ }
+
+ /*
+ * If the node has 2 children of height 1,
+ * then we allow the 2 leaves to be within
+ * 1, instead of 2.
+ */
+ if (((node->L != NULL && node->L->height == 1)
+ || (node->R != NULL && node->R->height == 1))&&delta>4) {
+ delta--;
+ }
+
+ node->edge_length = ((delta+1)/2) - 1;
+ }
+
+ /* Now fill in the height of node */
+ h = 1;
+ if (node->L != NULL) {
+ h = MAX(node->L->height + node->edge_length + 1, h);
+ }
+ if (node->R != NULL) {
+ h = MAX(node->R->height + node->edge_length + 1, h);
+ }
+ node->height = h;
+}
+
+
+/******************************************************************************
+ * EXTERNAL FUNCTIONS
+ ******************************************************************************/
+
+
+void ynode_print(struct ynode_t *n, const char *fmt)
+{
+ struct asciinode_t *proot;
+ int xmin;
+ int i;
+
+ if (n == NULL) {
+ return;
+ }
+
+ proot = build_ascii_tree(n, fmt);
+
+ compute_edge_lengths(proot);
+
+ for (i=0; i<proot->height && i<MAX_HEIGHT; i++) {
+ lprofile[i] = MY_INFINITY;
+ }
+
+ compute_lprofile(proot, 0, 0);
+ xmin = 0;
+
+ for (i=0; i<proot->height && i<MAX_HEIGHT; i++) {
+ xmin = MIN(xmin, lprofile[i]);
+ }
+
+ for (i=0; i<proot->height; i++) {
+ print_next = 0;
+ print_level(proot, -xmin, i);
+ printf("\n");
+ }
+
+ if (proot->height >= MAX_HEIGHT) {
+ printf("(This tree is taller than %d, and may be drawn incorrectly.)\n", MAX_HEIGHT);
+ }
+
+ free_ascii_tree(proot);
+}
diff --git a/tree/node_traverse.c b/tree/node_traverse.c
new file mode 100644
index 0000000..6c1f146
--- /dev/null
+++ b/tree/node_traverse.c
@@ -0,0 +1,90 @@
+#include "ytree.h"
+
+/******************************************************************************
+ * TRAVERSAL
+ ******************************************************************************/
+
+static int Traversal_index = 1;
+
+void __impl__ynode_traverse_inorder(struct ynode_t *n, ynode_traverse_cb visit)
+{
+ if (n != NULL) {
+ __impl__ynode_traverse_inorder(n->L, visit);
+ visit(n, Traversal_index++);
+ __impl__ynode_traverse_inorder(n->R, visit);
+ }
+ return;
+}
+
+void __impl__ynode_traverse_preorder(struct ynode_t *n, ynode_traverse_cb visit)
+{
+ if (n != NULL) {
+ visit(n, Traversal_index++);
+ __impl__ynode_traverse_preorder(n->L, visit);
+ __impl__ynode_traverse_preorder(n->R, visit);
+ }
+ return;
+}
+
+void __impl__ynode_traverse_postorder(struct ynode_t *n, ynode_traverse_cb visit)
+{
+ if (n != NULL) {
+ visit(n, Traversal_index++);
+ __impl__ynode_traverse_postorder(n->L, visit);
+ __impl__ynode_traverse_postorder(n->R, visit);
+ }
+ return;
+}
+
+/**
+ * ynode_traverse_inorder()
+ * ------------------------
+ * Perform inorder traversal, applying a callback at each node.
+ *
+ * @n : Node to begin traversal from
+ * @visit: Callback applied at each node in the traversal.
+ * Return: Nothing
+ */
+void ynode_traverse_inorder(struct ynode_t *n, ynode_traverse_cb visit)
+{
+ Traversal_index = 1; /* prevent division by 0 */
+ if (n != NULL) {
+ __impl__ynode_traverse_inorder(n, visit);
+ }
+}
+
+/**
+ * ynode_traverse_preorder()
+ * -------------------------
+ * Perform preorder traversal, applying a callback at each node.
+ *
+ * @n : Node to begin traversal from
+ * @visit: Callback applied at each node in the traversal.
+ * Return: Nothing
+ */
+void ynode_traverse_preorder(struct ynode_t *n, ynode_traverse_cb visit)
+{
+ Traversal_index = 1;
+ if (n != NULL) {
+ __impl__ynode_traverse_preorder(n, visit);
+ }
+ return;
+}
+
+/**
+ * ynode_traverse_postorder()
+ * --------------------------
+ * Perform postorder traversal, applying a callback at each node.
+ *
+ * @n : Node to begin traversal from
+ * @visit: Callback applied at each node in the traversal.
+ * Return: Nothing
+ */
+void ynode_traverse_postorder(struct ynode_t *n, ynode_traverse_cb visit)
+{
+ Traversal_index = 1;
+ if (n != NULL) {
+ __impl__ynode_traverse_postorder(n, visit);
+ }
+ return;
+}
diff --git a/tree/old/pnode.c b/tree/old/pnode.c
new file mode 100644
index 0000000..506eb14
--- /dev/null
+++ b/tree/old/pnode.c
@@ -0,0 +1,31 @@
+#include <stdlib.h>
+#include <float.h>
+#include <stdio.h>
+#include <string.h>
+#include <math.h>
+#include <stdbool.h>
+#include "../prng/mersenne.h"
+#include "../prng/coin.h"
+#include "../prng/dice.h"
+#include "../util/math.h"
+#include "ptree.h"
+
+extern int DATA_COUNT;
+static int ID;
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/tree/old/pnode.h b/tree/old/pnode.h
new file mode 100644
index 0000000..d028b07
--- /dev/null
+++ b/tree/old/pnode.h
@@ -0,0 +1,137 @@
+#ifndef _OCC_PTREE_YNODE_H
+#define _OCC_PTREE_YNODE_H
+
+#include <stdbool.h>
+#include <stdint.h>
+#include "../prng/alias.h"
+
+/******************************************************************************
+ * SWITCHES AND OPTIONS
+ ******************************************************************************/
+
+#define PTREE_PRINT_DEBUG
+
+/* Define this for random insert paths (you want it) */
+#define PTREE_RANDOM_INSERT
+
+/* Label of internal nodes (should be disjoint from input alphabet) */
+#define PTREE_INTERNAL_NODE_LABEL '.'
+
+/******************************************************************************
+ * DATA TYPES
+ ******************************************************************************/
+
+typedef int ynode_value_t; /* Node label value */
+typedef uint32_t ynode_ident_t; /* Node ident value */
+
+
+/* Represents a node in the tree. */
+struct ynode_t {
+ struct ynode_t *L;
+ struct ynode_t *R;
+ struct ynode_t *P;
+ ynode_ident_t id;
+ ynode_value_t value;
+ ynode_value_t key;
+};
+
+
+/* Function pointer used in traversal methods. */
+typedef void (*ynode_traverse_cb)(struct ynode_t *n, int i);
+
+/* Function pointer used in distance methods. */
+typedef void (*ynode_distance_cb)(struct ynode_t *a, struct ynode_t *b);
+
+
+/******************************************************************************
+ * PREDICATES
+ ******************************************************************************/
+
+int ynode_is_subtree_of (struct ynode_t *a, struct ynode_t *b);
+int ynode_is_disjoint (struct ynode_t *a, struct ynode_t *b);
+int ynode_is_sibling (struct ynode_t *a, struct ynode_t *b);
+int ynode_is_equal (struct ynode_t *a, struct ynode_t *b);
+int ynode_is_left_child (struct ynode_t *a);
+int ynode_is_right_child (struct ynode_t *a);
+int ynode_is_full (struct ynode_t *a);
+int ynode_is_internal (struct ynode_t *a);
+int ynode_is_leaf (struct ynode_t *a);
+int ynode_is_root (struct ynode_t *a);
+int ynode_is_ternary (struct ynode_t *a);
+
+/******************************************************************************
+ * COUNTING
+ ******************************************************************************/
+
+int ynode_count_leaves (struct ynode_t *a);
+int ynode_count_leaves_not (struct ynode_t *a, struct ynode_t *x);
+int ynode_count_internal (struct ynode_t *a);
+
+/******************************************************************************
+ * CREATION AND DELETION
+ ******************************************************************************/
+
+struct ynode_t *ynode_create (ynode_value_t value, ynode_value_t key);
+struct ynode_t *ynode_create_root (void);
+struct ynode_t *ynode_copy (struct ynode_t *a);
+void ynode_destroy (struct ynode_t *a);
+
+/******************************************************************************
+ * INSERTION
+ ******************************************************************************/
+
+struct ynode_t *ynode_add_left (struct ynode_t *a, ynode_value_t value, ynode_value_t key);
+struct ynode_t *ynode_add_right (struct ynode_t *a, ynode_value_t value, ynode_value_t key);
+struct ynode_t *ynode_add_left_level (struct ynode_t *a, ynode_value_t value, ynode_value_t key);
+struct ynode_t *ynode_add_right_level (struct ynode_t *a, ynode_value_t value, ynode_value_t key);
+struct ynode_t *ynode_add_before (struct ynode_t *a, ynode_value_t value, ynode_value_t key);
+
+struct ynode_t *ynode_insert (struct ynode_t *a, ynode_value_t value, ynode_value_t key);
+void ynode_insert_on_path (struct ynode_t *r, struct ynode_t *a, char *path);
+
+/******************************************************************************
+ * TRAVERSAL
+ ******************************************************************************/
+
+void ynode_traverse_inorder (struct ynode_t *a, ynode_traverse_cb visit);
+void ynode_traverse_preorder (struct ynode_t *a, ynode_traverse_cb visit);
+void ynode_traverse_postorder (struct ynode_t *a, ynode_traverse_cb visit);
+
+/******************************************************************************
+ * RANDOM SELECTION
+ ******************************************************************************/
+
+struct ynode_t *ynode_random (struct ynode_t *n);
+struct ynode_t *ynode_random_leaf (struct ynode_t *n);
+struct ynode_t *ynode_random_internal (struct ynode_t *n);
+
+/******************************************************************************
+ * VALUES
+ ******************************************************************************/
+
+struct ynode_t *ynode_get_sibling (struct ynode_t *n);
+struct ynode_t *ynode_get_root (struct ynode_t *n);
+char *ynode_get_path (struct ynode_t *n);
+ynode_value_t *ynode_get_values (struct ynode_t *n);
+ynode_value_t *ynode_get_values_not (struct ynode_t *n, struct ynode_t *x);
+
+/******************************************************************************
+ * NODE-LEVEL MUTATIONS
+ ******************************************************************************/
+struct ynode_t *ynode_contract (struct ynode_t *n);
+struct ynode_t *ynode_promote (struct ynode_t *n);
+
+void ynode_LEAF_INTERCHANGE (struct ynode_t *a, struct ynode_t *b);
+void ynode_SUBTREE_INTERCHANGE(struct ynode_t *a, struct ynode_t *b);
+void ynode_SUBTREE_TRANSFER (struct ynode_t *a, struct ynode_t *b);
+
+/******************************************************************************
+ * COST
+ ******************************************************************************/
+float ynode_get_cost (struct ynode_t *n, float **distance);
+float ynode_get_cost_max(struct ynode_t *a, float **d, int n);
+float ynode_get_cost_min(struct ynode_t *a, float **d, int n);
+float ynode_get_cost_scaled(float c, float M, float m);
+
+
+#endif
diff --git a/tree/old/print.c b/tree/old/print.c
new file mode 100644
index 0000000..e9736a5
--- /dev/null
+++ b/tree/old/print.c
@@ -0,0 +1,406 @@
+/*
+ * Based on an original from:
+ * http://www.openasthra.com/c-tidbits/printing-binary-trees-in-ascii/
+ *
+ * which appears to be down, although still in the Internet Archive.
+ */
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <limits.h>
+#include "print.h"
+
+
+struct asciinode_t {
+ struct asciinode_t *L;
+ struct asciinode_t *R;
+
+ /*
+ * Length of the edge to draw from
+ * this node to its children
+ */
+ int edge_length;
+
+ /* Height of the node (its rank) */
+ int height;
+
+ /* Length of the label string */
+ int label_length;
+
+ /*
+ * -1: I am the left child
+ * 0: I am the root
+ * 1: I am the right child
+ */
+ int parent_dir;
+
+ /*
+ * Label value. Max supported is uint32_t in
+ * decimal, 10 digits maximum.
+ */
+ char label[11];
+};
+
+
+#define MAX_HEIGHT 10000
+#define INFINITY (1<<20)
+#define MIN(x, y) ((x) <= (y)) ? x : y
+#define MAX(x, y) ((x) <= (y)) ? y : x
+
+
+int lprofile[MAX_HEIGHT];
+int rprofile[MAX_HEIGHT];
+
+
+/* Adjust the gap between the left and right nodes */
+int gap = 3;
+
+/*
+ * X coordinate of the next character printed; used
+ * for printing the next node in the same rank.
+ */
+int print_next;
+
+
+
+
+
+/******************************************************************************
+ * CONVERT a normal binary tree to the ASCIItree
+ ******************************************************************************/
+
+/**
+ * build_ascii_tree_recursive()
+ * ----------------------------
+ * Recursively build out the ASCII tree.
+ *
+ * @n : Native binary tree node.
+ * @fmt : Format specifier for printing each value.
+ * Return: ASCII tree
+ */
+struct asciinode_t *build_ascii_tree_recursive(struct pnode_t *n, const char *fmt)
+{
+ struct asciinode_t *node;
+
+ if (n == NULL) {
+ return NULL;
+ }
+
+ node = malloc(sizeof(struct asciinode_t));
+ node->L = build_ascii_tree_recursive(n->L, fmt);
+ node->R = build_ascii_tree_recursive(n->R, fmt);
+
+ if (node->L != NULL) {
+ node->L->parent_dir = -1;
+ }
+
+ if (node->R != NULL) {
+ node->R->parent_dir = 1;
+ }
+
+ if (n->value == '.') {
+ sprintf(node->label, ".");
+ } else {
+ sprintf(node->label, fmt, n->key);
+ }
+
+ node->label_length = strlen(node->label);
+
+ return node;
+}
+
+
+/**
+ * build_ascii_tree()
+ * ------------------
+ * Copy the tree into the ASCII tree structure.
+ *
+ * @n : Native binary tree node.
+ * @fmt : Format specifier for printing each value.
+ * Return: ASCII tree
+ */
+struct asciinode_t *build_ascii_tree(struct pnode_t *n, const char *fmt)
+{
+ struct asciinode_t *node;
+
+ if (n == NULL) {
+ return NULL;
+ }
+ node = build_ascii_tree_recursive(n, fmt);
+ node->parent_dir = 0;
+ return node;
+}
+
+
+/**
+ * free_ascii_tree()
+ * -----------------
+ * Free the nodes of the ASCII tree.
+ *
+ * @node : ASCII tree
+ * Return: nothing
+ */
+void free_ascii_tree(struct asciinode_t *node)
+{
+ if (node == NULL) {
+ return;
+ }
+ free_ascii_tree(node->L);
+ free_ascii_tree(node->R);
+ free(node);
+}
+
+
+
+
+
+/******************************************************************************
+ * Utilities used internally by the ASCIItree
+ ******************************************************************************/
+
+
+/**
+ * print_level()
+ * -------------
+ * Print the given level of the ASCII tree
+ *
+ * @node : ASCII tree
+ * @x : x-coordinate of the node
+ * @level: level of the tree
+ * Return: nothing
+ */
+void print_level(struct asciinode_t *node, int x, int level)
+{
+ int i;
+ int is_left;
+
+ if (node == NULL) {
+ return;
+ }
+
+ is_left = (node->parent_dir == -1);
+
+ if (level == 0) {
+ for (i=0; i<(x-print_next-((node->label_length-is_left)/2)); i++) {
+ printf(" ");
+ }
+ print_next += i;
+ printf("%s", node->label);
+ print_next += node->label_length;
+ } else if (node->edge_length >= level) {
+ if (node->L != NULL) {
+ for (i=0; i<(x-print_next-(level)); i++) {
+ printf(" ");
+ }
+ print_next += i;
+ printf("/");
+ print_next += 1;
+ }
+ if (node->R != NULL) {
+ for (i=0; i<(x-print_next+(level)); i++) {
+ printf(" ");
+ }
+ print_next += i;
+ printf("\\");
+ print_next += 1;
+ }
+ } else {
+ print_level(node->L, x-node->edge_length-1, level-node->edge_length-1);
+ print_level(node->R, x+node->edge_length+1, level-node->edge_length-1);
+ }
+}
+
+
+
+/**
+ * compute_lprofile()
+ * ------------------
+ * Construct the lprofile array for the tree.
+ *
+ * NOTE
+ * Assumes the center of the root label is at (x,y).
+ * Assumes that edge_length fields have been computed.
+ *
+ * @node : ASCII tree
+ * @x : x-coordinate of the root node label
+ * @y : x-coordinate of the root node label
+ * Return: nothing
+ */
+void compute_lprofile(struct asciinode_t *node, int x, int y)
+{
+ int is_left;
+ int i;
+
+ if (node == NULL) {
+ return;
+ }
+
+ is_left = (node->parent_dir == -1);
+
+ lprofile[y] = MIN(lprofile[y], x-((node->label_length-is_left)/2));
+
+ if (node->L != NULL) {
+ for (i=1; (i<=node->edge_length) && ((y+i)<MAX_HEIGHT); i++) {
+ lprofile[y+i] = MIN(lprofile[y+i], x-i);
+ }
+ }
+ compute_lprofile(node->L, x-node->edge_length-1, y+node->edge_length+1);
+ compute_lprofile(node->R, x+node->edge_length+1, y+node->edge_length+1);
+}
+
+
+/**
+ * compute_rprofile()
+ * ------------------
+ * Construct the rprofile array for the tree.
+ *
+ * NOTE
+ * Assumes the center of the root label is at (x,y).
+ * Assumes that edge_length fields have been computed.
+ *
+ * @node : ASCII tree
+ * @x : x-coordinate of the root node label
+ * @y : x-coordinate of the root node label
+ * Return: nothing
+ */
+void compute_rprofile(struct asciinode_t *node, int x, int y)
+{
+ int not_left;
+ int i;
+
+ if (node == NULL) {
+ return;
+ }
+
+ not_left = (node->parent_dir != -1);
+
+ rprofile[y] = MAX(rprofile[y], x+((node->label_length-not_left)/2));
+
+ if (node->R != NULL) {
+ for (i=1; (i<=node->edge_length) && ((y+i)<MAX_HEIGHT); i++) {
+ rprofile[y+i] = MAX(rprofile[y+i], x+i);
+ }
+ }
+
+ compute_rprofile(node->L, x-node->edge_length-1, y+node->edge_length+1);
+ compute_rprofile(node->R, x+node->edge_length+1, y+node->edge_length+1);
+}
+
+
+
+/**
+ * compute_edge_lengths()
+ * ----------------------
+ * Construct the edge lengths for the tree.
+ *
+ * @node : ASCII tree
+ * Return: nothing
+ */
+void compute_edge_lengths(struct asciinode_t *node)
+{
+ int h;
+ int hmin;
+ int i;
+ int delta;
+
+ if (node == NULL) {
+ return;
+ }
+
+ compute_edge_lengths(node->L);
+ compute_edge_lengths(node->R);
+
+ /* First, fill in the edge_length of node */
+ if (node->R == NULL && node->L == NULL) {
+ node->edge_length = 0;
+ } else {
+ if (node->L != NULL) {
+ for (i=0; i<node->L->height && i < MAX_HEIGHT; i++) {
+ rprofile[i] = -INFINITY;
+ }
+ compute_rprofile(node->L, 0, 0);
+ hmin = node->L->height;
+ } else {
+ hmin = 0;
+ }
+ if (node->R != NULL) {
+ for (i=0; i<node->R->height && i < MAX_HEIGHT; i++) {
+ lprofile[i] = INFINITY;
+ }
+ compute_lprofile(node->R, 0, 0);
+ hmin = MIN(node->R->height, hmin);
+ } else {
+ hmin = 0;
+ }
+ delta = 4;
+ for (i=0; i<hmin; i++) {
+ delta = MAX(delta, gap + 1 + rprofile[i] - lprofile[i]);
+ }
+
+ /*
+ * If the node has 2 children of height 1,
+ * then we allow the 2 leaves to be within
+ * 1, instead of 2.
+ */
+ if (((node->L != NULL && node->L->height == 1)
+ || (node->R != NULL && node->R->height == 1))&&delta>4) {
+ delta--;
+ }
+
+ node->edge_length = ((delta+1)/2) - 1;
+ }
+
+ /* Now fill in the height of node */
+ h = 1;
+ if (node->L != NULL) {
+ h = MAX(node->L->height + node->edge_length + 1, h);
+ }
+ if (node->R != NULL) {
+ h = MAX(node->R->height + node->edge_length + 1, h);
+ }
+ node->height = h;
+}
+
+
+/******************************************************************************
+ * EXTERNAL FUNCTIONS
+ ******************************************************************************/
+
+
+void pnode_print(struct pnode_t *n, const char *fmt)
+{
+ struct asciinode_t *proot;
+ int xmin;
+ int i;
+
+ if (n == NULL) {
+ return;
+ }
+
+ proot = build_ascii_tree(n, fmt);
+
+ compute_edge_lengths(proot);
+
+ for (i=0; i<proot->height && i<MAX_HEIGHT; i++) {
+ lprofile[i] = INFINITY;
+ }
+
+ compute_lprofile(proot, 0, 0);
+ xmin = 0;
+
+ for (i=0; i<proot->height && i<MAX_HEIGHT; i++) {
+ xmin = MIN(xmin, lprofile[i]);
+ }
+
+ for (i=0; i<proot->height; i++) {
+ print_next = 0;
+ print_level(proot, -xmin, i);
+ printf("\n");
+ }
+
+ if (proot->height >= MAX_HEIGHT) {
+ printf("(This tree is taller than %d, and may be drawn incorrectly.)\n", MAX_HEIGHT);
+ }
+
+ free_ascii_tree(proot);
+}
diff --git a/tree/old/print.h b/tree/old/print.h
new file mode 100644
index 0000000..3e70d78
--- /dev/null
+++ b/tree/old/print.h
@@ -0,0 +1,9 @@
+#ifndef _OCC_PTREE_PRINT_H
+#define _OCC_PTREE_PRINT_H
+
+#include "ptree.h"
+#include "pnode.h"
+
+void pnode_print(struct pnode_t *n, const char *fmt);
+
+#endif
diff --git a/tree/old/ptree.c b/tree/old/ptree.c
new file mode 100644
index 0000000..723901f
--- /dev/null
+++ b/tree/old/ptree.c
@@ -0,0 +1,17 @@
+#include <stdlib.h>
+#include <float.h>
+#include <stdio.h>
+#include <string.h>
+#include <math.h>
+#include <stdbool.h>
+#include "../prng/prng.h"
+#include "../prng/coin.h"
+#include "../prng/dice.h"
+#include "../util/math.h"
+#include "ptree.h"
+
+
+
+
+
+
diff --git a/tree/old/ptree.h b/tree/old/ptree.h
new file mode 100644
index 0000000..0860427
--- /dev/null
+++ b/tree/old/ptree.h
@@ -0,0 +1,73 @@
+#ifndef _OCK_PTREE_H
+#define _OCK_PTREE_H
+
+#include <stdbool.h>
+#include <stdint.h>
+#include "../prng/alias.h"
+#include "pnode.h"
+
+/******************************************************************************
+ * DATA TYPES
+ ******************************************************************************/
+
+/* Handle for an entire tree */
+struct ptree_t {
+ pnode_ident_t count;
+ float **data;
+ int num_leaves;
+ int num_internal;
+ float max_cost;
+ float min_cost;
+ struct pnode_t *root;
+ struct pnode_t **memory;
+};
+
+
+/******************************************************************************
+ * HOUSEKEEPING
+ ******************************************************************************/
+struct ptree_t *ptree_create (int n, float **data);
+void ptree_free (struct ptree_t *tree);
+struct ptree_t *ptree_copy (struct ptree_t *tree);
+
+/******************************************************************************
+ * QTC
+ ******************************************************************************/
+float ptree_cost (struct ptree_t *tree);
+float ptree_cost_scaled (struct ptree_t *tree);
+
+/******************************************************************************
+ * MUTATIONS
+ ******************************************************************************/
+int ptree_mutate (struct ptree_t *tree, struct alias_t *alias);
+int ptree_mutate_mmc (struct ptree_t *tree, struct alias_t *alias);
+struct ptree_t *ptree_mutate_mmc2 (struct ptree_t *tree, struct alias_t *alias, int *num_mutations);
+
+
+/******************************************************************************
+ * SPECIAL PROPERTY FUNCTIONS
+ ******************************************************************************/
+
+/**
+ * sufficient_k()
+ * ``````````````
+ * @n : Number of leaf nodes in the phylogenetic tree
+ * Return: Minimum k allowing complete mutation.
+ *
+ * NOTE
+ * For the class {T} of phylogenetic trees with n leaf nodes,
+ * this function returns the minimal k such that for an
+ * arbitrary tree t in {T}, there exists a k-mutation which
+ * can transform t into arbitrary t' also in {T}.
+ *
+ * That is, given at least k mutations, any tree in the class
+ * can be transformed into any other tree also in the class.
+ */
+static inline int sufficient_k(int n)
+{
+ return (5 * n) - 16;
+}
+
+
+#endif
+
diff --git a/tree/tree_alloc.c b/tree/tree_alloc.c
new file mode 100644
index 0000000..cfaf13d
--- /dev/null
+++ b/tree/tree_alloc.c
@@ -0,0 +1,157 @@
+#include "ytree.h"
+
+/******************************************************************************
+ * TREE CREATE
+ ******************************************************************************/
+
+/**
+ * ytree_create()
+ * --------------
+ * Create an entire tree from an NxN data matrix
+ *
+ * @n : Number of data points
+ * @data : @nx@n data matrix.
+ * Return: Pointer to a tree structure.
+ */
+struct ytree_t *ytree_create(int n, float **data)
+{
+ struct ytree_t *tree;
+ int i;
+
+ tree = calloc(1, sizeof(struct ytree_t));
+
+ tree->root = ynode_create_root();
+ tree->data = data;
+
+ for (i=0; i<n; i++) {
+ ynode_insert(tree->root, i, i);
+ if (!ynode_is_ternary(tree->root)) {
+ fprintf(stderr, "Malformed tree.\n");
+ exit(1);
+ }
+ }
+
+ tree->count = n;
+ tree->num_leaves = ynode_count_leaves(tree->root);
+ tree->num_internal = ynode_count_internal(tree->root);
+ tree->max_cost = ynode_get_cost_max(tree->root, tree->data, n);
+ tree->min_cost = ynode_get_cost_min(tree->root, tree->data, n);
+
+ return tree;
+}
+
+
+
+/******************************************************************************
+ * TREE FREE
+ ******************************************************************************/
+
+/*
+ * This could be done much more efficiently, so
+ * do that at some point.
+ */
+struct ynode_t **Mem;
+int Mem_index;
+int Mem_max;
+
+void __impl__ytree_free(struct ynode_t *n, int i)
+{
+ if (n != NULL && Mem_index < Mem_max) {
+ Mem[Mem_index++] = n;
+ }
+}
+
+/**
+ * ytree_free()
+ * ------------
+ * Remove a tree from memory.
+ *
+ * @tree : Pointer to tree structure.
+ * Return: Nothing.
+ */
+void ytree_free(struct ytree_t *tree)
+{
+ int i;
+
+ Mem_max = tree->count + (tree->count-2);
+
+ Mem = calloc(Mem_max, sizeof(struct ytree_t *));
+ Mem_index = 0;
+
+ ynode_traverse_preorder(tree->root, __impl__ytree_free);
+
+ for (i=0; i<Mem_index; i++) {
+ if (Mem[i] != NULL) {
+ ynode_destroy(Mem[i]);
+ }
+ }
+
+ if (Mem != NULL) {
+ free(Mem);
+ }
+ if (tree != NULL) {
+ free(tree);
+ }
+}
+
+
+/******************************************************************************
+ * TREE COPY
+ ******************************************************************************/
+
+struct ynode_t *This_root;
+
+void __impl__ytree_copy(struct ynode_t *a, int i)
+{
+ struct ynode_t *copy;
+
+ if (a != NULL && !ynode_is_root(a)) {
+ copy = ynode_copy(a);
+
+ char *path = ynode_get_path(a);
+
+ ynode_insert_on_path(This_root, copy, path);
+
+ if (path != NULL) {
+ free(path);
+ }
+ }
+
+ return;
+}
+
+
+/**
+ * ytree_copy()
+ * ------------
+ * Copy an entire tree structure.
+ *
+ * @tree : Pointer to tree being copied.
+ * Return: Pointer to new tree, identical to @tree.
+ */
+struct ytree_t *ytree_copy(struct ytree_t *tree)
+{
+ if (tree == NULL) {
+ fprintf(stderr, "Cannot copy NULL tree\n");
+ return NULL;
+ }
+
+ struct ytree_t *copy;
+
+ copy = calloc(1, sizeof(struct ytree_t));
+
+ copy->data = tree->data;
+ copy->count = tree->count;
+ copy->num_leaves = tree->num_leaves;
+ copy->num_internal = tree->num_internal;
+ copy->max_cost = tree->max_cost;
+ copy->min_cost = tree->min_cost;
+
+ copy->root = ynode_copy(tree->root);
+
+ This_root = copy->root;
+
+ ynode_traverse_preorder(tree->root, __impl__ytree_copy);
+
+ return copy;
+}
diff --git a/tree/tree_cost.c b/tree/tree_cost.c
new file mode 100644
index 0000000..2e22247
--- /dev/null
+++ b/tree/tree_cost.c
@@ -0,0 +1,49 @@
+#include "ytree.h"
+
+/******************************************************************************
+ * TREE COST
+ ******************************************************************************/
+
+float **Cost_data;
+float Tree_cost;
+
+void __impl__ytree_cost(struct ynode_t *n, int i)
+{
+ if (ynode_is_internal(n)) {
+ Tree_cost += ynode_get_cost(n, Cost_data);
+ }
+ return;
+}
+
+/**
+ * ytree_cost()
+ * ------------
+ * Compute the un-normalized tree cost.
+ *
+ * @tree : Pointer to a tree structure.
+ * Return: Cost C(T) of the tree at @tree.
+ */
+float ytree_cost(struct ytree_t *tree)
+{
+ Cost_data = tree->data;
+ Tree_cost = 0;
+
+ ynode_traverse_preorder(tree->root, __impl__ytree_cost);
+
+ return Tree_cost;
+}
+
+
+/**
+ * ytree_cost_scaled()
+ * -------------------
+ * Compute the normalized tree cost.
+ *
+ * @tree : Pointer to a tree structure.
+ * Return: Cost S(T) of the tree at @tree.
+ */
+float ytree_cost_scaled(struct ytree_t *tree)
+{
+ float Ct = ytree_cost(tree);
+ return (tree->max_cost - Ct) / (tree->max_cost - tree->min_cost);
+}
diff --git a/tree/tree_mutate.c b/tree/tree_mutate.c
new file mode 100644
index 0000000..56dc69c
--- /dev/null
+++ b/tree/tree_mutate.c
@@ -0,0 +1,267 @@
+#include "ytree.h"
+
+/******************************************************************************
+ * TREE MUTATE
+ ******************************************************************************/
+
+/**
+ * ytree_mutate()
+ * --------------
+ * Perform various mutations on a tree structure, preserving shape invariants.
+ *
+ * @tree : Pointer to a tree structure.
+ * Return: Number of mutations made
+ */
+int ytree_mutate(struct ytree_t *tree, struct alias_t *alias)
+{
+ struct ynode_t *a;
+ struct ynode_t *b;
+ int r;
+ int m;
+ int i;
+
+ m = alias_sample(alias)+1;
+
+ for (i=0; i<m; i++) {
+
+ r = dice_roll(3);
+
+ switch (r) {
+ case 0:
+ a = ynode_get_random_leaf(tree->root);
+ b = ynode_get_random_leaf(tree->root);
+ ynode_LEAF_INTERCHANGE(a, b);
+ break;
+ case 1:
+ a = ynode_get_random(tree->root);
+ b = ynode_get_random(tree->root);
+ ynode_SUBTREE_INTERCHANGE(a, b);
+ break;
+ case 2:
+ a = ynode_get_random(tree->root);
+ b = ynode_get_random(tree->root);
+ ynode_SUBTREE_TRANSFER(a, b);
+ break;
+ }
+
+ if (!ynode_is_ternary(tree->root)) {
+ fprintf(stderr, "Malformed tree.\n");
+ exit(1);
+ }
+
+ if (tree->num_leaves != ynode_count_leaves(tree->root)) {
+ fprintf(stderr, "Malformed tree.\n");
+ exit(1);
+ }
+ }
+
+ return m;
+}
+
+
+/**
+ * ytree_mutate_mmc()
+ * ------------------
+ * Perform various mutations on a tree structure, preserving shape invariants.
+ *
+ * @tree : Pointer to a tree structure.
+ * Return: Number of mutations made
+ */
+int ytree_mutate_mmc(struct ytree_t *tree, struct alias_t *alias)
+{
+ struct ynode_t *a;
+ struct ynode_t *b;
+ float cost;
+ float best;
+ int r;
+ int m;
+ int i;
+
+ m = alias_sample(alias)+1;
+
+ best = ytree_cost(tree);
+
+ for (i=0; i<m; i++) {
+
+ r = dice_roll(3);
+
+ switch (r) {
+ case 0:
+ a = ynode_get_random_leaf(tree->root);
+ b = ynode_get_random_leaf(tree->root);
+ ynode_LEAF_INTERCHANGE(a, b);
+
+ /* Check cost and undo if bad step */
+ cost = ytree_cost(tree);
+
+ if (prng_uniform_random() < min_float(2, 1.0, cost/best)) {
+ /*if (cost > best) {*/
+ best = cost;
+ } else {
+ /* Undo mutation */
+ ynode_LEAF_INTERCHANGE(b, a);
+ }
+ break;
+ case 1:
+ a = ynode_get_random(tree->root);
+ b = ynode_get_random(tree->root);
+ ynode_SUBTREE_INTERCHANGE(a, b);
+
+ /* Check cost and undo if bad step */
+ cost = ytree_cost(tree);
+ /*if (cost > best) {*/
+ if (prng_uniform_random() < min_float(2, 1.0, cost/best)) {
+ best = cost;
+ } else {
+ /* Undo mutation */
+ ynode_SUBTREE_INTERCHANGE(b, a);
+ }
+ break;
+ case 2:
+ a = ynode_get_random(tree->root);
+ b = ynode_get_random(tree->root);
+ ynode_SUBTREE_TRANSFER(a, b);
+
+ /* Check cost and undo if bad step */
+ cost = ytree_cost(tree);
+ /*if (cost > best) {*/
+ if (prng_uniform_random() < min_float(2, 1.0, cost/best)) {
+ best = cost;
+ } else {
+ /* Undo mutation */
+ ynode_SUBTREE_INTERCHANGE(b, a);
+ }
+ break;
+ }
+
+ if (!ynode_is_ternary(tree->root)) {
+ fprintf(stderr, "Malformed tree.\n");
+ exit(1);
+ }
+
+ if (tree->num_leaves != ynode_count_leaves(tree->root)) {
+ fprintf(stderr, "Malformed tree.\n");
+ exit(1);
+ }
+ }
+
+ return m;
+}
+
+
+
+/**
+ * ytree_mutate_mmc2()
+ * -------------------
+ * Perform various mutations on a tree structure, preserving shape invariants.
+ *
+ * @tree : Pointer to a tree structure.
+ * Return: Number of mutations made
+ */
+struct ytree_t *ytree_mutate_mmc2(struct ytree_t *tree, struct alias_t *alias, int *num_mutations)
+{
+ struct ynode_t *a;
+ struct ynode_t *b;
+ float cost;
+ float init;
+ int r;
+ int m;
+ int i;
+
+ struct ytree_t *test;
+
+ m = alias_sample(alias)+1;
+
+ if (num_mutations != NULL) {
+ *num_mutations = m;
+ }
+
+ init = ytree_cost(tree);
+ test = ytree_copy(tree);
+
+ for (i=0; i<m; i++) {
+
+ r = dice_roll(3);
+
+ switch (r) {
+ case 0:
+ a = ynode_get_random_leaf(test->root);
+ b = ynode_get_random_leaf(test->root);
+ ynode_LEAF_INTERCHANGE(a, b);
+
+ /* Check cost and undo if bad step */
+ /*cost = tt_tree_cost(test);*/
+
+ /*if (uniform_random() < min_float(2, 1.0, cost/running)) {*/
+ /*[>[>if (cost > best) {<]<]*/
+ /*[>best = cost;<]*/
+ /*running = cost;*/
+ /*} else {*/
+ /*[>[> Undo mutation <]<]*/
+ /*tt_LEAF_INTERCHANGE(b, a);*/
+ /*}*/
+ break;
+ case 1:
+ a = ynode_get_random(test->root);
+ b = ynode_get_random(test->root);
+ ynode_SUBTREE_INTERCHANGE(a, b);
+
+ /* Check cost and undo if bad step */
+ /*cost = tt_tree_cost(test);*/
+ /*[>[>if (cost > best) {<]<]*/
+ /*if (uniform_random() < min_float(2, 1.0, cost/running)) {*/
+ /*[>best = cost;<]*/
+ /*running = cost;*/
+ /*} else {*/
+ /*[>[> Undo mutation <]<]*/
+ /*tt_SUBTREE_INTERCHANGE(b, a);*/
+ /*}*/
+ break;
+ case 2:
+ a = ynode_get_random(test->root);
+ b = ynode_get_random(test->root);
+ ynode_SUBTREE_TRANSFER(a, b);
+
+ /* Check cost and undo if bad step */
+ /*cost = tt_tree_cost(test);*/
+ /*[>[>if (cost > best) {<]<]*/
+ /*if (uniform_random() < min_float(2, 1.0, cost/running)) {*/
+ /*[>best = cost;<]*/
+ /*running = cost;*/
+ /*} else {*/
+ /*[>[> Undo mutation <]<]*/
+ /*tt_SUBTREE_INTERCHANGE(b, a);*/
+ /*}*/
+ break;
+ }
+
+ if (!ynode_is_ternary(test->root)) {
+ fprintf(stderr, "Malformed tree.\n");
+ exit(1);
+ }
+
+ if (test->num_leaves != ynode_count_leaves(test->root)) {
+ fprintf(stderr, "Malformed tree.\n");
+ exit(1);
+ }
+ }
+
+ cost = ytree_cost(test);
+
+ /*printf("cost(%f)/init(%f):%f\n", cost, init, cost/init);*/
+
+ /*if (cost > init) {*/
+ /*ytree_free(tree);*/
+ /*return test;*/
+ if (prng_uniform_random() < 1.0 - (cost/init)) {
+ /*if (uniform_random() < 1.0-(cost/init)) {*/
+ ytree_free(tree);
+ return test;
+ } else {
+ ytree_free(test);
+ return tree;
+ }
+}
+
+
+
diff --git a/tree/ytree.h b/tree/ytree.h
new file mode 100644
index 0000000..d0f43e7
--- /dev/null
+++ b/tree/ytree.h
@@ -0,0 +1,188 @@
+#ifndef __YTREE_H
+#define __YTREE_H
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <stdbool.h>
+#include <stdint.h>
+#include <string.h>
+#include "../prng/prng.h"
+#include "../prng/coin.h"
+#include "../prng/dice.h"
+#include "../prng/alias.h"
+#include "../util/math.h"
+
+extern int DATA_COUNT;
+
+/******************************************************************************
+ * SWITCHES AND OPTIONS
+ ******************************************************************************/
+
+#define YTREE_PRINT_DEBUG
+
+/* Define this for random insert paths (you want it) */
+#define YTREE_RANDOM_INSERT
+
+/* Label of internal nodes (should be disjoint from input alphabet) */
+#define YTREE_INTERNAL_NODE_LABEL '.'
+
+/******************************************************************************
+ * DATA TYPES
+ ******************************************************************************/
+
+typedef int ynode_value_t; /* Node label value */
+typedef uint32_t ynode_ident_t; /* Node ident value */
+
+
+/* Represents a node in the tree. */
+struct ynode_t {
+ struct ynode_t *L;
+ struct ynode_t *R;
+ struct ynode_t *P;
+ ynode_ident_t id;
+ ynode_value_t value;
+ ynode_value_t key;
+};
+
+
+/* Handle for an entire tree */
+struct ytree_t {
+ ynode_ident_t count;
+ float **data;
+ int num_leaves;
+ int num_internal;
+ float max_cost;
+ float min_cost;
+ struct ynode_t *root;
+ struct ynode_t **memory;
+};
+
+
+/* Function pointer used in traversal methods. */
+typedef void (*ynode_traverse_cb)(struct ynode_t *n, int i);
+
+/* Function pointer used in distance methods. */
+typedef void (*ynode_distance_cb)(struct ynode_t *a, struct ynode_t *b);
+
+
+/* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
+/* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
+/* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
+/* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
+/* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
+/* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
+
+
+/******************************************************************************
+ * ALLOCATION node_alloc.c
+ ******************************************************************************/
+struct ynode_t *ynode_create (ynode_value_t value, ynode_value_t key);
+struct ynode_t *ynode_create_root (void);
+struct ynode_t *ynode_copy (struct ynode_t *a);
+void ynode_destroy (struct ynode_t *a);
+
+/******************************************************************************
+ * INSERTION node_insert.c
+ ******************************************************************************/
+struct ynode_t *ynode_add_left (struct ynode_t *a, ynode_value_t value, ynode_value_t key);
+struct ynode_t *ynode_add_right (struct ynode_t *a, ynode_value_t value, ynode_value_t key);
+struct ynode_t *ynode_add_left_level (struct ynode_t *a, ynode_value_t value, ynode_value_t key);
+struct ynode_t *ynode_add_right_level (struct ynode_t *a, ynode_value_t value, ynode_value_t key);
+struct ynode_t *ynode_add_before (struct ynode_t *a, ynode_value_t value, ynode_value_t key);
+struct ynode_t *ynode_insert (struct ynode_t *a, ynode_value_t value, ynode_value_t key);
+void ynode_insert_on_path (struct ynode_t *r, struct ynode_t *a, char *path);
+
+/******************************************************************************
+ * SELECTION node_get.c
+ ******************************************************************************/
+struct ynode_t *ynode_get_random (struct ynode_t *n);
+struct ynode_t *ynode_get_random_leaf (struct ynode_t *n);
+struct ynode_t *ynode_get_random_internal(struct ynode_t *n);
+struct ynode_t *ynode_get_sibling (struct ynode_t *n);
+struct ynode_t *ynode_get_root (struct ynode_t *n);
+char *ynode_get_path (struct ynode_t *n);
+ynode_value_t *ynode_get_values (struct ynode_t *n);
+ynode_value_t *ynode_get_values_not (struct ynode_t *n, struct ynode_t *x);
+
+/******************************************************************************
+ * MUTATIONS node_mutate.c
+ ******************************************************************************/
+struct ynode_t *ynode_contract (struct ynode_t *n);
+struct ynode_t *ynode_promote (struct ynode_t *n);
+void ynode_LEAF_INTERCHANGE (struct ynode_t *a, struct ynode_t *b);
+void ynode_SUBTREE_INTERCHANGE(struct ynode_t *a, struct ynode_t *b);
+void ynode_SUBTREE_TRANSFER (struct ynode_t *a, struct ynode_t *b);
+
+/******************************************************************************
+ * COST node_cost.c
+ ******************************************************************************/
+float ynode_get_cost (struct ynode_t *n, float **distance);
+float ynode_get_cost_max (struct ynode_t *a, float **d, int n);
+float ynode_get_cost_min (struct ynode_t *a, float **d, int n);
+float ynode_get_cost_scaled (float c, float M, float m);
+
+/******************************************************************************
+ * COUNTING node_count.c
+ ******************************************************************************/
+int ynode_count_leaves (struct ynode_t *a);
+int ynode_count_leaves_not (struct ynode_t *a, struct ynode_t *x);
+int ynode_count_internal (struct ynode_t *a);
+
+/******************************************************************************
+ * TRAVERSAL node_traverse.c
+ ******************************************************************************/
+void ynode_traverse_inorder (struct ynode_t *a, ynode_traverse_cb visit);
+void ynode_traverse_preorder (struct ynode_t *a, ynode_traverse_cb visit);
+void ynode_traverse_postorder (struct ynode_t *a, ynode_traverse_cb visit);
+
+/******************************************************************************
+ * PREDICATES node_check.c
+ ******************************************************************************/
+int ynode_is_subtree_of (struct ynode_t *a, struct ynode_t *b);
+int ynode_is_disjoint (struct ynode_t *a, struct ynode_t *b);
+int ynode_is_sibling (struct ynode_t *a, struct ynode_t *b);
+int ynode_is_equal (struct ynode_t *a, struct ynode_t *b);
+int ynode_is_left_child (struct ynode_t *a);
+int ynode_is_right_child (struct ynode_t *a);
+int ynode_is_full (struct ynode_t *a);
+int ynode_is_internal (struct ynode_t *a);
+int ynode_is_leaf (struct ynode_t *a);
+int ynode_is_root (struct ynode_t *a);
+int ynode_is_ternary (struct ynode_t *a);
+
+/******************************************************************************
+ * PRINT node_print.c
+ ******************************************************************************/
+void ynode_print (struct ynode_t *n, const char *fmt);
+
+
+/* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
+/* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
+/* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
+/* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
+/* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
+/* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
+
+
+/******************************************************************************
+ * TREE ALLOCATION
+ ******************************************************************************/
+struct ytree_t *ytree_create (int n, float **data);
+void ytree_free (struct ytree_t *tree);
+struct ytree_t *ytree_copy (struct ytree_t *tree);
+
+/******************************************************************************
+ * TREE QTC
+ ******************************************************************************/
+float ytree_cost (struct ytree_t *tree);
+float ytree_cost_scaled (struct ytree_t *tree);
+
+/******************************************************************************
+ * TREE MUTATIONS
+ ******************************************************************************/
+int ytree_mutate (struct ytree_t *tree, struct alias_t *alias);
+int ytree_mutate_mmc (struct ytree_t *tree, struct alias_t *alias);
+struct ytree_t *ytree_mutate_mmc2 (struct ytree_t *tree, struct alias_t *alias, int *num_mutations);
+
+
+#endif
diff --git a/util/list.c b/util/list.c
new file mode 100644
index 0000000..4e9a03d
--- /dev/null
+++ b/util/list.c
@@ -0,0 +1,297 @@
+/*
+ "... A pointer to a structure object, suitably converted,
+ points to its initial member. There may be unnamed padding
+ within a structure object, but not at its beginning".
+
+ WG14/N1124 Section 6.7.2.1.13
+ ANSI C Standard
+*/
+#include <stdlib.h>
+#include "list.h"
+
+
+/**
+ * list_entry()
+ * ````````````
+ * Convert a list_node back into the structure containing it.
+ *
+ * @n : the list_node
+ * @type : the type of the entry
+ * @member: the list_node member of the type
+ */
+#define list_entry(n, type, member) container_of(n, type, member)
+
+
+/**
+ * list_alloc()
+ * ````````````
+ * Dynamically allocate and initialize a list object.
+ *
+ * Return: A new list object.
+ */
+struct list_t *list_alloc(void)
+{
+ struct list_t *new;
+
+ new = calloc(1, sizeof(struct list_t));
+ list_init(new);
+
+ return new;
+}
+
+
+/**
+ * list_init()
+ * ```````````
+ * Initialize a dynamically-allocated list object.
+ *
+ * @h : Address of the list object to be initialized.
+ * Return: Nothing.
+ */
+void list_init(struct list_t *head)
+{
+ head->n.next = head->n.prev = &head->n;
+}
+
+
+
+/******************************************************************************
+ * LIST STATUS
+ ******************************************************************************/
+
+/**
+ * list_empty()
+ * ````````````
+ * Is a list empty?
+ *
+ * @head : pointer to the list head
+ * Return: true if list is empty, else false.
+ */
+bool list_empty(const struct list_t *head)
+{
+ return head->n.next == &head->n;
+}
+
+
+/**
+ * list_length()
+ * `````````````
+ * Compute the length of a given list.
+ *
+ * @head : pointer to the list head
+ * Return: length of the list
+ */
+int list_length(const struct list_t *head)
+{
+ int count;
+ void *tmp;
+
+ count = 0;
+
+ list_foreach(head, tmp) {
+ count++;
+ }
+
+ return count;
+}
+
+
+
+/******************************************************************************
+ * ADD ITEMS TO LIST
+ ******************************************************************************/
+
+/**
+ * list_add()
+ * ``````````
+ * Add a new node to the start of a linked list.
+ *
+ * @head : address of the list to be operated on
+ * @structure: address of the list_node to add to the list at @head.
+ *
+ * NOTE
+ * The list_node does not need to be initialized; it will be overwritten.
+ */
+void list_add(struct list_t *head, void *structure)
+{
+ struct list_node *node = (struct list_node *)structure;
+
+ node->next = head->n.next;
+ node->prev = &head->n;
+ head->n.next->prev = node;
+ head->n.next = node;
+}
+
+
+/**
+ * list_append()
+ * `````````````
+ * Add an entry at the end of a linked list.
+ *
+ * @list : address of the list to be operated on
+ * @structure: address of the list_node to add to the list at @list.
+ *
+ * NOTE
+ * The list_node does not need to be initialized; it will be overwritten.
+ */
+void list_append(struct list_t *head, void *structure)
+{
+ struct list_node *node = (struct list_node *)structure;
+
+ node->next = &head->n;
+ node->prev = head->n.prev;
+ head->n.prev->next = node;
+ head->n.prev = node;
+}
+
+
+
+/******************************************************************************
+ * REMOVE ITEMS FROM LIST
+ ******************************************************************************/
+
+/**
+ * list_remove()
+ * `````````````
+ * Remove an entry from the linked list it is associated with.
+ *
+ * @structure: address of the list_node to remove from the (unknown) list.
+ * Return : Nothing
+ *
+ * NOTE
+ * This leaves @structure in an undefined state; it can be added to
+ * another list, but not deleted again.
+ */
+void list_remove(void *structure)
+{
+ struct list_node *node = (struct list_node *)structure;
+
+ node->next->prev = node->prev;
+ node->prev->next = node->next;
+
+ return;
+}
+
+
+/**
+ * list_remove_from()
+ * ``````````````````
+ * Remove an entry from a specific linked list.
+ *
+ * @list : address of the list holding @structure
+ * @structure: address of the list_node to remove from @head
+ *
+ * NOTE
+ * This explicitly indicates which list a node is expected to be in,
+ * which is better documentation and can catch more bugs.
+ */
+void list_remove_from(struct list_t *list, void *structure)
+{
+ struct list_node *node = (struct list_node *)structure;
+
+ #ifdef CCAN_LIST_DEBUG
+ {
+ /* Thorough check: make sure it was in list! */
+ struct list_node *i;
+ for (i = list->n.next; i != node; i = i->next) {
+ assert(i != &list->n);
+ }
+ }
+ #endif
+
+ /* Quick test that catches a surprising number of bugs. */
+ assert(!list_empty(list));
+ list_remove(node);
+}
+
+
+
+/******************************************************************************
+ * REFERENCE ITEMS IN LIST
+ ******************************************************************************/
+
+/**
+ * list_head()
+ * ```````````
+ * Get the first entry in a list.
+ *
+ * @head : address of the list head
+ * Return : Address of first entry in list, NULL if list is empty.
+ */
+const void *list_head(struct list_t *list)
+{
+ if (list_empty(list)) {
+ return NULL;
+ }
+
+ return (const char *)list->n.next;
+}
+
+
+
+/**
+ * list_tail()
+ * ```````````
+ * Get the last entry in a list
+ *
+ * @head : address of the list head
+ * Return: list item at the end of the list, NULL if list is empty.
+ */
+const void *list_tail(struct list_t *list)
+{
+ if (list_empty(list)) {
+ return NULL;
+ }
+
+ return (const char *)list->n.prev;
+}
+
+
+
+/**
+ * list_rotate_forward()
+ * `````````````````````
+ * Append the current head to the tail, then return the new head.
+ *
+ * @list : address of the list
+ * Return: the new list head (after rotation)
+ */
+const void *list_rotate_forward(struct list_t *list)
+{
+ struct list_node *node;
+
+ /* Quick test that catches a surprising number of bugs. */
+ assert(!list_empty(list));
+
+ node = (struct list_node *)list_head(list);
+ list_remove_from(list, node);
+ list_append(list, node);
+
+ return list_head(list);
+}
+
+
+/**
+ * list_rotate_backward()
+ * ``````````````````````
+ * Append the current tail to the head, then return the new head.
+ *
+ * @list : the list_t
+ * Return: the new list head (after rotation)
+ */
+const void *list_rotate_backward(struct list_t *list)
+{
+ struct list_node *node;
+
+ /* Quick test that catches a surprising number of bugs. */
+ assert(!list_empty(list));
+
+ node = (struct list_node *)list_tail(list);
+ list_remove_from(list, node);
+ list_add(list, node);
+
+ return list_head(list);
+}
+
+
+
diff --git a/util/list.h b/util/list.h
new file mode 100644
index 0000000..ae23b8d
--- /dev/null
+++ b/util/list.h
@@ -0,0 +1,153 @@
+#ifndef _LIST_HEADER_H
+#define _LIST_HEADER_H
+#include <stdbool.h>
+#include <assert.h>
+
+
+/******************************************************************************
+ * DATA STRUCTURES
+ ******************************************************************************/
+
+/**
+ * struct list_node
+ * ````````````````
+ * The object that makes any struct into a linked list node.
+ *
+ * @next: next entry (self if empty)
+ * @prev: previous entry (self if empty)
+ */
+struct list_node {
+ struct list_node *next;
+ struct list_node *prev;
+};
+
+
+
+/**
+ * struct list_t
+ * `````````````
+ * The head of a doubly linked list.
+ */
+struct list_t {
+ struct list_node n;
+};
+
+
+
+/******************************************************************************
+ * INITIALIZATION
+ ******************************************************************************/
+
+#define LIST(name, members) \
+ name { \
+ struct list_node node; \
+ members \
+ }
+
+
+/* Initialize a static list object. */
+#define LIST_INIT(list) { { &list.n, &list.n } }
+
+/* Declare and initialize a statically-allocated list object. */
+#define LIST_HEAD(name) struct list_t name = LIST_INIT(name)
+
+
+
+/******************************************************************************
+ * FUNCTIONS
+ ******************************************************************************/
+
+struct list_t *list_alloc (void);
+void list_init (struct list_t *list);
+void list_add (struct list_t *list, void *structure);
+void list_append (struct list_t *list, void *structure);
+void list_remove (void *structure);
+void list_remove_from (struct list_t *list, void *structure);
+const void *list_head (struct list_t *list);
+const void *list_tail (struct list_t *list);
+const void *list_rotate_forward (struct list_t *list);
+const void *list_rotate_backward(struct list_t *list);
+
+bool list_empty(const struct list_t *list);
+int list_length(const struct list_t *list);
+
+
+
+static inline struct list_node *to_node(void *ptr)
+{
+ return (struct list_node *)ptr;
+}
+
+static inline void *from_node(struct list_node *ptr)
+{
+ return (void *)ptr;
+}
+
+
+
+/******************************************************************************
+ * CONTROL STRUCTURES
+ ******************************************************************************/
+
+/**
+ * list_foreach
+ * ````````````
+ * Iterate over a list.
+ *
+ * @h : the list_t (warning: evaluated multiple times!)
+ * @i : the structure containing the list_node
+ *
+ * USAGE
+ * Iterate @i over the entire list rooted at @h. The macro wraps
+ * a 'for' loop, so you can break and continue as normal.
+ */
+#define list_foreach(h, i) \
+ for (i = from_node((h)->n.next); \
+ to_node((void *)i) != &(h)->n; \
+ i = from_node(to_node((void *)i)->next))
+
+
+
+/**
+ * list_foreach_rev
+ * ````````````````
+ * Iterate backwards over a list.
+ *
+ * @h: the list_t
+ * @i: the structure containing the list_node
+ *
+ * USAGE
+ * Iterate @i over the entire list rooted at @h. The macro wraps
+ * a 'for' loop, so you can break and continue as normal.
+ */
+#define list_foreach_rev(h, i) \
+ for (i = from_node((h)->n.prev); \
+ to_node((void *)i) != &(h)->n; \
+ i = from_node(to_node((void *)i)->prev))
+
+
+/**
+ * list_foreach_safe
+ * `````````````````
+ * Iterate through a list, maybe during deletion
+ *
+ * @h : the list_t
+ * @i : the structure containing the list_node
+ * @nxt : the structure containing the list_node
+ *
+ * USAGE
+ * This is a convenient wrapper to iterate @i over the
+ * entire list. It's a for loop, so you can break and
+ * continue as normal.
+ *
+ * The extra variable @nxt is used to hold the next
+ * element, so you can delete @i from the list.
+ */
+#define list_foreach_safe(h, i, nxt) \
+ for (i = from_node((h)->n.next), nxt = from_node(to_node((void *)i)->next); \
+ to_node((void *)i) != &(h)->n; \
+ i = nxt, nxt = from_node(to_node((void *)i)->next))
+
+
+
+#endif
diff --git a/util/math.c b/util/math.c
new file mode 100644
index 0000000..3713e7e
--- /dev/null
+++ b/util/math.c
@@ -0,0 +1,204 @@
+#include <stdlib.h>
+#include "math.h"
+
+
+static unsigned long gcd_ui(unsigned long x, unsigned long y)
+{
+ unsigned long t;
+
+ if (y < x) {
+ t = x;
+ x = y;
+ y = t;
+ }
+
+ /*
+ * y1<-(x0%y0), x1<-y0
+ */
+ while (y > 0) {
+ t = y;
+ y = x%y;
+ x = t;
+ }
+
+ return x;
+}
+
+unsigned long binomial(unsigned long n, unsigned long k)
+{
+ unsigned long d = 1;
+ unsigned long g = 1;
+ unsigned long r = 1;
+
+ if (k == 0) {
+ return 1;
+ }
+ if (k == 1) {
+ return n;
+ }
+ if (k >= n) {
+ return (k == n);
+ }
+ if (k > n/2) {
+ k = n-k;
+ }
+
+ for (d=1; d<=k; d++) {
+ /* Possible overflow */
+ if (r >= ULONG_MAX/n) {
+ /* Reduced numerator */
+ unsigned long nr;
+ /* Reduced denominator */
+ unsigned long dr;
+
+ g = gcd_ui(n, d);
+ nr = n/g;
+ dr = d/g;
+
+ g = gcd_ui(r, dr);
+ r = r/g;
+ dr = dr/g;
+
+ if (r >= ULONG_MAX/nr) {
+ /* Unavoidable overflow */
+ return 0;
+ }
+
+ r *= nr;
+ r /= dr;
+ n--;
+ } else {
+ r *= n--;
+ r /= d;
+ }
+ }
+ return r;
+}
+
+
+
+
+float max_float(int n, ...)
+{
+ va_list arg_ptr;
+ float max;
+ float num;
+ int i;
+
+ /* Requires the last positional parameter (to get the address) */
+ va_start(arg_ptr, n);
+
+ /* Provisionally set max to the first argument. */
+ max = va_arg(arg_ptr, double);
+
+ /* Loop over the other arguments to find the true min. */
+ for (i=1; i<n; i++) {
+ /*
+ * Requires the type to cast the argument into.
+ *
+ * Also increments the argument pointer up to the
+ * next argument.
+ */
+ num = va_arg(arg_ptr, double);
+ max = (max < num) ? num : max;
+ }
+
+ /* Free memory and reset. */
+ va_end(arg_ptr);
+
+ return max;
+}
+
+
+float min_float(int n, ...)
+{
+ va_list arg_ptr;
+ float min;
+ float num;
+ int i;
+
+ /* Requires the last positional parameter (to get the address) */
+ va_start(arg_ptr, n);
+
+ /* Provisionally set min to the first argument. */
+ min = va_arg(arg_ptr, double);
+
+ /* Loop over the other arguments to find the true min. */
+ for (i=1; i<n; i++) {
+ /*
+ * Requires the type to cast the argument into.
+ *
+ * Also increments the argument pointer up to the
+ * next argument.
+ */
+ num = va_arg(arg_ptr, double);
+ min = (min > num) ? num : min;
+ }
+
+ /* Free memory and reset. */
+ va_end(arg_ptr);
+
+ return min;
+}
+
+
+int max_int(int n, ...)
+{
+ va_list arg_ptr;
+ int max;
+ int num;
+ int i;
+
+ /* Requires the last positional parameter (to get the address) */
+ va_start(arg_ptr, n);
+
+ max = 0;
+
+ for (i=0; i<n; i++) {
+ /*
+ * Requires the type to cast the argument into.
+ *
+ * Also increments the argument pointer up to the
+ * next argument.
+ */
+ num = va_arg(arg_ptr, int);
+ max = (max < num) ? num : max;
+ }
+
+ /* Free memory and reset. */
+ va_end(arg_ptr);
+
+ return max;
+}
+
+
+int min_int(int n, ...)
+{
+ va_list arg_ptr;
+ int min;
+ int num;
+ int i;
+
+ /* Requires the last positional parameter (to get the address) */
+ va_start(arg_ptr, n);
+
+ min = 0;
+
+ for (i=0; i<n; i++) {
+ /*
+ * Requires the type to cast the argument into.
+ *
+ * Also increments the argument pointer up to the
+ * next argument.
+ */
+ num = va_arg(arg_ptr, int);
+ min = (min > num) ? num : min;
+ }
+
+ /* Free memory and reset. */
+ va_end(arg_ptr);
+
+ return min;
+}
+
+
diff --git a/util/math.h b/util/math.h
new file mode 100644
index 0000000..2b35d2a
--- /dev/null
+++ b/util/math.h
@@ -0,0 +1,19 @@
+#ifndef __MATH_COMBINATORICS_H
+#define __MATH_COMBINATORICS_H
+
+#include <stdarg.h>
+#include <stdint.h>
+#include <limits.h>
+#include <math.h>
+
+
+unsigned long binomial(unsigned long n, unsigned long k);
+
+
+float max_float(int n, ...);
+float min_float(int n, ...);
+int max_int (int n, ...);
+int min_int (int n, ...);
+
+#endif
+