Exponential Generating Function of rooted minimal directed acyclic graphs I am trying to find the exponential generating function 
of directed minimal acyclic graphs (which I now call dag), 
where every non-leaf node has two outgoing edges.  
Context: A simple tree compression algorithm consists of saving repeated subtrees only once. 
Further occurrences of repeated trees are simply linked to the first occurrence. 
This way one gets a unique minimal directed acyclic graph, and since we started with a tree it's also rooted. 
For simplicity I would like to treat binary trees, hence two outgoing edges per non-leaf node.
A natural question is how big dags of binary trees of size $n$ are, 
and the question has been answered 
here 
(the paper is Analytic Variations on the Common Subexpression Problem, by Flajolet et al).  
I would like to ask a different question, namely how many different dags of size $n$ there are, or equivalently, how many rooted plane binary trees have a dag of size $n$?
As an example, for $n=3$, we have three such trees, namely 
$a(a(a,a),a(a,a))$, $a(a(a,a),a)$ and $a(a,a(a,a))$.
For $n=4$ there are $15$ trees, for $n=5$ there are $111$.
A promising sequence from OEIS is A001063, but I can neither make sense of the differential equation mentioned there, nor do I have a combinatorial explanation for the formula there that calculates $a_{n+1}$, given $a_1,\dots,a_n$:
$$
a_{n+1} = \sum_{k=0..n} \frac{n!}{k!} \cdot \binom{n-1}{k-1}\cdot a_k
$$
If requested, I could add where I got stuck (I mostly tried to make sense of the formula), but I think my post is already too long.
Addendum. This question has generated its own OEIS-series (A254789)! Thanks to everyone involved!
 A: OK, now it looks better: I'm quite confident the sequence starts
(1,1)
(2,3)
(3,15)
(4,111)
(5,1119)
(6,14487)
(7,230943)
(8,4395855)
(9,97608831)
(10,2482988079)
(11,71321533887)
(12,2286179073663)
(13,80984105660415)
(14,3144251526824991)
(15,132867034410319359)

and that's computed within a few seconds, using the following approach:
based on function count :: [[Bool]] -> Int where count xss is the number of dags with map length xss nodes at the respective level, and in each level, coded by an element xs :: [Bool] of xss, the entries of xs mark whether this node should have a predecessor.
In more detail, here's the specification of count:
We define a function (just for specification, it is not in the source below) shape :: DAG -> [[Bool]] that takes a DAG (any DAG, may have several roots), computes the list of level sets, then for each set, a canonical ordering (a list) of its nodes (lexicographic by left-child, right child, using the ordering in the lower levels), then for each node, whether it has a predecessor (a node higher up that points here).
Now count s gives the number of DAGs d that have shape d == s.
The point is that we can define count recursively (induction by the number of levels), and we never really construct DAGs - we just count.
And while we count, we avoid recomputations, using memoFix (a fixpoint combinator with a cache, really). You may simply think count arg = case arg ... return $ count ...
To run this with ghc, you need packages lens and memoize.
You can load the source code in ghci and evaluate expressions like count [[False],[True],[True]]. (It seems the code indentation here is broken. Watch out that expressions inside do are aligned properly.)
import Control.Monad ( guard, forM_ )
import Control.Applicative
import Control.Lens
import Data.List (tails, sort)
import Data.Function.Memoize
import System.IO

type Shape = [[Bool]]

main = forM_ [ 1 .. ] $ \ s -> do
       print ( s, sum $ map count $ shapes s )
   hFlush stdout

shapes s = do  sh <- deep_shapes (s-1) ;  return $ [False] : sh

deep_shapes :: Int -> [Shape]
deep_shapes 0 = return []
deep_shapes s = do
  x <- [ 1 .. s ] ; xs <- deep_shapes (s-x)
  return $ (replicate x True) : xs

count :: Shape -> Int
count = memoFix $ \ self arg -> case arg of
       [] -> 1
       (sh : ape) -> sum $ do
      guard $ and $ map not sh
      top <- pairs (length sh) ape
      return $ self $ apply top ape

type Node = (Int,Int)
type Pair = (Node,Node)

apply :: [Pair] -> Shape -> Shape
apply top shape = 
    foldr ( \ (h,k) sh -> sh & ix (length shape - h) . ix k .~ False ) 
      shape $ do (p,q) <- top ; [p,q]

pairs s shape = pick s $ sort $ do
    let cs = candidates shape
        lower = concat $ drop 1 cs
            top = concat $ take 1 cs
    (left,right) <- [(lower,top),(top,top),(top,lower)]
    (,) <$> left <*> right

candidates :: [[Bool]] -> [[(Int,Int)]]
candidates shape = ( do
   (h,ops) <- zip [length shape, length shape-1 ..] shape
   return $ do (n, _) <- zip [0..] ops ; return (h,n) ) ++ [[(0,0)]]

pick :: Int -> [a] -> [[a]]
pick 0 _ = return []
pick s xs = do
  z : ys <- tails xs ; guard $ length ys >= s-1
      zs <- pick (s-1) ys ; return $ z : zs

A: The  following C  program  uses  POSIX threads  to  create a  parallel
implementation of  the single-threaded C program I  posted earlier. It
can  be used to  compute the  distributions for  $n=8.$ This  had peak
memory allocation $2.1$GB and took $23$ minutes on a machine with $24$
processors at $2$GHz. It would appear that the distributions for $n=9$
could perhaps be  computed on a machine with  at least $32$ processors
and at least $16$GB of memory, which the reader is invited to try.

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <assert.h>
#include <time.h>
#include <unistd.h>
#include <pthread.h>

