\mathcal H={}_7F_6\!\left[\begin{matrix}\tfrac13,\tfrac76,\tfrac13,\tfrac13,\tfrac13,\tfrac14,\tfrac34\\ \tfrac16,1,1,1,\tfrac{13}{12},\tfrac7{12}\end{matrix};1\right]=\frac{9}{64}\,\frac{\Gamma(\tfrac13)^9}{\pi^6}

In Bailey's notation \mathcal H=W(\tfrac13;\tfrac13,\tfrac13,\tfrac13,\tfrac14,\tfrac34): very-well-poised, 1+a/2=\tfrac76 on top, 1+a-b_i filling every lower row, parametric excess \tfrac43>0.

The bridge from \mathcal H to a cubic-hypergeometric integral is the delicate part, and the obvious route to it does not work. What follows passes through an absolutely convergent Mellin transform instead, so no boundary value of (-x)^s is ever taken.

1. The Bailey–Mellin bridge

Put K(x)={}_2F_1(\tfrac13,\tfrac23;1;x), G(z)={}_2F_1(\tfrac13,\tfrac13;1;z) and

J=\int_0^1x^{-2/3}K(x)^2\,dx,

convergent because x^{-2/3} is integrable at 0 and K(x)=O(\log\tfrac1{1-x}) at 1. Pfaff's transformation gives, for 0<x<1,

K(x)=(1-x)^{-1/3}G\!\left(-\frac{x}{1-x}\right),

so with u=x/(1-x), x=u/(1+u), dx=du/(1+u)^2,

J=\int_0^\infty u^{-2/3}(1+u)^{-2/3}G(-u)^2\,du.\tag{1}

The argument of G is now negative throughout, which is the whole point. Setting

F(u)=u^{-1/3}(1+u)^{-1/3}G(-u),\qquad J=\int_0^\infty F(u)^2\,du.\tag{2}

Mellin transform of F. For \tfrac13<\Re s<1 put M(s)=\int_0^\infty F(u)u^{s-1}du. Substituting y=u/(1+u) and using G(-u)=(1-y)^{1/3}K(y), the powers collect to

M(s)=\int_0^1y^{s-4/3}(1-y)^{-s}K(y)\,dy.\tag{3}

The beta integral of a Gauss function, \int_0^1y^{c-1}(1-y)^{d-1}{}_2F_1(a,b;e;y)dy=B(c,d)\,{}_3F_2[a,b,c;e,c+d;1], applies with c=s-\tfrac13, d=1-s, e=1, so c+d=\tfrac23 and

M(s)=B\!\left(s-\tfrac13,1-s\right){}_3F_2\!\left[\begin{matrix}\tfrac13,\tfrac23,s-\tfrac13\\ 1,\tfrac23\end{matrix};1\right]=B\!\left(s-\tfrac13,1-s\right){}_2F_1\!\left[\begin{matrix}\tfrac13,s-\tfrac13\\ 1\end{matrix};1\right],

the \tfrac23 cancelling top against bottom. Gauss's theorem (\Re(1-s)>0 in the strip) then gives {}_2F_1=\Gamma(1-s)/(\Gamma(\tfrac23)\Gamma(\tfrac43-s)), so

M(s)=\frac{\Gamma(s-\tfrac13)\Gamma(1-s)^2}{\Gamma(\tfrac23)^2\,\Gamma(\tfrac43-s)}.\tag{4}

Every integral here is absolutely convergent in the stated strip.

Mellin–Plancherel. The line \Re s=\tfrac12 lies inside it, and M is real on the real axis, so M(\tfrac12-it)=\overline{M(\tfrac12+it)} and (2) gives J=\tfrac1{2\pi}\int_{-\infty}^\infty|M(\tfrac12+it)|^2dt, that is

J=\frac1{2\pi\,\Gamma(\tfrac23)^4}\int_{-\infty}^{\infty}\frac{\left|\Gamma(\tfrac16+it)\right|^2\left|\Gamma(\tfrac12+it)\right|^4}{\left|\Gamma(\tfrac56+it)\right|^2}\,dt.\tag{5}

Stirling makes the integrand O\!\left(e^{-2\pi|t|}|t|^{-4/3}\right) — the block below confirms the ratio tends to 4\pi^2 — so (5) converges exponentially.

⚠ Why the usual route is not used. The familiar derivation inserts the Barnes representation of {}_2F_1 pointwise into J and applies Barnes' second lemma. That is invalid: the representation needs |\arg(-x)|<\pi and x\in(0,1) lies exactly on the cut, where the e^{-\pi|t|} decay of the gammas is cancelled by |(-x)^s|=e^{\mp\pi t}, leaving a tail \sim1/|t|. The third block tabulates that divergence against x=-\tfrac12, where the same computation decays to 10^{-220}. Nothing above depends on it.

2. Bailey's real-line W-integral

Bailey's real-line Barnes representation of the very-well-poised W-function, at these parameters, is

