How many connected graphs over V vertices and E edges?

Is there a way to calculate the number of simple connected graphs possible over given edges and vertices?

Eg: 3 vertices and 2 edges will have 3 connected graphs But 3 vertices and 3 edges will have 1 connected graph

Then 4 edges and 3 will have 4 connected graphs

Till such values...it is easy to see its

V choose E

But what about when the number of vertices are less than number of edges...how to calculate then?

I am not able to visualize that

Can it be a variation of the Stars and Bars problem

Like...number of ways 7 edges(balls) can be connected to (kept in) 5 vertices(bags) such that no vertex(bag) is isolated (is empty)

Here maybe we might have to consider the number of edges twice as each edge needs two vertices...


Solution 1:

First of all, let me state my preconditions. Since you write that there are three graphs with two edges on three vertices it seems you are talking about the labelled case, which is what I will work with from now on. As this is truly a vast field of investigation I will just show you how to calculate these numbers (connected graphs on $n$ nodes having $k$ edges). This should enable you to consult the relevant entries of the OEIS and decide what course to take in your research.

Method I. Let $\mathcal{Q}$ be the combinatorial class of connected graphs and $\mathcal{G}$ the combinatorial class of labelled graphs, all of them. The relation between these two classes is the set-of relation: a graph is a set of connected components. This gives the relation between the two classes: $$\def\textsc#1{\dosc#1\csod} \def\dosc#1#2\csod{{\rm #1{\small #2}}}\mathcal{G} = \textsc{SET}(\mathcal{Q}).$$

Translating to generating functions this means that $$G(z) = \exp Q(z)$$ and hence $$Q(z) = \log G(z).$$

Now observe that the mixed exponential generating function of $\mathcal{G}$ by vertices and edge count is not difficult to calculate, it is simply $$G(z) = \sum_{m\ge 0} (1+u)^{m(m-1)/2} \frac{z^m}{m!} = 1 + \sum_{m\ge 1} (1+u)^{m(m-1)/2} \frac{z^m}{m!}.$$ This yields for $Q(z)$ that $$Q(z) = \log\left(1+ \sum_{m\ge 1} (1+u)^{m(m-1)/2} \frac{z^m}{m!}\right) = \sum_{q\ge 1} (-1)^{q+1} \frac{1}{q} \left(\sum_{m\ge 1} (1+u)^{m(m-1)/2} \frac{z^m}{m!}\right)^q.$$ We are interested in the count for $n$ vertices and $k$ edges where $k\ge n-1$. Note that the term in the parenthesis has minimum degree $q$ in $z$, so we can omit the rest of the series where $q>n.$ This finally yields the formula for the number $q_{n,k}$ of connected labelled graphs with $k$ edges: $$q_{n,k} = n! [z^n] [u^k] \sum_{q=1}^n (-1)^{q+1} \frac{1}{q} \left(\sum_{m=1}^n (1+u)^{m(m-1)/2} \frac{z^m}{m!}\right)^q.$$ Substituting concrete values into this formula and entering it into your favorite CAS yields for $k=n-1$ the sequence $$1, 1, 3, 16, 125, 1296, 16807, 262144, 4782969, 100000000,\ldots$$ which is OEIS A000272, the tree sequence with value $n^{n-2}.$ Setting $k=n$, we get $$0, 0, 1, 15, 222, 3660, 68295, 1436568, 33779340, 880107840,\ldots$$ which is OEIS A057500. Continuing with $k=n+1$ we have $$0, 0, 0, 6, 205, 5700, 156555, 4483360, 136368414, 4432075200, 154060613850,\ldots$$ which is OEIS A061540.

We could keep going like this but the pattern should be clear.

