# Expansion of a hypergeometric function takes ages with Mathematica 9 and 10 (regression?)

Mathematica 8 (Linux version) can evaluate

AbsoluteTiming[Series[Hypergeometric2F1[1, 1 - eps/2, 3 - eps, 1/2], {eps, 0, 0}]]


in no time. On one of the university machines the result is

{0.175840, (4 - 4 Log[2]) + O[eps]}


For some reason both Mathematica 9 and 10 seem to have a lot of troubles with this expansion. On Mathematica 10 I had to stop the evaluation after 4 minutes without getting any result. Mathematica 9 tried harder and run out of memory after a couple of minutes:

No more memory available.
Mathematica kernel has shut down.
Try quitting other applications and then retry.


Unless I'm missing something, this looks pretty much like a regression as compared to the version 8.

Can someone confirm this?

• I confirm it seems to take a long time in mathematica 10.0.2 Mar 29, 2015 at 18:42
• It is nearly instantaneous if you use the equivalent Hypergeometric2F1[1, 1 - eps/2, 3 - eps, 1/2] /. eps -> 0 or plot with Plot[ Hypergeometric2F1[1, 1 - eps/2, 3 - eps, 1/2], {eps, -2, 2}] Mar 29, 2015 at 18:48
• {eps, 0, 0} seems to be a Typo.
– user9660
Mar 29, 2015 at 19:07
• @BobHanlon Unfortunately, this kind of workaround doesn't really help me. The code I provided it just the smallest working example to reproduce the problem. In reality (1-loop integrals in QFT), I have expressions where such functions are multiplied by other functions that also depend on eps and are divergent if one just puts eps to 0 (e.g Gamma[eps] etc.). This is why I must use Series and it is very unfortunate that for this particular case it fails in MMA 10, while in the version 8 everything is fine.
– vsht
Mar 29, 2015 at 19:29
• @Lou If you like, you can replace it by {eps,0,1} or {eps,0,2}. This still takes ages, unfortunately.
– vsht
Mar 29, 2015 at 19:39

Your reply to Jens's answer suggests that you do not want to manually replace variables. This function just automates the method in Jens answer...