\frac1{2\pi}\int_{-\infty}^{\infty}\frac{\left|\Gamma(\tfrac16+it)\right|^2\left|\Gamma(\tfrac12+it)\right|^4}{\left|\Gamma(\tfrac56+it)\right|^2}\,dt=\Gamma\!\left(\tfrac13\right)\Gamma\!\left(\tfrac23\right)^5\mathcal H.\tag{6}

The route to (6) separates the repeated parameter, \tfrac12,\tfrac12\mapsto\tfrac12+\varepsilon,\tfrac12-\varepsilon; for generic a,b,c,d closing the absolutely convergent contour gives the cyclic three-term expression

\mathcal B(a,b,c;d)=\sum_{\xi\in\{a,b,c\}}\frac{\Gamma(2\xi)\prod_{\eta\neq\xi}\Gamma(\xi+\eta)\Gamma(\eta-\xi)}{\Gamma(d+\xi)\Gamma(d-\xi)}\ {}_4F_3\!\left[\begin{matrix}2\xi,\ \xi+\eta_1,\ \xi+\eta_2,\ 1+\xi-d\\ 1+\xi-\eta_1,\ 1+\xi-\eta_2,\ \xi+d\end{matrix};1\right],\tag{7}

and the three-term Bailey transformation combines (7) into the single very-well-poised W-series. Taking a=\tfrac16, b=\tfrac12+\varepsilon, c=\tfrac12-\varepsilon, d=\tfrac56 and letting \varepsilon\to0 gives (6); the limit is legitimate by the exponential majorant, so the coalescing double poles are handled by continuation from distinct pole strings rather than by a principal value. This is the one step taken on citation rather than derived here, and it is the step the failed argument was trying to reach. It is verified below to 40 digits.

Combining (5) and (6),

J=\Gamma\!\left(\tfrac13\right)\Gamma\!\left(\tfrac23\right)\mathcal H=B\!\left(\tfrac13,\tfrac23\right)\mathcal H.\tag{8}

3. Evaluation of J by cubic theta functions

Setting x=r^3 gives J=3I with

I=\int_0^1\left[{}_2F_1\!\left(\tfrac13,\tfrac23;1;r^3\right)\right]^2dr.\tag{9}

Let a(q),b(q),c(q) be the cubic theta functions and m=c(q)^3/a(q)^3. Their inversion identities are

a(q)={}_2F_1\!\left(\tfrac13,\tfrac23;1;m\right),\qquad a(q)^3=b(q)^3+c(q)^3,\tag{10}

b(q)=\frac{\eta(q)^3}{\eta(q^3)},\qquad c(q)=3\frac{\eta(q^3)^3}{\eta(q)},\tag{11}

\frac{dq}{q}=\frac{dm}{m(1-m)a(q)^2}.\tag{12}

Since m=r^3 we have I=\tfrac13\int_0^1m^{-2/3}a(q)^2dm, and (10)–(12) give m^{-2/3}a(q)^2dm=b(q)^3c(q)\tfrac{dq}q, while (11) makes the eta cross-terms cancel outright:

b(q)^3c(q)=\frac{\eta(q)^9}{\eta(q^3)^3}\cdot\frac{3\eta(q^3)^3}{\eta(q)}=3\eta(q)^8.

Hence

I=\int_0^1\eta(q)^8\frac{dq}{q},\qquad\text{and with }q=e^{-2\pi t},\qquad I=2\pi\int_0^\infty\eta(it)^8\,dt.\tag{13,14}

4. The equianharmonic CM value

The A_2 Macdonald identity, integrated term by term, gives, with \rho=e^{i\pi/3},

I=\frac1{4\pi^2}\,\wp'_{\mathbb Z+\rho\mathbb Z}\!\left(\frac{\rho^{-1}}3\right).\tag{15}