Method II. A different approach uses the exponential generating function of the set of rooted labelled trees given by $$T(z) = \sum_{m\ge 1} m^{m-1} \frac{z^m}{m!}.$$ The method is to use a combinatorial decomposition of the relevant classes of graphs in terms of a reduced structure consisting of cycles to which trees are attached at the nodes. Roughly speaking this structure is what you obtain by removing vertices of degree one from the graph until there are none left, and merging adjacent vertices of degree two. The result is a multigraph. The connected graphs that fall into the class of graphs corresponding to this multigraph are obtained by placing sequences of trees on the multigraph edges such that no self-loops or multi-edges remain.

Note the the so-called excess of a connected graph is the number of edges plus one minus the number of vertices. That means that a tree has excess zero. We will now compute the generating functions for graphs of excess zero, one, and two.

For starters, we have $$q_{n, n-1} = \frac{1}{n} \times n! [z^n] T(z)$$ which produces the correct sequence (the division accounts for the difference between rooted and unrooted trees).

Next we have that the graphs with excess one counted by $q_{n,n}$ consist of a cycle with trees attached. Summing over the size $m$ of the cycle this yields $$q_{n, n} = n! [z^n] \sum_{m=3}^n \frac{T(z)^m}{2m}$$ where the two in the denominator accounts for the fact that the cycle is not directed (dihedral group with $2m$ permutations).

Finally in calculating $q_{n,n+1}$ it becomes evident that the underlying structure consists of two cycles joined at a common node or by a path, or a cycle with a chord, which in fact turns out to be two vertices joined by three edges.

Start the inventory. If we have two cycles that are joined at a common node the resulting operator has four or eight automorphisms depending on whether the two cycles have the same size, for a contribution of $$T(z) \left(\sum_{m_1\ge 2} \sum_{m_2\ge m_1+1} \frac{T(z)^{m_1+m_2}}{4} + \sum_{m\ge 2} \frac{T(z)^{2m}}{8}\right).$$ If they are joined by a path we must place a sequence of trees on that path and the contribution is $$ \frac{T(z)^2}{1-T(z)}\left( \sum_{m_1\ge 2} \sum_{m_2\ge m_1+1} \frac{T(z)^{m_1+m_2}}{4} + \sum_{m\ge 2} \frac{T(z)^{2m}}{8}\right).$$ With the chord there are two cases -- one chord stays empty or all three chords are occupied. For the first case we get $$T(z)^2 \left( \sum_{m_1\ge 1} \sum_{m_2\ge m_1+1} \frac{T(z)^{m_1+m_2}}{2} + \sum_{m\ge 1} \frac{T(z)^{2m}}{4}\right).$$ The second case produces $$T(z)^2 \left( \sum_{m_1\ge 1}\sum_{m_2\ge m_1+1}\sum_{m_3\ge m_2+1} \frac{T(z)^{m_1+m_2+m_3}}{2} + \sum_{m_1\ge 1}\sum_{m_2\ge m_1+1} \frac{T(z)^{m_1+2m_2}}{4}\\ + \sum_{m_1\ge 1}\sum_{m_2\ge m_1+1} \frac{T(z)^{2m_1+m_2}}{4} + \sum_{m\ge 1} \frac{T(z)^{3m}}{12} \right).$$

Adding these four contributions yields the generating function $$\frac{T(z)^4}{24} \frac{6-T(z)}{(1-T(z))^3}$$ and the formula $$q_{n, n+1} = n! [z^n] \frac{T(z)^4}{24} \frac{6-T(z)}{(1-T(z))^3}.$$

As for the practical utility of this formula it can be treated by Lagrange inversion to give a closed form suitable for computation. The species of labelled trees has the specification $$\mathcal{T} = \mathcal{Z} \times \mathfrak{P}(\mathcal{T})$$ which gives the functional equation $$T(z) = z \exp T(z).$$

Extracting coefficients via Lagrange inversion we have $$q_{n,n+1} = n! \frac{1}{2\pi i} \int_{|z|=\epsilon} \frac{1}{z^{n+1}} \frac{T(z)^4}{24} \frac{6-T(z)}{(1-T(z))^3} dz.$$

