Using two different approaches to evaluating the dimensionally regularized ($d=4-2\epsilon$ dimensional Euclidean space), equal mass ($x=m^2$), 2-loop vacuum Feynman diagram
$$
\begin{align}
I(x) &= \int\frac{\mathrm{d}^dp}{\pi^{d/2}}\frac{\mathrm{d}^dk}{\pi^{d/2}}
      \frac1{(k^2+x)(p^2+x)((k+p)^2+x)} \\\
     &= \int _0^{\infty }\int _0^{\infty }\int _0^{\infty }
      \frac{e^{-x(s_1+s_2+s_3)}}{\left(s_1s_2+s_2s_3+s_3s_1\right)^{d/2}}
      \mathrm{d}s_1\mathrm{d}s_2\mathrm{d}s_3 \\\ 
     &= x^{d-3}\\,\Gamma\left(2-\frac{d}{2}\right)\Gamma\left(1-\frac{d}{2}\right)
        \\,f(d)\\,,
\end{align}
$$
the following hypergeometric identity arises
$$
\begin{align}
f(d) &=f_1(d) = 
2\\, {}_2F_1\left(1,\frac{d-1}{2};\frac{3}{2};-\frac{1}{3}\right)
-4^{2-d} 3^{(d - 3)/2} B\left(\frac{3 - d}2, \frac{3 - d}2\right) \\\
&= f_2(d) = \frac{4}{3} \left(
{}_2F_1\left(1,\frac{d-1}{2};\frac{3}{2};-\frac{1}{3}\right)+\frac{1 }{d-3}  {}_2F_1\left(1,\frac{d-1}{2};\frac{5-d}{2};-\frac{1}{3}\right) \right)\\,,
\end{align}
$$
where the second term in $f_1(d)$ can be reduced with the identity 
$B(x,x)=2^{1-2x}B(x,1/2)$.

The identity $f_1(d)=f_2(d)$ checks out numerically and (provided no mistakes have been made in the calculations) it should be identically true. So far I have been unable to find a direct proof of the identity.

Can anyone here prove this identity or suggest a good approach?
A computer proof (using human checkable code/steps) is acceptable.

----

For convenience I've provided the Mathematica `InputForm` of the two functions:

    f1[d_] := 2 Hypergeometric2F1[1, (d - 1)/2, 3/2, -1/3] - 
              2^(4 - 2 d) 3^((d - 3)/2) Beta[(3 - d)/2, (3 - d)/2]

    f2[d_] := 4/3 (Hypergeometric2F1[1, (d - 1)/2, 3/2, -1/3] + 
              1/(d - 3) Hypergeometric2F1[1, (d - 1)/2, (5 - d)/2, -1/3])

----

**Aside:**    
$f_1(d)$ comes from direct integration using Feynman parameters (my own calculation, I don't know of a reference that includes it).   
$f_2(d)$ comes from direct integration using the Mellin-Barnes representation (the result presented above is a version of eqn(33) of [hep-ph/9304303](http://arxiv.org/abs/hep-ph/9304303), see also references within) .

----

**Edit:**
I just noticed [this MO answer](http://mathoverflow.net/questions/37787/proving-a-hypergeometric-function-identity/37836#37836) that mentions the [HolonomicFunctions](http://www.risc.jku.at/research/combinat/software/HolonomicFunctions/index.php) package for Mathematica. It shows that both functions obey the recursion     
$(4+4 d)f_i(d+4)+(4-7 d)f_i(d+2)+(-6+3 d)f_i(d)=0$,     
but of course, the integral diverges for integer $d\geq4$ and I need to prove the relation for all $d<4$ (dimensional reduction) or for all complex $d$ (dimensional regularization).