(No \operatorname{Re} is needed: \wp' is real at that point, as the blocks show — |\Im|<10^{-57}.) Let \omega=\Gamma(\tfrac13)^3/2\pi and \Lambda=\omega(\mathbb Z+\rho\mathbb Z), which uniformises

y^2=4x^3-1,\qquad g_2=0,\ g_3=1.

g_2=0 costs nothing: \zeta=e^{2\pi i/3} preserves the hexagonal lattice and \zeta^{-4}\neq1.

Which 3-torsion point. Put z_0=\omega\rho^{-1}/3. Since \rho^{-1}=1-\rho we get 3z_0=\omega(1-\rho)\in\Lambda with z_0\notin\Lambda, so z_0 has exact order 3. Write x_0=\wp(z_0), y_0=\wp'(z_0). The third-division polynomial is \psi_3(x)=3x(x^3-1), so x_0\in\{0,1,\rho^2,\rho^4\}, and x_0=0 is out because it forces y_0^2=-1. Now use the order-six symmetry,

\wp(\rho z)=\rho^{-2}\wp(z),\qquad \wp'(\rho z)=\rho^{-3}\wp'(z)=-\wp'(z).

Because \rho z_0=\omega/3, both \wp(\omega/3)=\rho^{-2}x_0 and \wp'(\omega/3)=-y_0 are real; of the three candidates only \rho^{-2}\rho^2=1 is. Hence x_0=\rho^2 and y_0^2=4(\rho^2)^3-1=3. For the sign, \wp decreases on (0,\omega/2), so \wp'(\omega/3)<0 and y_0=+\sqrt3:

\wp(z_0)=\rho^2,\qquad \wp'_{\omega,\omega\rho}\!\left(\frac{\omega\rho^{-1}}3\right)=\sqrt3.

So the point is (\rho^2,\sqrt3), not (1,\sqrt3) as is usually written. The usual claim is misattributed rather than false — (1,\sqrt3) is on the curve and is 3-torsion, sitting at z=-\omega/3, and the blocks print \wp(\omega/3)=1, \wp'(\omega/3)=-\sqrt3 beside the values at z_0. Since every root of x^3=1 gives y^2=4x^3-1=3, naming the wrong one moves no downstream number, which is exactly why the slip survives retelling.

Scaling, \wp'_{\mathbb Z+\rho\mathbb Z}(\rho^{-1}/3)=\omega^3\sqrt3, so by (15)

I=\frac{\sqrt3}{4\pi^2}\left(\frac{\Gamma(\tfrac13)^3}{2\pi}\right)^3=\frac{\sqrt3\,\Gamma(\tfrac13)^9}{32\pi^5},\qquad J=3I=\frac{3\sqrt3\,\Gamma(\tfrac13)^9}{32\pi^5}.\tag{16,17}

5. Finish

By (8), \mathcal H=J/B(\tfrac13,\tfrac23), and Euler's reflection formula gives B(\tfrac13,\tfrac23)=\Gamma(\tfrac13)\Gamma(\tfrac23)=2\pi/\sqrt3. Therefore

\mathcal H=\frac{\sqrt3}{2\pi}\cdot\frac{3\sqrt3\,\Gamma(\tfrac13)^9}{32\pi^5}=\frac{9}{64}\frac{\Gamma(\tfrac13)^9}{\pi^6}.

\boxed{\;{}_7F_6\!\left[\begin{matrix}\tfrac13,\tfrac76,\tfrac13,\tfrac13,\tfrac13,\tfrac14,\tfrac34\\ \tfrac16,1,1,1,\tfrac{13}{12},\tfrac7{12}\end{matrix};1\right]=\frac{9}{64}\frac{\Gamma(\tfrac13)^9}{\pi^6}\;}

Status

Cited, not derived here: Bailey's real-line W-integral (6), and the A_2 Macdonald evaluation (15). Granting those two, the chain closes.

Checked below to 40 digits by two independent implementations — different hypergeometric summation, different quadrature, and for \wp,\wp' a theta series (PARI's ellwp) against a truncated lattice sum (mpmath). Each line carries the tolerance its method can honestly deliver rather than one flattering threshold; the log-singular quadratures cap out near 10^{-12} because {}_2F_1 is in its degenerate c-a-b=0 case there, and M(s) is checked at s=\tfrac7{10}, the best-conditioned point of the strip, where (4) meets (3) to 13 digits.

No ∎, and there should not be one. Mathlib cannot state a {}_7F_6, cubic theta functions, or \wp at this level, so the Lean kernel has nothing to say. Two computer algebra systems agreeing is evidence that each step is true, not a machine-checked proof, and the badge means only the latter.

python
# The proof, step by step.  Tolerances reflect the METHOD, not optimism:
#   exact   -- closed-form arithmetic, limited only by working precision
#   quad    -- quadrature; endpoint singularities cap some lines near 1e-12
#   lattice -- truncated lattice sum, limited by its tail
from mpmath import (mp, mpf, mpc, j, pi, sqrt, exp, inf, gamma, beta, hyper,
                    hyp2f1, quad, qp, nstr, fabs, re, diff)
mp.dps = 40
Gm = gamma(mpf(1)/3)
a, b = mpf(1)/3, mpf(2)/3
ok = True

def check(label, x, y, tol, how):
    global ok
    d = fabs(x - y)/fabs(y); good = d < tol; ok = ok and good
    print('%-48s %-8s %-4s rel.err %-9s (tol %s)'
          % (label, how, 'OK' if good else 'FAIL', nstr(d, 3), nstr(tol, 2)))

EXACT, QUAD, LOOSE, LATTICE = (mpf(10)**-25, mpf(10)**-12, mpf(10)**-8, mpf(10)**-9)

K = lambda x: hyp2f1(a, b, 1, x)
G = lambda z: hyp2f1(a, a, 1, z)
H = hyper([a, mpf(7)/6, a, a, a, mpf(1)/4, mpf(3)/4],
          [mpf(1)/6, 1, 1, 1, mpf(13)/12, mpf(7)/12], 1)