Put $T(z)=w$ so that $z=w/\exp(w) = w\exp(-w)$ and $dz = (\exp(-w) - w\exp(-w)) \; dw$ to get $$n! \frac{1}{2\pi i} \int_{|w|=\epsilon} \frac{\exp(w(n+1))}{w^{n+1}} \times \frac{w^4}{24} \frac{6-w}{(1-w)^3} \times (\exp(-w) - w\exp(-w)) dw \\ = \frac{1}{24} n! \frac{1}{2\pi i} \int_{|w|=\epsilon} \frac{\exp(wn)}{w^{n-3}} \frac{6-w}{(1-w)^2} dw \\ = \frac{1}{24} n! \frac{1}{2\pi i} \int_{|w|=\epsilon} \frac{\exp(wn)}{w^{n-3}} \left(\frac{1}{1-w} + 5\frac{1}{(1-w)^2}\right) dw.$$

Extracting coefficients we obtain

$$q_{n,n+1} = \frac{1}{24} n! \sum_{q=0}^{n-4} \frac{n^q}{q!} (1 + 5 (n-4-q+1)) \\ = \frac{1}{24} n! \sum_{q=0}^{n-4} \frac{n^q}{q!} (5 (n-q) - 14).$$

Concluding observation. We have made repeated use of the labelled enumeration formula $$B(z) = \frac{A(z)^n}{|G|}$$ which produces the exponential generating function $B(z)$ of objects enumerated by $A(z)$ being distributed into $n$ slots acted on by a permutation group $G,$ where the size of the compound object is the sum of the constituent sizes. This is the lablled counterpart of the Polya Enumeration Theorem.

Addendum. What follows is the Maple code for the above closed formula in terms of the mixed generating function, which is intended for sequence identification and not necessarily for computation, there are recurrences for that, consult e.g. Graphical Enumeration by Harary and Palmer.

with(combinat);

gf :=
proc(n)
option remember;
local subgf;

    subgf := add((1+u)^(m*(m-1)/2)*z^m/m!, m=1..n);

    expand(n!*coeftayl(add((-1)^(q+1)/q*subgf^q, q=1..n),
                       z=0, n));
end;

qval :=
proc(n, k)
option remember;
    coeff(gf(n), u, k);
end;

Addendum. Thu Nov 27 01:12:06 CET 2014. It was pointed out to me in a personal communication that the above closed formula performs extremely poorly where memory and time complexity are concerned. We can however create a very fast equivalent by performing the coefficient extraction for $[z^n]$ before the main computation.

Let us recall the formula: $$q_{n,k} = n! [z^n] [u^k] \sum_{q=1}^n (-1)^{q+1} \frac{1}{q} \left(\sum_{m=1}^n (1+u)^{m(m-1)/2} \frac{z^m}{m!}\right)^q.$$

The key here is to recognize that we are iterating over integer partitions $\lambda\vdash n$ of length $l(\lambda) = q.$ The coefficient on $z^n$ is given by $$\frac{1}{n!} {n\choose \lambda}.$$

The exponent of $(1+u)$ is given by $$\sum_{\lambda_i\in\lambda} \lambda_i(\lambda_i-1)/2.$$

Finally the multiplicity of each partition i.e. the number of corresponding compositions when $\lambda = 1^{f_1} 2^{f_2} 3^{f_3} \cdots$ is $${q\choose f}.$$

This gives the generating function for $n$ fixed which is $$\large{g_n(u) = \sum_{\lambda\vdash n} \frac{(-1)^{q+1}}{q} {n\choose\lambda} {q\choose f} (1+u)^{\sum_{\lambda_i\in\lambda} \lambda_i (\lambda_i-1)/2}}$$ where $q=l(\lambda)$ and $f$ represents the multiplicities in the partition so that $\lambda = 1^{f_1} 2^{f_2} 3^{f_3}\cdots$