cme[hgm_, symb_: c] := Module[
{ places, symbs },
places = Position[hgm, eps];
symbs = Array[symb, Length[places]];
FullSimplify[Normal[Series[
ReplacePart[hgm, (Rule @@ # &) /@ Transpose[{places, symbs}]],
Sequence @@ ({#, 0, 0} & /@ symbs)
]]]
];


On the machine I'm sitting at, it takes 3.4 ms to do your example in Mathematica 10.0.2.0, Linux x64.

AbsoluteTiming[
cme[Hypergeometric2F1[1, 1 - eps/2, 3 - eps, 1/2]]
]


yields

{0.003400, 4 - 2 Log[4]}


Edit 13 May 2015:

cme[hgm_, symb_: eps, ord_: 0] := Module[
{ places, symbs, c },
places = Position[hgm, symb];
symbs = Array[c, Length[places]];
Series[
FullSimplify[Normal[Series[
ReplacePart[hgm, (Rule @@ # &) /@ Transpose[{places, symbs}]],
Sequence @@ ({#, 0, ord} & /@ symbs)] /. c[_] -> symb
]]
, {symb, 0, ord}]
]


(If one does not care about having a series as the result, drop "Series[" and ", {symb, 0, ord}]" each appearing on their own lines.) On the same machine as the prior timing data,

AbsoluteTiming[
cme[Hypergeometric2F1[1, 1 - eps/2, 3 - eps, 1/2]]
]


yields

{0.003157, 4 - 2 Log[4]}


and

AbsoluteTiming[
cme[Hypergeometric2F1[1, 1 - eps/2, 3 - eps, 1/2], eps, 3]
]


yields

{0.064078,
(4-2 Log[4]) +
(2-\[Pi]^2/3+Log[4]-1/2 (Hypergeometric2F1^(0,1,0,0))[1,1,3,1/2]) eps +
1/24 (48+4 \[Pi]^2-72 Zeta[3]+12 (Hypergeometric2F1^(0,1,1,0))[1,1,3,1/2]+3 (Hypergeometric2F1^(0,2,0,0))[1,1,3,1/2]) eps^2 +
1/720 (1440-28 \[Pi]^4+1080 Zeta[3]-180 (Hypergeometric2F1^(0,1,2,0))[1,1,3,1/2]-90 (Hypergeometric2F1^(0,2,1,0))[1,1,3,1/2]-15 (Hypergeometric2F1^(0,3,0,0))[1,1,3,1/2]) eps^3 +
O[eps]^4
}


We could attempt the same sort of expansion and collapse for the derivatives of the hypergeometric functions, but this temds to just produce higher order derivatives and a lot of large (eventually) cancelling terms. Instead, I refer you to the HypExp package, discussed at the question "Expanding derivatives of hypergeometric functions".

HypExp will also do epsilon expansions of exactly the form you are interested. Invocation is via

HypExp[Hypergeometric2F1[...], eps, order]


as described in section 5.1 of the HypExp paper for the package.

• It's indeed a nice workaround for people who are stuck with version 10. As I said in the comment to Jens' answer, for now I'll better stick with Mathematica 8 (luckily my university provides this possibility).
– vsht
Mar 30, 2015 at 12:34
• To make a general answer, this shouldn't have the order of expansion hard-coded. As it is, the result could also be obtained with Limit. I understood the comment to mean that he needs higher orders, and that was the whole point of my answer.
– Jens
Mar 30, 2015 at 16:45

Expand@Normal@Series[Hypergeometric2F1[1, 1 - eps/2, 3 - eps, x], {eps,0,0}]/.x->(1/2)


giving

4 - 4 Log[2]

• I wonder why this wasn't the accepted answer - it seems the most straightforward (+1). Maybe the reason this works is the same as why my answer works: the series expansions for the hypergeometric function for undefined arguments always seem to work best.
– Jens
Mar 30, 2015 at 16:17
• Oh, actually, this doesn't really work for higher orders, so it's as limited as the approach of taking the Limit. Yes, it appears to spit out a result if you do the series with 1 instead of 0, but it's in terms of derivatives of the hypergeometric function, which when I try to evaluate numerically doesn't return a result. So this doesn't get us back to the behavior in version 8, as my answer does.
– Jens
Mar 30, 2015 at 16:42

This does indeed work flawlessly with Mathematica version 8. But in version 10, I had to resort to the following workaround:

AbsoluteTiming[
Series[Normal[
Series[Hypergeometric2F1[1, 1 - eps1, 3 - eps, 1/2], {eps1, 0,
1}, {eps, 0, 1}]] /. eps1 -> (eps/2), {eps, 0, 0}]]

(*
==> {0.016225, SeriesData[eps, 0, {2 (EulerGamma + PolyGamma[0,
Rational[3, 2]])}, 0, 1, 1]}
*)


To check that this yields the correct result, I converted to numbers and compared to the analogous result in version 8, getting complete agreement:

N[%]

(*
==> {0.016225, SeriesData[
eps, 0., {1.2274112777602189}, 0, 1, 1]}
*)


This can now be continued to higher order as well.

AbsoluteTiming[
Series[Normal[
Series[Hypergeometric2F1[1, 1 - eps1, 3 - eps, 1/2], {eps1, 0,
1}, {eps, 0, 1}]] /. eps1 -> (eps/2), {eps, 0, 1}]]

(*
==> {0.013317, SeriesData[eps, 0, {2 (EulerGamma + PolyGamma[0,
Rational[3, 2]]), 4 - EulerGamma + Rational[-1, 3] Pi^2 - PolyGamma[0,
Rational[3, 2]] + Rational[-1, 2] Derivative[0, 1, 0, 0][
Hypergeometric2F1][1, 1, 3,
Rational[1, 2]]}, 0, 2, 1]}
*)

N[%]

(*
==> {0.013317, SeriesData[
eps, 0., {1.2274112777602189, -0.03682639753846961}, 0, 2, 1]}
*)


Again, the result agrees numerically with what version 8 yielded. The comparison was easiest numerically because the exact form of the results differs between the versions.

To explain what I did: I simply introduced two different expansion variables, one for each argument in the hypergeometric function, and expanded first in both those variables individually. This works without problems. Afterwards, I replace one of the two (eps1) by eps/2 as it was desired in the original expansion. Then I just have to make sure that we throw out all terms that aren't consistent with the desired order of the expansion. This is done by adding a second Series command. In that expansion, we are only dealing with a polynomial in eps, so there will never be any difficulties.

• Nice, but not very practical for my case, where I get 13 different Hypergeometric2F1 that need to be expanded in eps (that only for one particular loop integral). I guess I'll just stick to Mathematica 8 for now. So far, to me it seems to the version with the least number of bugs.
– vsht
Mar 29, 2015 at 21:51
• @vsht Yes, that's why I also kept version 8... but in principle it would be possible to automate the method I used here. You "just" have to wrap it in a rule that picks out all hypergeometric functions. The implementation will of course depend on how different the various terms look. So I agree, it's in many cases easier to not upgrade from version 8 if you want to avoid having to fix stuff that will eventually change again anyway.
– Jens
Mar 29, 2015 at 22:00

There is an approach that works to all orders—and I expect would be useful for other expressions arising in 1-loop integrals in QFT.

Using a Hypergeometric2F1 integral definition, e.g.

Gamma[c]/(Gamma[b] Gamma[c - b]) *
Integrate[(t^(b - 1)*(1 - t)^(c - b - 1))/(1 - z*t)^a, {t, 0, 1}]


the sum you are examining,

Hypergeometric2F1[1, 1 - eps, 3 - eps, 1/2]


can be computed from the integrand of the integral definition,

int = FullSimplify[
(Gamma[c]/(Gamma[b]*Gamma[c - b]))*
((t^(b - 1)*(1 - t)^(-b + c - 1))/(1 - t*z)^a) /.
{a -> 1, b -> 1 - eps/2, c -> 3 - eps, z -> 1/2}]


and then using (ordinary) series expansion in $\epsilon$ to any desired order,

intseries = ExpandAll /@ Apart /@ (int + O[eps]^2)


followed by integration wrt $t$ over $0<t<1$,

(FullSimplify[Integrate[#, {t, 0, 1}]] & ) /@ intseries


yields closed-form coefficients, rather than as derivatives of Hypergeometric2F1 functions. This is what Series should yield automatically.