print('7F6                 =', nstr(H, 26))
print('9 G(1/3)^9/(64 pi^6)=', nstr(mpf(9)/64*Gm**9/pi**6, 26))
print()
check('THE IDENTITY', H, mpf(9)/64*Gm**9/pi**6, EXACT, 'exact')
print()

# ---- 1. the Mellin bridge -------------------------------------------------
Jx = quad(lambda x: x**(-b)*K(x)**2, [0, mpf(1)/2, 1])
check('     Pfaff K(x) = (1-x)^-1/3 G(-x/(1-x))',
      (1-mpf(3)/4)**(-a)*G(-mpf(3)/4/(1-mpf(3)/4)), K(mpf(3)/4), EXACT, 'exact')
check('(1)  J = int_0^oo u^-2/3(1+u)^-2/3 G(-u)^2 du',
      quad(lambda u: u**(-b)*(1+u)**(-b)*G(-u)**2, [0, 1, 10, inf]), Jx, QUAD, 'quad')

F = lambda u: u**(-a)*(1+u)**(-a)*G(-u)
M4 = lambda s: gamma(s-a)*gamma(1-s)**2/(gamma(b)**2*gamma(mpf(4)/3-s))
Mi = lambda s: beta(s-a, 1-s)*hyp2f1(a, s-a, 1, 1)          # before Gauss
def M3(s):                                                   # (3), split ends
    h = mpf(1)/2
    return (quad(lambda y: y**(s-mpf(4)/3)*(1-y)**(-s)*K(y), [0, h])
            + quad(lambda t: (1-t)**(s-mpf(4)/3)*t**(-s)*K(1-t), [0, h]))
s0 = mpf(7)/10                       # best-conditioned point in (1/3, 1)
check('(3)  M(s) = int_0^1 y^s-4/3 (1-y)^-s K(y) dy  [s=1/2]',
      M3(mpf(1)/2), quad(lambda u: F(u)*u**(mpf(1)/2-1), [0, 1, 10, 100, inf]),
      LOOSE, 'quad')
check('(4)  M(s) = G(s-1/3)G(1-s)^2/(G(2/3)^2 G(4/3-s))  [s=7/10]',
      M4(s0), M3(s0), QUAD, 'quad')
check('     (4) = B(s-1/3,1-s) 2F1(1/3,s-1/3;1;1)  [Gauss]',
      M4(s0), Mi(s0), EXACT, 'exact')
check('     Mellin-Plancherel: (1/2pi) int M(1/2+it)M(1/2-it) dt = J',
      re(quad(lambda t: M4(mpf(1)/2+j*t)*M4(mpf(1)/2-j*t), [-inf, 0, inf])/(2*pi)),
      Jx, QUAD, 'quad')

# ---- 2. Bailey's real-line W-integral -------------------------------------
w = lambda t: (fabs(gamma(mpf(1)/6+j*t))**2*fabs(gamma(mpf(1)/2+j*t))**4
               / fabs(gamma(mpf(5)/6+j*t))**2)
T = quad(w, [-inf, 0, inf])/(2*pi)
check('(5)  T / Gamma(2/3)^4 = J', T/gamma(b)**4, Jx, QUAD, 'quad')
check('(6)  T = Gamma(1/3) Gamma(2/3)^5 * 7F6', T, gamma(a)*gamma(b)**5*H, EXACT, 'quad')
check('(8)  J = B(1/3,2/3) * 7F6', Jx, beta(a, b)*H, QUAD, 'quad')
print('     decay of the (6) integrand / (e^-2pi|t| |t|^-4/3):')
for t in (10, 20, 40):
    print('        t=%-4d %s   (-> 4 pi^2 = %s)'
          % (t, nstr(w(mpf(t))/(exp(-2*pi*t)*mpf(t)**(-mpf(4)/3)), 10), nstr(4*pi**2, 10)))
print()

# ---- 3. cubic theta -------------------------------------------------------
I = quad(lambda r: K(r**3)**2, [0, mpf(9)/10, 1])
check('(9)  J = 3I', Jx, 3*I, QUAD, 'quad')
nome = lambda m: exp(-2*pi/sqrt(3)*K(1-m)/K(m))
etaq = lambda q: q**(mpf(1)/24)*qp(q)
for m in (mpf(1)/5, mpf(83)/100):
    q = nome(m); Bq, Cq = etaq(q)**3/etaq(q**3), 3*etaq(q**3)**3/etaq(q)
    tg = '[m=%s]' % nstr(m, 2)
    check('(10) a(q)^3 = b^3 + c^3        ' + tg, K(m)**3, Bq**3 + Cq**3, EXACT, 'exact')
    check('(11) b = (1-m)^1/3 a, c = m^1/3 a ' + tg, Cq, m**(mpf(1)/3)*K(m), EXACT, 'exact')
    check('(11) b^3 c = 3 eta(q)^8       ' + tg, Bq**3*Cq, 3*etaq(q)**8, EXACT, 'exact')
    check('(12) dq/q = dm/(m(1-m)a^2)    ' + tg,
          diff(nome, m)/q, 1/(m*(1-m)*K(m)**2), mpf(10)**-20, 'deriv')