This formula would appear to be practical even for large values. Here is the sequence of the enumeration of connected labeled graphs OEIS A001187 up to $n=30,$ computed by Maple almost instantly. These values are obtained by setting $u=1$ in the above formula, which avoids the cost of computing with those polynomials in $u,$ giving $$\large{q_n = \sum_{\lambda\vdash n} \frac{(-1)^{q+1}}{q} {n\choose\lambda} {q\choose f} 2^{\sum_{\lambda_i\in\lambda} \lambda_i (\lambda_i-1)/2}}$$

1,
1,
4,
38,
728,
26704,
1866256,
251548592,
66296291072,
34496488594816,
35641657548953344,
73354596206766622208,
301272202649664088951808,
2471648811030443735290891264,
40527680937730480234609755344896,
1328578958335783201008338986845427712,
87089689052447182841791388989051400978432,
11416413520434522308788674285713247919244640256,
2992938411601818037370034280152893935458466172698624,
1569215570739406346256547210377768575765884983264804405248,
1645471602537064877722485517800176164374001516327306287561310208,
34508369722950116062601714914260936851437546115328069963470233458442...
24,
14473931784581530777452916362195345689326195578125463551466449404195...
748970496,
12141645838784034832247737828641414668703840762841807733278352921867...
1227143860518912,
20370329409143419676922561585800800631483979568699568444273558936889...
94716051486372603625472,
68351532186533737864736355381396298734910952426503780423683990730318...
777915378756861378792989392896,
45869953864873439868450361909803259294922972126320661426113608442339...
62960637520118252235915249481987129344,
61565621838274124223450863197683805128241193119763036274703372417422...
2395343543109861028695816566950855890811486208,
16526397434352809199623091939881315484783346104710447766695225793956...
4080953537482898938408257044203946031706125367800496128,
88725425253946309579607515290733826999038832348034303708272765654674...
479763074364231597119435621862686597717341418971119460584259584,

Maple took $181.311$ seconds (three minutes) to compute the generating function for $n=38$ and $1346$ MB of used memory using the initial unsimplified version of the formula and $4.392$ seconds and $90$ MB using the fast version, a stunning improvement. Only the recurrences by Harary and Palmer and by E.M. Wright will do better.

This was the Maple code.

with(combinat, partition);

gf2 :=
proc(n)
local s, q, p, mcf1, mcf2, li;

    s := 0;

    for p in partition(n) do
        q := nops(p);

        mcf1 := n!/mul(li!, li in p);
        mcf2 :=
        q!/mul(li!, li in
               map(ent -> ent[2], convert(p, multiset)));

        s := s + (-1)^(q+1)/q *
        mcf1*mcf2* (1+u)^add(li*(li-1)/2, li in p);
    od;

    expand(s);
end;

Concluding remark. At this point there is nothing to stop us from extracting the coefficient of $[u^k]$ as well, giving the formula

$$\large{q_{n,k} = \sum_{\lambda\vdash n} \frac{(-1)^{q+1}}{q} {n\choose\lambda} {q\choose f} {\sum_{\lambda_i\in\lambda} \lambda_i (\lambda_i-1)/2 \choose k}}$$

Another optimization. The memory usage of the above is not quite optimal and can be improved by allocating partitions one at a time instead of all at once. This gives the following.

gf2 :=
proc(n)
local s, q, p, mcf1, mcf2, li;
option remember;

    s := 0;

    p := firstpart(n);
    while type(p, list) do
        q := nops(p);

        mcf1 := n!/mul(li!, li in p);
        mcf2 :=
        q!/mul(li!, li in
               map(ent -> ent[2], convert(p, multiset)));

        s := s + (-1)^(q+1)/q *
        mcf1*mcf2* (1+u)^add(li*(li-1)/2, li in p);

        p := nextpart(p);
    od;

    expand(s);
end;

