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

Here a one-liner twitterable solution:

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

giving

4 - 4 Log[2]

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.

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:

To address Jens comment about hard-coded order, and also to produce a series instead of just the leading order terms:

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.