Weird behavior of `HypergeometricPFQ`

Question

Consider the following function:

fun = HypergeometricPFQ[{1, x1-y1, 1-x1-y1-2u}, {x2-y1, 3-x2-y1-2u}, 1];

Let us try and evaluate this function on the following specific values of the arguments:

sub = {x1 -> 3, x2 -> 2, y1 -> 7, y2 -> 4, u -> 51/10};
fun/.sub

43/11

Alternatively, we can functionally expand fun first to obtain

funEx = FunctionExpand[fun, Element[x1 | y1 | x2 | y2, Integers] && Element[u, Rationals] && y1 >= x1 && y1 >= x2 && y2 >= x2 && y1 >= y2 && x1 >= x2]

-((2 (-1 + x2 - y1))/((-1 - x1 + x2) (-2 + 2 u + x1 + x2))) + ( 2 u (-1 + x2 - y1))/((-1 - x1 + x2) (-2 + 2 u + x1 + x2)) + ( x2 (-1 + x2 - y1))/((-1 - x1 + x2) (-2 + 2 u + x1 + x2)) + ((-1 + x2 - y1) y1)/((-1 - x1 + x2) (-2 + 2 u + x1 + x2)) - ( Gamma[3 - 2 u - x2 - y1] Gamma[ x2 - y1])/((1 + x1 - x2) (-2 + 2 u + x1 + x2) Gamma[ 1 - 2 u - x1 - y1] Gamma[x1 - y1])

First of all we see that even though FunctionExpand was explicitly told that y1>=x1 and y1>=x2 and x1,x2,y1 are integers, it still produced the ill behaved ratio Gamma[x2 - y1]/Gamma[x1 - y1], such that trying to evaluate funEx/.sub leads to the error:

Indeterminate

What is worse, even if we try to regularize by setting y1->7+e and expanding around small e, we find a discrepancy with the original evaluation above:

Series[funEx /. y1 -> 7 + e /. sub, {e, 0, 0}] // FullSimplify// Normal

-(5203/125)

Which is clearly different from 43/11. Why does this function behave this way? I need to simplify the function analytically and I am only evaluating it to check if the functional expansion was correct. Clearly it is not correct and I am at a loss on how to proceed. Any suggestions on how to obtain the correct simplification?

EDIT:

J. M. suggested to use the Thomae transformation. In the above case this leads to

subv={a -> 1, b -> x1 - y1, c -> 1 - x1 - y1 - 2 u, e -> x2 - y1, f -> 3 - x2 - y1 - 2 u};
s = e + f - a - b - c;
fun2=(Gamma[e] Gamma[f] Gamma[s])/(Gamma[a] Gamma[s + b] Gamma[s + c])HypergeometricPFQ[{s, e - a, f - a}, {s + b, s + c}, 1]/.subv

Now, trying to functionally expand this as before, gives

funEx2=FunctionExpand[fun2, Element[x1 | y1 | x2 | y2, Integers] && Element[u, Rationals] && y1 >= x1 && y1 >= x2 && y2 >= x2 && y1 >= y2 && x1 >= x2] //FullSimplify

-(((-1 + x2 - y1) (-2 + 2 u + x2 + y1) + (Gamma[3 - 2 u - x2 - y1] Gamma[x2 - y1])/(Gamma[1 - 2 u - x1 - y1] Gamma[x1 - y1]))/((1 + x1 - x2) (-2 + 2 u + x1 + x2)))

which still produces the same ill behaved ratio of gamma functions.

Resolution:

Just trying to evaluate the function fun a few times on some different y1 and x1 values like

fun/. y1 -> 3 /. x1 -> 2 // Simplify

((-4 + x2) (1 + 2 u + x2))/((-3 + x2) (2 u + x2))

One can infer that the whole hypergeometric function actually in general reduces to just

((-1 + x2 - y1) (-2 + 2 u + x2 + y1))/((-1 + x2 - x1) (-2 + 2 u + x1 + x2))

This is basically what both the functional expansions above produced, but with all the gamma functions dropped. Even though I am almost certain now that this is the correct expansion, this way to do things is highly unsatisfying. Mathematica should do better than this.

Answer 1

Just trying to evaluate the function fun a few times on some different y1 and x1 values like

fun/. y1 -> 3 /. x1 -> 2 // Simplify

((-4 + x2) (1 + 2 u + x2))/((-3 + x2) (2 u + x2))

One can infer that the whole hypergeometric function actually in general reduces to just

((-1 + x2 - y1) (-2 + 2 u + x2 + y1))/((-1 + x2 - x1) (-2 + 2 u + x1 + x2))

This is basically what both the functional expansions in the question above produced, but with all the gamma functions dropped. Even though I am almost certain now that this is the correct expansion, this way to do things is highly unsatisfying. Mathematica should do better than this.

EDIT

With a little bit of hindsight, it is possible to interpret mathematica output a bit better. As already mentioned, the first part which is independent of the gamma functions:

((-1 + x2 - y1) (-2 + 2 u + x2 + y1))/((1 + x1 - x2) (2 - 2 u - x1 - x2))

is the result we are after. But why do we get some extra terms with gamma functions - is something going wrong or does it have an interpretation? Well, it turns out that if we regularize the gamma functions, we get the extra contribution:

(-1)^(x1 - x2) (Gamma[3 - 2 u - x2 - y1] Gamma[y1 - x1 + 1] )/( Gamma[1 - 2 u - x1 - y1] Gamma[ y1 - x2 + 1] (1 + x1 - x2) (2 - 2 u - x1 - x2))

This contribution makes up for the difference in results when the arguments are non-integers. Having a negative integer for one of the first arguments in the hypergeometric function makes the sum terminate, which gives back just the first part of the result (and which we are interested in here). If this one argument is even slightly different from integer, the summation continues until infinity - and picks up exactly the second part of the result involving the gamma functions. From this perspective it is also clear why the gamma functions needed to be regularized - they only arise in the case when the argument actually is not integer so that the sum does not terminate.