#define MXDST 12

typedef struct th {
  struct th *distsub[MXDST];
  int size;
} tree_inst, *tree_ptr;

tree_ptr tree_new(tree_ptr left, tree_ptr right, int mx)
{
  tree_ptr mergebuf[2*mx];
  int pos, idx;

  memcpy(mergebuf, left->distsub,
         left->size*sizeof(tree_ptr));

  pos = left->size;
  for(idx = 0; idx < right->size; idx++){
    tree_ptr item = right->distsub[idx];

    int cmp;
    for(cmp=0; cmp < left->size; cmp++){
      if(mergebuf[cmp] == item) break;
    }

    if(cmp == left->size){
      mergebuf[pos] = item;
      pos++;
    }
  }

  if(pos >= mx) return NULL;

  tree_ptr item = malloc(sizeof(tree_inst));
  assert(item != NULL);

  memcpy(item->distsub, mergebuf,
         pos*sizeof(tree_ptr));

  item->distsub[pos] = item;
  item->size = pos+1;

  return item;
}

typedef struct {
  tree_ptr *data;
  int alloc;
  int count;
  pthread_mutex_t mutex;
} coll_inst, *coll_ptr;

#define COLL_CHUNK 512

coll_ptr coll_new(void)
{
  coll_ptr item = malloc(sizeof(coll_inst));
  assert(item != NULL);

  item->data = malloc(COLL_CHUNK*sizeof(tree_ptr));
  assert(item->data != NULL);

  item->alloc = COLL_CHUNK;
  item->count = 0;

  item->mutex = 
    (pthread_mutex_t)PTHREAD_MUTEX_INITIALIZER;

  return item;
}

coll_ptr coll_record(coll_ptr coll, tree_ptr entry)
{
  if(coll->count == coll->alloc){
    coll->alloc += COLL_CHUNK;

    coll->data = realloc(coll->data, 
                        coll->alloc*sizeof(tree_ptr));
    assert(coll->data != NULL);
  }

  coll->data[coll->count++] = entry;
  return coll;
}

typedef struct {
  coll_ptr *table_ptr;
  int mx;
  int n, m;
} state_info;

void *compute(void *ptr)
{
  state_info *st = (state_info *)ptr;

  int mx = st->mx;
  int n = st->n, m = st->m;
  int dst1, dst2;

  for(dst1 = 0; dst1 < st->mx; dst1++){
    for(dst2 = 0; dst2 < st->mx; dst2++){
      coll_ptr 
        ca = st->table_ptr[m*(mx+1) + dst1],
        cb = st->table_ptr[(n-1-m)*(mx+1) + dst2];

      int c1, c2;
      for(c1 = 0; c1 < ca->count; c1++){
        for(c2 = 0; c2 < cb->count; c2++){
          tree_ptr 
            t1 = ca->data[c1],
            t2 = cb->data[c2];

          tree_ptr tree = tree_new(t1, t2, mx);
          if(tree != NULL){
            coll_ptr targ =
              st->table_ptr[n*(mx+1)+tree->size];

            pthread_mutex_lock(&(targ->mutex));
            coll_record(targ, tree);
            pthread_mutex_unlock(&(targ->mutex));
          }
        }
      }
    }
  }

  return NULL;
}

int main(int argc, char **argv)
{
  int mx = 1;

  if(argc>1){
    int mxcmd = atoi(argv[1]);

    if(1 <= mxcmd && mxcmd <= MXDST){
      mx = mxcmd;
    }
    else{
      fprintf(stderr, "invalid maxdist value, "
              "got %d\n", mxcmd);
      exit(-1);
    }
  }

  int nmx = 1 << mx;
  coll_ptr table[nmx+1][mx+1];

  int n, dst;
  for(n=0; n <= nmx; n++){
    for(dst = 0; dst <= mx; dst++){
      table[n][dst] = coll_new();
    }
  }

  int grand[mx+1];
  for(dst = 0; dst <= mx; dst++){
    grand[dst] = 0;
  }


  tree_inst base;
  base.distsub[0] = NULL;
  base.size = 0;
  coll_record(table[0][0], &base);

  for(n=1; n <= nmx; n++){
    int m; time_t time_begin, time_end;

    time_begin = time(NULL);

    pthread_t threads[n];

    state_info state[n];

    for(m=0; m <= n-1; m++){
      state[m].table_ptr = (coll_ptr *)table;
      state[m].mx = mx;
      state[m].n = n;
      state[m].m = m;

      pthread_create(threads+m, NULL, compute, state+m);
    }

    for(m=0; m <= n-1; m++){
      pthread_join(threads[m], NULL);
    }

    time_end = time(NULL);

    printf("%d: ", n);

    int total = 0, ents = 0;
    for(dst = 1; dst <= mx; dst++){
      int cval = table[n][dst]->count;

      if(cval > 0){
        if(ents > 0) printf(" + ");

        if(cval > 1) printf("%d ", cval);
        printf("u");
        if(dst > 1) printf("^%d", dst);

        total += cval; ents++;
        grand[dst] += cval;
      }
    }

    if(ents == 0) printf("0");
    printf(" (%d)\n", total);

    int secs = (int)difftime(time_end, time_begin);
    if(!isatty(fileno(stdout)))
      fprintf(stderr, "%d [%ds]\n", n, secs);
  }

  printf("-\n");
  for(dst = 1; dst <= mx; dst++){
    if(dst > 1) printf(" + ");

    if(grand[dst] > 1) printf("%d ", grand[dst]);
    printf("u");
    if(dst > 1) printf("^%d", dst);
  }

  printf("\n");

  exit(0);
}

