Exponential of a specific hypergeometric series This is motivated by this question.
Let $f$ be the hypergeometric series
$ f(x) = 2 x \, _{4}F_3([1, 1, 4/3, 5/3], [2, 2, 2], 27 x) $
which is explictly given by
$ f(x) = \sum_{n \geq 1} \frac{(3n-1)!}{n!^3} x^n $
and appears in some studies of Mahler measures.

Is is true that $\exp(f)$ has only positive integer coefficients ?

The first few coefficients are $1, 2, 17, 218, 3404, 59644$.
If true, is there a combinatorial meaning for these coefficients ?
In the same vein, what can one say on the coefficients of $\exp(-f)$ which start
$1, -2, -13, -158, -2431, -42454, -802790$ ?
 A: The integrality is a re-incarnation of the problem discussed by Jan Stienstra in "Mahler Measure Variations, Eisenstein Series and Instanton Expansions" (particularly eq. (6) there). The point is that the normalised derivative of the function is an explicit modular form after a suitable modular-function parametrisation, and all this integrality and postivity properties can be extracted from this modular structure on demand. I do not see an easy way to produce an explicit answer to the OP, as this really requires some technical work. But if the author (F.C.) is interested, he can certainly reconstruct the technicalities. I would say that this particular integrality/positivity is less challenging than usual, because we can relate the expansions to modular forms.
A: This is really a comment, but it's too long.
Here's another way to look at this phenomenon that explains the connection with Lyndon words.
Let $L_{\bf n}$ be the number of Lyndon words that are permutations of the multiset $\{1^{n_1}, 2^{n_2},\dots \}$ and let $L(x)=\sum_{\bf n} L_{\bf n} {\bf x^n}$, where 
${\bf x^n} = x_1^{n_2}x_2^{n_2}\cdots$. It is well known that 
$$L(x) = \sum_{n=1}^\infty \frac 1n \sum_{d|n}\mu(d) p_{d}^{n/d},$$
where $p_m=\sum_i x_i^m$ is the power sum symmetric function, though we don't need this formula.
It is also well known that 
$$\frac{1}{1-x_1-x_2-\cdots} = \sum_{n=0}^\infty h_n[L(x)],$$
where $h_n$ is the $n$th complete symmetric function and the square brackets denote composition (plethysm) of symmetric functions, and this yields the formula 
$$\frac{1}{1-x_1-x_2-\cdots} = \prod_{\bf n\ne 0}\frac{1}{(1-{\bf x^n})^{L_{\bf n} }},$$
which may also be seen directly from the Lyndon factorization of words.  From this formula we may extract the diagonal of the logarithm (in as many variables as we like) as above. 
In the case of two variables, it's possible to give a nice combinatorial interpretation in terms of lattice paths to a factorization of $1/(1-x_1-x_2)$ that explains why
$$\exp\biggl(\sum_{n=1}^\infty \frac{1}{2n}\binom{2n}{n}x^n\biggr)$$ is the generating function for the Catalan numbers, and this extends to any power series of the form $1/(1-f(x_1, x_2))$ where $f$ has nonnegative coefficients. (See my paper ￼A factorization for formal power series and lattice path enumeration, J. Combin. Theory Ser. A 28 (1980), 321-337.) But I don't know how to generalize this to more than two variables. 
A: Here's another approach. It resembles Ira's solutions, although less combinatorial. Maybe it will be of use.
Let $a_n := \frac{1}{3} \binom{3n}{n,n,n}$. ($a_1=1$ and $a_n \in \mathbb{Z}$, as $a_n = \binom{3n-1}{n-1,n,n}$.)
Note that $f(x) = \sum_{n \ge 1} \frac{a_n}{n} x^n$. Hence, the function we're dealing with is $\zeta_{\{a_n\}}(x) := \exp( \sum_{n \ge 1} \frac{a_n}{n} x^n)$, which we'll refer to as "the zeta function of $\{a_n\}$".
Lemma 1: Consider a sequence $\{ a_n \}_{n\ge1}$ of rational numbers, with $a_1=1$. $\zeta_{\{a_n\}}(x) \in \mathbb{Z}[[x]]$ iff $\forall n \in \mathbb{N}: n\mid \sum_{d\mid n} \mu(n/d) a_d$.
This lemma reduces our problem to proving a congruence: $\forall n \in \mathbb{N}: 3n\mid \sum_{d\mid n} \mu(n/d) \binom{3d}{d,d,d}$.
Lemma 2: The condition $\forall n \in \mathbb{N}: n\mid \sum_{d\mid n} \mu(n/d) a_d$ is equivalent to $a_{np^k} \equiv a_{np^{k-1}} \mod {p^k}$ for any prime $p$ and $n,k \in \mathbb{N}$.
So we need to prove the congruence $\binom{3np}{np,np,np} \equiv \binom{3n}{n,n,n} \mod {3np}\mathbb{Z}_{p}$ (note how I turned to $p$-adics, as it simplifies computations and presentation).
This congruence is a simple consequence of either:


*

*A combinatorial argument, reminiscent of the combinatorial proof of (a weak version of) Wolstenholme's theorem.

*Properties of the $p$-adic Gamma function. The congruence is equivalent to $\frac{\Gamma_{p}(3np)}{\Gamma_{p}^3(np)} \equiv 1 \mod {3np \mathbb{Z}_{p}}$. Kazandzidis congruences give something much stronger (although a small modification is necessary for the even prime and the second prime):
$$\frac{\Gamma_{p}(3np)}{\Gamma_{p}^3(np)} \equiv 1 \mod {n^3p^3 \mathbb{Z}_{p}}$$


Proofs of both lemmas are not hard. The first lemma follows by writing $\zeta_{a_n}(x)$ as $\prod_{n\ge1} (1-x^n)^{-\frac{b_n}{n}}$, where taking the logarithmic derivative shows $b_n = \sum_{d\mid n} \mu(n/d) a_d$. This gives a proof of one direction. For the other direction, notice that 
$$[x^n]\prod_{i=1}^{\infty} (1-x^i)^{-\frac{b_i}{i}} = [x^n]\prod_{i=1}^{n-1} (1-x^i)^{-\frac{b_i}{i}} + b_n$$
 and use an induction argument. The second lemma is just a "coupling" argument (note how $\mu(ap) = -\mu(a)$ for $p\nmid a$).
A: Here's a sketch of a proof of a much more general result. For simplicity I'll consider only the 3-variable case, but clearly a similar result holds for any number of variables. 
Let $\nabla$ be the operator on formal power series in $x$, $y$, and $z$ defined 
by 
$$\nabla\biggl(\sum_{i,j,k}a_{i,j,k}x^iy^jz^k\biggr) = \sum_i a_{i,i,i}x^i.$$
Then 
$$f(x)  = \nabla\left(\log\frac{1}{1-x-y-z}\right).$$
We have the following result. 
Theorem. Let $U(x,y,z)$ be a formal power series with constant term 1 and integer coefficients. Then $\exp(\nabla(\log U(x,y,z)))$ has integer coefficients.
To prove this we first prove a lemma. 
Lemma. Let $U(x,y,z)$ be a formal power series with constant term 1 and integer coefficients. Then $U(x,y,z)$ can be expressed as an infinite product
$$\prod_{i,j,k}(1+x^iy^jz^k)^{b_{i,j,k}},$$
where the $b_{i,j,k}$ are integers and the product is over all triples $(i,j,k)$ of nonnegative integers, not all zero.
Proof of the lemma (sketch). Equating coefficients of powers of the variables in
$U(x,y,z) = \prod_{i,j,k}(1+x^iy^jz^k)^{b_{i,j,k}}$ gives a recurrence for the $b_{i,j,k}$ that shows that they are integers.
Proof of the theorem (sketch).
Let $$U(x,y,z) = \prod_{i,j,k}(1+x^iy^jz^k)^{b_{i,j,k}}.$$ Then
$$\nabla \log U(x,y,z) = \sum_{i=1}^\infty \log (1+x^i)^{b_{i,i,i}},$$
so
$$
\exp(\nabla(\log U(x,y,z))) = \prod_{i=1}^\infty  (1+x^i)^{b_{i,i,i}},
$$
which has integer coefficients.
Instead of powers of $1+x^iy^jz^k$, we could use powers of $1/(1-x^iy^jz^k)$. In fact, if 
define $c_i$ by 
$$e^{f(x)}=\prod_{i=1}^\infty (1-x^i)^{-c_i},$$
then $c_i$ is the number of Lyndon words with $3n$ beads of 3 colors, $n$ beads of each color (OEIS sequence  A074655.) This gives, in principle, a combinatorial interpretation to the coefficients of $e^f$, though perhaps not a very natural one, and I don't see how it gives a combinatorial interpretation to the coefficients of $1-e^{-f}$.