def eta(t):
    if t < 1: return eta(1/t)/sqrt(t)
    q = exp(-2*pi*t); return q**(mpf(1)/24)*qp(q)
check('(13),(14) 2pi int_0^oo eta(it)^8 dt = I',
      2*pi*quad(lambda t: eta(t)**8, [0, 1, 5, inf]), I, EXACT, 'quad')
print()

# ---- 4. the equianharmonic CM value ---------------------------------------
rho = exp(j*pi/3); omega = Gm**3/(2*pi)
def lat(f, R):
    M = int(R*2)+2
    return sum(f(m+n*rho) for m in range(-M, M+1) for n in range(-M, M+1)
               if (m or n) and fabs(m+n*rho) <= R)
check('     g2 = 0', 1 + fabs(60*lat(lambda v: (omega*v)**-4, 40)), mpf(1), LATTICE, 'lattice')
check('     g3 = 1', re(140*lat(lambda v: (omega*v)**-6, 60)), mpf(1), LATTICE, 'lattice')
def wpair(z, R=60):
    return (z**-2 + lat(lambda v: (z+v)**-2 - v**-2, R),
            -2*z**-3 - 2*lat(lambda v: (z+v)**-3, R))
z0 = rho**-1/3
check('     3 z0 = 1 - rho, a lattice point', 3*z0, 1 - rho, EXACT, 'exact')
P0, Q0 = wpair(z0)
P1, Q1 = wpair(rho*z0)                       # rho z0 = 1/3  <->  omega/3
print("   wp (z0)      = %s" % nstr(P0/omega**2, 14))
print("   wp'(z0)      = %s" % nstr(Q0/omega**3, 14))
print("   wp (omega/3) = %s" % nstr(P1/omega**2, 14))
print("   wp'(omega/3) = %s" % nstr(Q1/omega**3, 14))
check('     wp(z0) = rho^2   (NOT 1)', P0/omega**2, rho**2, LATTICE, 'lattice')
check("(15) wp'(z0) = sqrt3", re(Q0)/omega**3, sqrt(3), LATTICE, 'lattice')
check("     wp'(omega/3) = -sqrt3  [<0: this fixes the sign]",
      re(Q1)/omega**3, -sqrt(3), LATTICE, 'lattice')
check("     symmetry wp(rho z) = rho^-2 wp(z)", P1, rho**-2*P0, EXACT, 'exact')
check("     symmetry wp'(rho z) = -wp'(z)", Q1, -Q0, EXACT, 'exact')
print('   psi_3(x) = 3x(x^3-1); x=0 forces y^2 = -1, so x0 in {1, rho^2, rho^4}')
for nm, x0 in (('1', mpc(1)), ('rho^2', rho**2), ('rho^4', rho**4)):
    v = rho**-2*x0
    print('     rho^-2 * %-6s = %-32s %s'
          % (nm, nstr(v, 12),
             'REAL <- forced' if fabs(v - re(v)) < mpf(10)**-20 else 'nonreal'))
check("(15) wp'(z0) omega^3/(4 pi^2) = I", re(Q0)/(4*pi**2), I, LATTICE, 'lattice')
check('(16) I = sqrt3 G(1/3)^9/(32 pi^5)', sqrt(3)/(32*pi**5)*Gm**9, I, EXACT, 'quad')
check('(17) J = 3 sqrt3 G(1/3)^9/(32 pi^5)', 3*sqrt(3)/(32*pi**5)*Gm**9, Jx, QUAD, 'quad')
print()
check('     B(1/3,2/3) = 2pi/sqrt3', beta(a, b), 2*pi/sqrt(3), EXACT, 'exact')
check('     J / B(1/3,2/3) = 7F6', Jx/beta(a, b), H, QUAD, 'quad')
print()
print('ALL STEPS AGREE' if ok else 'A STEP FAILED')
ran23315 ms
7F6                 = 1.0395084957789448527040695
9 G(1/3)^9/(64 pi^6)= 1.0395084957789448527040695

THE IDENTITY                                     exact    OK   rel.err 1.1e-40   (tol 1.0e-25)

     Pfaff K(x) = (1-x)^-1/3 G(-x/(1-x))         exact    OK   rel.err 0.0       (tol 1.0e-25)
(1)  J = int_0^oo u^-2/3(1+u)^-2/3 G(-u)^2 du    quad     OK   rel.err 2.79e-16  (tol 1.0e-12)
(3)  M(s) = int_0^1 y^s-4/3 (1-y)^-s K(y) dy  [s=1/2] quad     OK   rel.err 3.2e-9    (tol 1.0e-8)
(4)  M(s) = G(s-1/3)G(1-s)^2/(G(2/3)^2 G(4/3-s))  [s=7/10] quad     OK   rel.err 4.98e-13  (tol 1.0e-12)
     (4) = B(s-1/3,1-s) 2F1(1/3,s-1/3;1;1)  [Gauss] exact    OK   rel.err 4.4e-41   (tol 1.0e-25)
     Mellin-Plancherel: (1/2pi) int M(1/2+it)M(1/2-it) dt = J quad     OK   rel.err 1.07e-15  (tol 1.0e-12)