With this version we are only limited by time and not space. Here is the total count for $n=50:$

57775629806264131981532128463353986108213291999872288565750767218860631769...
    6301924134068233518707877841769252356274834883678320922291785288952259...
    3249600859338855724814764410440416662456329476306676699006233890696555...
    2334495222211417966008667425130052344927925607827177068266427605834927...
    5922600493471476178420154378012048571333436567365397136152469165480980...
    158369042006016

Addendum Thu Dec 4 00:48:39 CET 2014. For the purpose of collecting everything in one place let me just briefly comment on those recurrences, following Harary and Palmer.

They use an extremely simple observation namely that if we have two formal power series related by a log-relationship as in

$$\sum_{q\ge 0} a_q z^q = \log\left(\sum_{q\ge 0} A_q z^q \right)$$

then differentiation (which is a standard trick and is often used on recurrences involving the tree function) gives

$$\sum_{q\ge 1} q a_q z^{q-1} = \frac{\sum_{q\ge 1} q A_q z^{q-1}}{\sum_{q\ge 0} A_q z^q}.$$

Re-write this so that the Cauchy product appears more clearly: $$\left(\sum_{q\ge 1} q a_q z^{q-1}\right)\times \left(\sum_{q\ge 0} A_q z^q\right) = \sum_{q\ge 1} q A_q z^{q-1}.$$

Comparing coefficients we obtain for the coefficient on $z^q$ $$\sum_{m=0}^q A_{q-m} (m+1) a_{m+1} = (q+1) A_{q+1}$$

This yields $$(q+1) a_{q+1} = (q+1) A_{q+1} - \sum_{m=0}^{q-1} A_{q-m} (m+1) a_{m+1}$$ or finally for $q\ge 1$ $$a_{q+1} = A_{q+1} - \frac{1}{q+1} \sum_{m=0}^{q-1} A_{q-m} (m+1) a_{m+1}.$$

Note that we are working with exponential generating functions so we are using the values $b_q = q! \times a_q$ and $B_q = q! \times A_q.$ Multiply the recurrence by $(q+1)!$ to get $$b_{q+1} = B_{q+1} - \sum_{m=0}^{q-1} q! \frac{B_{q-m}}{(q-m)!} \frac{b_{m+1}}{m!}.$$ which finally yields $$b_{q+1} = B_{q+1} - \sum_{m=0}^{q-1} {q\choose m} B_{q-m} b_{m+1}.$$

In the present case we have $B_q = (1+u)^{q(q-1)/2}.$ It is important for these recurrences to work that we put $B_0=B_1=b_0=b_1 = 1.$

This recurrence is amazingly fast. With memoization turned on Maple took $58$ seconds to compute $g_{50}(u)$ using the partition iteration method and $1.25$ seconds using the recurrence. For $g_{55}(u)$ the timings were $145$ seconds vs. $1$ second. It does not get any faster than this.

gf3 :=
proc(n)
    option remember;
    local s, B, q;

    if n < 2 then return 1 fi;

    B :=
    proc(m)
        if m < 2 then return 1 fi;
        (1+u)^(m*(m-1)/2);
    end;

    q := n-1;

    s := B(q+1)
    - add(binomial(q,m)*B(q-m)*gf3(m+1),
          m=0..q-1);

    expand(s);
end;

Definitely the last addendum. We can perform coefficient extraction on the terms of this last formula to get

$$q_{n, k} = \begin{cases} 0 \quad\text{if}\quad k\lt n-1 \quad\text{or}\quad k\gt n(n-1)/2 \\ n^{n-2} \quad\text{if}\quad k = n-1, \quad\text{and otherwise}\\ {n(n-1)/2\choose k} - \sum_{m=0}^{n-2} {n-1\choose m} \sum_{p=0}^k {(n-1-m)(n-2-m)/2 \choose p} q_{m+1, k-p}. \end{cases}$$

