Why the results are different using 'hypergeom' function and 'pochhammer' function to calculate the Gauss hypergeometric function

First, I calculated the the Gauss hypergeometric function using hypergeometric function. For example,

>> hypergeom([-30,30.5],1,0.64)

ans =

      0.0379

And

>> F(-30,30.5,1,0.64)

ans =

    -86.1497

The 'F' function I wrote using pochhammer function is as follows.

function f = F(n,a,c,x)

k = 0:1:abs(n);
vector = pochhammer(n,k).*pochhammer(a,k).*(x.^k)./pochhammer(c,k)./factorial(k);
f = sum(vector);

end

However, they are same for some other parameters. For example,

>> F(-10,10.5,1,0.64)

ans =

      0.1592

>> hypergeom([-10,10.5],1,0.64)

ans =

      0.1592

Both the results are 0.1592. Why is that?

https://au.mathworks.com/help/symbolic/hypergeom.html#bt1nxzg-5

http://au.mathworks.com/help/symbolic/pochhammer.html

Above are the links of the two functions.

Answer 1

[David Goodmanson]

Hi huanyu,

The answer is pretty simple: lack of precision. I modified F in the code below to output the vector that you sum to get the answer.

When there is a wide range of values in a vector, the standard Matlab display scales by the largest elements and effectively does not display the small ones. There are better ways to do an improved display but here I settled for just reading out all the values separately.

format long
format compact
R = 10;
[f vector] = F(-R,R+1/2,1,0.64);
for k = 1:length(vector)
  disp(vector(k))
end

function [f vector] = F(n,a,c,x)
k = 0:1:abs(n);
vector = pochhammer(n,k).*pochhammer(a,k).*(x.^k)./pochhammer(c,k)./factorial(k);
f = sum(vector);
end

       1
    -67.200000000000003
     1.112832000000000e+03
    -7.913472000000001e+03
     2.991292416000000e+04
    -6.662206468915202e+04
     9.179040023838721e+04
    -7.912707155243828e+04
     4.154171256503009e+04
    -1.214454016715942e+04
     1.515638612861496e+03

So far so good. Now try it with R = 30:

     1
    -5.856000000000000e+02
     8.559129600000000e+04
    -5.538707865600000e+06
     2.003904505774080e+08
    -4.601605994699137e+09
     7.260311680525304e+10
    -8.306981919530012e+11
     7.164771905594635e+12
    -4.794913030351531e+13
     2.545523429553021e+14
    -1.090577971801890e+15
     3.821869914514622e+15
    -1.107211544465301e+16
     2.673576937574172e+16
    -5.414646839616882e+16
     9.238741170096304e+16
    -1.331913190695130e+17
     1.624605225193565e+17
    -1.676268571413018e+17
     1.460365179415021e+17
    -1.070272167770827e+17
     6.559618509246656e+16
    -3.333129405076562e+16
     1.386952180223525e+16
    -4.644181636434868e+15
     1.220128193240877e+15
    -2.420841452268320e+14
     3.408940004214570e+13
    -3.035213169031953e+12
     1.284232416408186e+11

You are adding and subtracting large numbers on the order of 1e17 with around 16 digits of precision, trying to come up with a small answer. So it can't be done.

[Walter Roberson]

You could switch to the symbolic pochhammer for higher precision.

Alternately you could use https://www.mathworks.com/matlabcentral/fileexchange/26800-xsum for the Kahan Compensated Sum

[Huanyu Zuo]

Thanks David for your answer and Walter for your advice.

I've got the precise result after revising the 'F' function by using symbolic function.

function f = F(n,a,c,x)

syms y(n1,a1,c1,x1,k)
y(n1,a1,c1,x1,k) = pochhammer(n1,k).*pochhammer(a1,k).*(x1.^k)./pochhammer(c1,k)./factorial(k);
S = sum(y(n,a,c,x,0:1:abs(n)));
f = vpa(S,16);

end

The result is as follows

>> F(-30,30.5,1,0.64)
ans =
0.03787177451622411

Thanks again!