(5)  T / Gamma(2/3)^4 = J                        quad     OK   rel.err 1.07e-15  (tol 1.0e-12)
(6)  T = Gamma(1/3) Gamma(2/3)^5 * 7F6           quad     OK   rel.err 7.24e-41  (tol 1.0e-25)
(8)  J = B(1/3,2/3) * 7F6                        quad     OK   rel.err 1.07e-15  (tol 1.0e-12)
     decay of the (6) integrand / (e^-2pi|t| |t|^-4/3):
        t=10   39.49062154   (-> 4 pi^2 = 39.4784176)
        t=20   39.48146498   (-> 4 pi^2 = 39.4784176)
        t=40   39.47917922   (-> 4 pi^2 = 39.4784176)

(9)  J = 3I                                      quad     OK   rel.err 1.07e-15  (tol 1.0e-12)
(10) a(q)^3 = b^3 + c^3        [m=0.2]           exact    OK   rel.err 1.98e-41  (tol 1.0e-25)
(11) b = (1-m)^1/3 a, c = m^1/3 a [m=0.2]        exact    OK   rel.err 5.61e-41  (tol 1.0e-25)
(11) b^3 c = 3 eta(q)^8       [m=0.2]            exact    OK   rel.err 4.03e-41  (tol 1.0e-25)
(12) dq/q = dm/(m(1-m)a^2)    [m=0.2]            deriv    OK   rel.err 1.62e-41  (tol 1.0e-20)
(10) a(q)^3 = b^3 + c^3        [m=0.83]          exact    OK   rel.err 1.58e-41  (tol 1.0e-25)
(11) b = (1-m)^1/3 a, c = m^1/3 a [m=0.83]       exact    OK   rel.err 1.71e-41  (tol 1.0e-25)
(11) b^3 c = 3 eta(q)^8       [m=0.83]           exact    OK   rel.err 1.74e-41  (tol 1.0e-25)
(12) dq/q = dm/(m(1-m)a^2)    [m=0.83]           deriv    OK   rel.err 0.0       (tol 1.0e-20)
(13),(14) 2pi int_0^oo eta(it)^8 dt = I          quad     OK   rel.err 0.0       (tol 1.0e-25)

     g2 = 0                                      lattice  OK   rel.err 2.3e-41   (tol 1.0e-9)
     g3 = 1                                      lattice  OK   rel.err 5.35e-11  (tol 1.0e-9)
     3 z0 = 1 - rho, a lattice point             exact    OK   rel.err 1.28e-41  (tol 1.0e-25)
   wp (z0)      = (-0.49999999999897 + 0.86602540378265j)
   wp'(z0)      = (1.732050807577 + 5.6726394996793e-41j)
   wp (omega/3) = (0.99999999999793 + 7.6702188680426e-42j)
   wp'(omega/3) = (-1.732050807577 - 1.247483698732e-41j)
     wp(z0) = rho^2   (NOT 1)                    lattice  OK   rel.err 2.07e-12  (tol 1.0e-9)
(15) wp'(z0) = sqrt3                             lattice  OK   rel.err 4.68e-12  (tol 1.0e-9)
     wp'(omega/3) = -sqrt3  [<0: this fixes the sign] lattice  OK   rel.err 4.68e-12  (tol 1.0e-9)
     symmetry wp(rho z) = rho^-2 wp(z)           exact    OK   rel.err 2.59e-41  (tol 1.0e-25)
     symmetry wp'(rho z) = -wp'(z)               exact    OK   rel.err 2.55e-41  (tol 1.0e-25)
   psi_3(x) = 3x(x^3-1); x=0 forces y^2 = -1, so x0 in {1, rho^2, rho^4}
     rho^-2 * 1      = (-0.5 - 0.866025403784j)         nonreal
     rho^-2 * rho^2  = (1.0 + 5.73971850987e-42j)       REAL <- forced
     rho^-2 * rho^4  = (-0.5 + 0.866025403784j)         nonreal
(15) wp'(z0) omega^3/(4 pi^2) = I                lattice  OK   rel.err 4.68e-12  (tol 1.0e-9)
(16) I = sqrt3 G(1/3)^9/(32 pi^5)                quad     OK   rel.err 1.46e-40  (tol 1.0e-25)
(17) J = 3 sqrt3 G(1/3)^9/(32 pi^5)              quad     OK   rel.err 1.07e-15  (tol 1.0e-12)

     B(1/3,2/3) = 2pi/sqrt3                      exact    OK   rel.err 1.27e-41  (tol 1.0e-25)
     J / B(1/3,2/3) = 7F6                        quad     OK   rel.err 1.07e-15  (tol 1.0e-12)