The complexity of this formula is quite poor as it makes $O(nk)$ recursive calls with $k$ being on average quadratic in $n$. The reader is invited to do the combinatorics and produce a better recurrence. As it stands the fastest version is the one that uses the recurrence derived from the log-relationship, the procedure gf3.

q :=
proc(n, k)
    option remember;
    local res;

    if k < n-1 or k > n*(n-1)/2 then return 0 fi;

    if k = n-1 then return n^(n-2) fi;

    res := binomial(n*(n-1)/2, k)
    - add(binomial(n-1, m)*
          add(binomial((n-1-m)*(n-2-m)/2, p)*q(m+1, k-p),
              p=0..k), m=0..n-2);

    res;
end;

gf4 :=
proc(n)
    option remember;

    add(q(n,k)*u^k, k=n-1..n*(n-1)/2);
end;

As it turns out some optimizations that I thought Maple would do automatically must be implemented manually, giving the following optimized code. This is better than gf4 but still a far cry from what gf3 produces. For that it would appear to need a fundamental change in the mechanism of the recurrence.

qq :=
proc(n, k)
    option remember;
    local res, res1, m, p;

    if k < n-1 or k > n*(n-1)/2 then return 0 fi;

    if k = n-1 then return n^(n-2) fi;

    res := binomial(n*(n-1)/2, k);

    for m from 0 to n-2 do
        res1 := 0;
        for p from max(0, k-1/2*(m+1)*m) to k-m do
            res1 := res1 +
            binomial((n-1-m)*(n-2-m)/2, p)*qq(m+1, k-p)
        od;

        res := res - binomial(n-1, m)*res1;
    od;

    res;
end;

gf5 :=
proc(n)
    option remember;

    add(qq(n,k)*u^k, k=n-1..n*(n-1)/2);
end;

The interested reader may also want to consult this MSE link.

Solution 2:

Since this link has proved somewhat popular but not everybody has access to Maple I am providing a C implementation with a Rosetta stone effect in mind. It uses Linux, GCC and GMP, the GNU multiprecision libary, all of it free software.

This simple implementation cannot compete with Maple's fast polynomial multiplication routines but it easily and dramatically outperforms Maple where memory requirements are considered. This is an implementation of method three from the other answer.

Here is an excerpt from the count of labelled connected graphs on 100 vertices, from an all-inclusive list that was computed in seven minutes of computation time.

000099: 100000000000000000000000000000000000000000000000000000000...
00000000000000000000000000000000000000000000000000000000000000000...
00000000000000000000000000000000000000000000000000000000000000000...
0000000000
000100: 510998031510799015012656647840750649920635707969633111192...
74962605644046874114483034395313967828782400387155843473213972763...
99467183347288626786506867710921922978521653363154616320000000000...
000000000000
000101: 142104312990406781581856724418579090842592500736876137027...
51254362674861198019976979320546703753308719954831818261458369844...
20062161942816326874907532100392442319172473774298628096000000000...
000000000000000
000102: 285447614274319041715048588557494933600690703567446638297...
69028676129587231453229391167155803192891772461953364139194426826...
37181075598419404686535843240073366679296837828275433635840000000...
00000000000000000
000103: 463447763717575519689003062414703595381839972980633113575...
39624468095301976364461153148637638737169543440731383717246460520...
02029075049263459712354346265371921759545922872425765142528000000...
0000000000000000000
000104: 645027341544448563672840738480002899046044088823901830764...
81475379638063891987227335687086646777908138584098153398093224443...
15708495657071771542347812451362360998192317323814303705753190400...
000000000000000000000

And this is the source code.

#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>

typedef struct {
  mpz_t *data;
  int degree;
} pl, *plptr;

plptr pl_new(int degree)
{
  plptr item = (plptr)malloc(sizeof(pl));

  item->degree = degree;
  item->data = (mpz_t *)malloc((degree+1)*sizeof(mpz_t));

  int k;
  for(k = 0; k <= degree; k++){
    mpz_init(item->data[k]);
  }

  return item;
}

void pl_free(plptr item)
{
  int k;

  for(k = 0; k <= item->degree; k++){
    mpz_clear(item->data[k]);
  }

  free(item->data);
  free(item);
}

void pl_print(plptr item, FILE *stream)
{
  int k;

  for(k = 0; k <= item->degree; k++){
    fprintf(stream, "%06d: ", k);
    mpz_out_str(stream, 10, item->data[k]);
    fprintf(stream, "\n");
  }
}

plptr pl_binom(int degree)
{
  plptr item = pl_new(degree);

  int k; mpz_t val;
  mpz_init(val);
  mpz_set_si(val, 1);

  for(k = 0; k <= degree; k++){
    mpz_set(item->data[k], val);

    mpz_mul_si(val, val, degree-k);
    mpz_divexact_ui(val, val, k+1);
  }

  mpz_clear(val);

  return item;
}

void pl_scale(plptr item, mpz_t fact)
{
  int k;

  for(k = 0; k <= item->degree; k++){
    mpz_mul(item->data[k], item->data[k], fact);
  }
}

plptr pl_mul(plptr op1, plptr op2)
{
  int newdeg = op1->degree + op2->degree;
  plptr item = pl_new(newdeg);

  int q; mpz_t term; mpz_t cfprod;
  mpz_init(term); mpz_init(cfprod);
  for(q = 0; q <= newdeg; q++){
    mpz_set_si(term, 0);

    int p;
    for(p = 0; p <= q; p++){
      if(p <= op1->degree && q-p <= op2->degree){
        mpz_mul(cfprod, op1->data[p], op2->data[q-p]);
        mpz_add(term, term, cfprod);
      }
    }

    mpz_set(item->data[q], term);
  }
  mpz_clear(cfprod); mpz_clear(term);

  return item;
}

plptr pl_add(plptr op1, plptr op2)
{
  int newdeg = 
    (op1->degree > op2->degree ?
     op1->degree : op2->degree);
  plptr item = pl_new(newdeg);

  int q; mpz_t term;
  mpz_init(term);
  for(q = 0; q <= newdeg; q++){
    if(q > op1->degree){
      mpz_set(term, op2->data[q]);
    }
    else if(q > op2->degree){
      mpz_set(term, op1->data[q]);
    }
    else{
      mpz_add(term, op1->data[q], op2->data[q]);
    }

    mpz_set(item->data[q], term);
  }
  mpz_clear(term);

  return item;
}

void choose(mpz_t targ, int n, int k)
{
  mpz_set_si(targ, 1);

  int q;
  for(q = 0; q < k; q++){
    mpz_mul_si(targ, targ, n-q);
    mpz_divexact_ui(targ, targ, q+1);
  }
}

int main(int argc, char **argv)
{
  long n = 3, k = -1;

  if(argc >= 2){
    n = atol(argv[1]);

    if(n < 1){
      fprintf(stderr, "vertex count out of range, "
              "got %ld\n", n);
      exit(-1);
    }
  }

  if(argc >= 3){
    k = atol(argv[2]);

    if(k < n-1 || k > n*(n-1)/2){
      fprintf(stderr, "edge count out of range for %ld, "
              "got %ld\n", n, k);
      exit(-2);
    }
  }

  plptr *bintable = 
    (plptr *)malloc((n+1)*sizeof(plptr));

  plptr *gftable = 
    (plptr *)malloc((n+1)*sizeof(plptr));

  int deg;
  for(deg = 0; deg <= n; deg++){
    bintable[deg] = pl_binom(deg*(deg-1)/2);
    gftable[deg] = NULL;
  }

  gftable[0] = pl_binom(0);
  gftable[1] = pl_binom(0);

  int nval; mpz_t factor;

  mpz_init(factor);
  for(nval = 2; nval <= n; nval++){
    int q = nval-1;

    int m; plptr contrib; plptr all;
    all = bintable[q+1];

    for(m = 0; m <= q-1; m++){
      choose(factor, q, m);
      mpz_neg(factor, factor);

      contrib = pl_mul(bintable[q-m], gftable[m+1]);
      pl_scale(contrib, factor);

      plptr prev = all;
      all = pl_add(all, contrib);
      if(m > 0) pl_free(prev);

      pl_free(contrib);
    }

    gftable[nval] = all;
  }
  mpz_clear(factor);

  if(argc == 2){
    pl_print(gftable[n], stdout);
  }
  else{
    mpz_out_str(stdout, 10, gftable[n]->data[k]);
    printf("\n");
  }

  for(deg = 0; deg <= n; deg++){
    pl_free(gftable[deg]);
    pl_free(bintable[deg]);
  }
  free(gftable);
  free(bintable);

  return 0;
}