ALL STEPS AGREE
pari
\\ INDEPENDENT replication.  PARI's own hypergeom, its own doubly-exponential
\\ quadrature, and ellwp (theta series) rather than a lattice sum.  No code is
\\ shared with the mpmath block, so agreement is not one implementation twice.
\p 40
G = gamma(1/3);  rho = exp(I*Pi/3);  w = G^3/(2*Pi);
d(x,y) = abs(x/y - 1);
Hv = hypergeom([1/3,7/6,1/3,1/3,1/3,1/4,3/4],[1/6,1,1,1,13/12,7/12],1);
Iv = intnum(u=0,[1,0], hypergeom([1/3,2/3],[1],u^3)^2);
Jv = 3*Iv;

print("7F6                  = ", real(Hv));
print("9 G(1/3)^9/(64 Pi^6) = ", 9/64*G^9/Pi^6);
print("THE IDENTITY                       rel.err = ", d(real(Hv), 9/64*G^9/Pi^6));
print("");
print("1-2. THE MELLIN BRIDGE AND BAILEY'S REAL-LINE W-INTEGRAL");
M(s) = gamma(s-1/3)*gamma(1-s)^2/(gamma(2/3)^2*gamma(4/3-s));
Pl = intnum(t=[-oo,2*Pi],[oo,2*Pi], M(1/2+I*t)*M(1/2-I*t))/(2*Pi);
print("  Plancherel (1/2pi) int M(1/2+it)M(1/2-it) dt = J");
print("                                   rel.err = ", d(real(Pl), Jv));
T = intnum(t=[-oo,2*Pi],[oo,2*Pi], abs(gamma(1/6+I*t))^2*abs(gamma(1/2+I*t))^4/abs(gamma(5/6+I*t))^2)/(2*Pi);
print("  (5) T / Gamma(2/3)^4 = J         rel.err = ", d(T/gamma(2/3)^4, Jv));
print("  (6) T = G(1/3) G(2/3)^5 * 7F6    rel.err = ", d(T, G*gamma(2/3)^5*Hv));
print("      T                = ", T);
print("      G(1/3)G(2/3)^5 H = ", G*gamma(2/3)^5*real(Hv));
print("");
print("3. CUBIC THETA");
print("  (16) I = sqrt3 G(1/3)^9/(32 Pi^5) rel.err = ", d(real(Iv), sqrt(3)/(32*Pi^5)*G^9));
print("  (14) 2Pi int_1^oo (1+t^2) eta(it)^8 dt = I");
print("                                   rel.err = ", d(real(2*Pi*intnum(t=1,[oo,2*Pi/3], (1+t^2)*eta(I*t,1)^8)), Iv));
print("");
print("4. THE EQUIANHARMONIC CM VALUE");
E = ellinit([0,0,0,0,-1/4]);   \\ y^2 = x^3 - 1/4  <=>  g2 = 0, g3 = 1
print("  omega = G(1/3)^3/(2 Pi)          rel.err = ", d(E.omega[1], w));
print("  omega2/omega = rho^-1            rel.err = ", d(E.omega[2]/E.omega[1], rho^-1));
z0 = E.omega[2]/3;  v0 = ellwp(E, z0, 1);      \\ z0 = omega rho^-1 / 3
z1 = E.omega[1]/3;  v1 = ellwp(E, z1, 1);      \\ rho z0 = omega/3
print("  wp (z0)      = ", v0[1]);
print("  wp'(z0)      = ", v0[2]);
print("  wp (omega/3) = ", v1[1]);
print("  wp'(omega/3) = ", v1[2]);
print("  (15) wp(z0) = rho^2              rel.err = ", d(v0[1], rho^2));
print("  (15) wp'(z0) = sqrt3             rel.err = ", d(real(v0[2]), sqrt(3)));
print("       wp'(omega/3) = -sqrt3       rel.err = ", d(real(v1[2]), -sqrt(3)));
print("  so (1, sqrt3) IS on the curve and IS 3-torsion -- it sits at z = -omega/3,");
print("  not at z0.  psi_3(x) = 3x(x^3-1), and every root of x^3=1 gives y^2 = 3,");
print("  which is why naming the wrong one moves no downstream number.");
print("  wp'^2 - (4 wp^3 - 1) at z0       = ", abs(v0[2]^2 - (4*v0[1]^3 - 1)));
print("  wp' is REAL at z0,          |Im| = ", abs(imag(v0[2])));
print("");
print("5. FINISH");
print("  B(1/3,2/3) = 2 Pi/sqrt3          rel.err = ", d(G*gamma(2/3), 2*Pi/sqrt(3)));
print("  J / B(1/3,2/3) = 7F6             rel.err = ", d(Jv/(2*Pi/sqrt(3)), real(Hv)));
ran524 ms
realprecision = 57 significant digits (40 digits displayed)
7F6                  = 1.039508495778944852704069517516990246343
9 G(1/3)^9/(64 Pi^6) = 1.039508495778944852704069517516990246343
THE IDENTITY                       rel.err = 3.186183822264904554 E-58