Solution 3:

Here is my attempt at a Python implementation (2 or 3) of gf3 and gf5. I am only using builtin libraries, so hopefully that will encourage others to play with this. My results agree with the above for $1\leq n \leq30$ and $n=50$ but this will obviously need verification.

The performance of gf5 is not great, over 400 seconds for $g_{55}(u)$ alone. I experimented with a few different external math libraries for calculating the binomial coefficients more quickly, but they actually had surprisingly little effect. I have written the code in such a way that it's easy to substitute but still take advantage of memoization. Perhaps the interested reader could find something better?

gf3, on the other hand, performs as expected. It finished $1\leq n \leq 542$ in 34 seconds. The output is rather cumbersome, but can be viewed here.

from __future__ import division, print_function

from math import factorial
binomial_coefficient_cache = dict()
qq_cache = dict()


def binomial_coefficient_naive(n, k):
    d = n - k
    if d < 0:
        return 0
    return factorial(n) // factorial(k) // factorial(d)
current_binomial = binomial_coefficient_naive


def binomial_memoized(n, k):
    if (n, k) in binomial_coefficient_cache:
        return binomial_coefficient_cache[n, k]
    res = current_binomial(n, k)
    binomial_coefficient_cache[n, k] = res
    return res
binomial = binomial_memoized


def qq(n, k):
    '''Number of labeled, simply connected Graphs of order n, size k '''
    if (n, k) in qq_cache:
        return qq_cache[n, k]
    s = n * (n - 1) // 2
    if k < n - 1 or k > s:
        res = 0
    elif k == n - 1:
        res = int(pow(n, (n - 2)))
    else:
        res = binomial(s, k)
        for m in range(0, n - 1):
            res1 = 0
            lb = max(0, k - (m + 1) * m // 2)
            for p in range(lb, k - m + 1):
                np = (n - 1 - m) * (n - 2 - m) // 2
                res1 += binomial(np, p) * qq(m + 1, k - p)

            res -= binomial(n - 1, m) * res1

    qq_cache[n, k] = res
    return res


def gf5(n):
    '''Number of labeled, simply connected Graphs of order n'''
    ub = (n * (n - 1)) // 2
    qn = sum([qq(n, k) for k in range(n - 1, ub + 1)])
    return(qn)

gf3_cache = dict()
B_cache = dict()


def B(m):
    if m in B_cache:
        return B_cache[m]
    B_cache[m] = int(pow(2, m * (m - 1) // 2)) if m >= 2 else 1
    return B_cache[m]


def gf3(n):
    '''Number of labeled, simply connected Graphs of order n, computed very quickly'''
    if n in gf3_cache:
        return gf3_cache[n]
    if n < 2:
        s = 1
    else:
        q = n - 1
        s = B(q + 1) - sum(binomial(q, m) * B(q - m) * gf3(m + 1)
                           for m in range(q))
    gf3_cache[n] = s
    return s

Any suggestions would be much appreciated, this math is quite a bit over my head ;)