1-2. THE MELLIN BRIDGE AND BAILEY'S REAL-LINE W-INTEGRAL
  Plancherel (1/2pi) int M(1/2+it)M(1/2-it) dt = J
                                   rel.err = 9.558551466794713662 E-58
  (5) T / Gamma(2/3)^4 = J         rel.err = 6.372367644529809108 E-58
  (6) T = G(1/3) G(2/3)^5 * 7F6    rel.err = 6.372367644529809108 E-58
      T                = 12.67865694069209400911976310976140784560
      G(1/3)G(2/3)^5 H = 12.67865694069209400911976310976140784560

3. CUBIC THETA
  (16) I = sqrt3 G(1/3)^9/(32 Pi^5) rel.err = 4.779275733397356831 E-58
  (14) 2Pi int_1^oo (1+t^2) eta(it)^8 dt = I
                                   rel.err = 3.186183822264904554 E-58

4. THE EQUIANHARMONIC CM VALUE
  omega = G(1/3)^3/(2 Pi)          rel.err = 0.E-57
  omega2/omega = rho^-1            rel.err = 0.E-57
  wp (z0)      = -0.5000000000000000000000000000000000000000 + 0.8660254037844386467637231707529361834714*I
  wp'(z0)      = 1.732050807568877293527446341505872366943 - 7.965459555662261385 E-58*I
  wp (omega/3) = 1.000000000000000000000000000000000000000
  wp'(omega/3) = -1.732050807568877293527446341505872366943
  (15) wp(z0) = rho^2              rel.err = 2.8719872860574887955 E-58
  (15) wp'(z0) = sqrt3             rel.err = 3.186183822264904554 E-58
       wp'(omega/3) = -sqrt3       rel.err = 3.186183822264904554 E-58
  so (1, sqrt3) IS on the curve and IS 3-torsion -- it sits at z = -omega/3,
  not at z0.  psi_3(x) = 3x(x^3-1), and every root of x^3=1 gives y^2 = 3,
  which is why naming the wrong one moves no downstream number.
  wp'^2 - (4 wp^3 - 1) at z0       = 5.309734641931228801 E-57
  wp' is REAL at z0,          |Im| = 7.965459555662261385 E-58

5. FINISH
  B(1/3,2/3) = 2 Pi/sqrt3          rel.err = 1.5930919111324522770 E-58
  J / B(1/3,2/3) = 7F6             rel.err = 6.372367644529809108 E-58
python
# The one caveat, measured rather than asserted: the 2F1 Barnes representation
# needs |arg(-z)| < pi, and x in (0,1) lies ON the cut.  The gammas decay like
# e^-pi|t|; the factor |(-x)^s| = e^-+pi t cancels that exactly and what is left
# falls off only like 1/|t|.  Divergent, either branch.  (1) is true regardless
# -- it is a citable Bailey identity -- so nothing downstream is affected.
from mpmath import mp, mpf, j, pi, gamma, exp, log, nstr, fabs
mp.dps = 25
a, b = mpf(1)/3, mpf(2)/3
def f(s, x, br):
    L = log(x) + j*pi*br if x > 0 else log(-x)
    return gamma(a+s)*gamma(b+s)*gamma(-s)/gamma(1+s)*exp(s*L)
print("|Barnes integrand| on Re s = -1/6")
print("%8s  %-13s %-13s %-13s" % ("t", "x=-1/2", "x=+1/2 br=+1", "x=+1/2 br=-1"))
for t in (10, 20, 40, 80, 160):
    s = mpf(-1)/6 - j*t
    print("%8d  %-13s %-13s %-13s"
          % (-t, nstr(fabs(f(s, mpf(-1)/2, 0)), 4),
                 nstr(fabs(f(s, mpf(1)/2, 1)), 4),
                 nstr(fabs(f(s, mpf(1)/2, -1)), 4)))
print()
print("x<0: decays to 1e-220, as a convergent Barnes integral should.")
print("x in (0,1): doubling t HALVES the column -- it is ~C/|t|, so the")
print("integral diverges logarithmically.  One direction or the other,")
print("for either branch of (-x)^s.")
ran83 ms
|Barnes integrand| on Re s = -1/6
       t  x=-1/2        x=+1/2 br=+1  x=+1/2 br=-1 
     -10  1.602e-14     0.7054        3.639e-28    
     -20  1.819e-28     0.3527        9.382e-56    
     -40  4.691e-56     0.1763        1.248e-110   
     -80  6.24e-111     0.08816       4.416e-220   
    -160  2.208e-220    0.04408       1.106e-438   

x<0: decays to 1e-220, as a convergent Barnes integral should.
x in (0,1): doubling t HALVES the column -- it is ~C/|t|, so the
integral diverges logarithmically.  One direction or the other,
for either branch of (-x)^s.