Dieter Kaiser wrote:
> Hello Ray,
>
> I had a look at the problem.
>
> The bug vanishes if I remove the 4 extra digits in the call to bffac in the
> routine simpgamma.
>
> Here again is the example with the 4 extra digits and some debugging info:
>
> That is the example with the wrong result. Simpgamma calculats two times the
> Gamma function in bigfloat precision. The results have extra digits but fpprec
> has after the call of bffac again the value 16.
>
> (%i2) bfloat(gamma(3/4)/gamma(1/4));
>
> in GAMMA bigfloat-numerical for
> ((%GAMMA) ((BIGFLOAT SIMP 56) 36028797018963968 -1))
> : result = ((BIGFLOAT SIMP 69) 535045584704600947817 2)
> : fpprec = 16
>
> in GAMMA bigfloat-numerical for
> ((%GAMMA) ((BIGFLOAT SIMP 56) 54043195528445952 0))
> : result = ((BIGFLOAT SIMP 69) 361679172704387465868 1)
> : fpprec = 16
>
> (%o2) 4.125843750410604b-5
>
> Here the correct calculation. At first the numerical values of the Gamma
> function are calculated in bigfloat precision and then divided. The results of
> simpgamma looks identically.
>
> (%i3) bfloat(gamma(3/4))/bfloat(gamma(1/4));
>
> in GAMMA bigfloat-numerical for
> ((%GAMMA) ((BIGFLOAT SIMP 56) 54043195528445952 0))
> : result = ((BIGFLOAT SIMP 69) 361679172704387465868 1)
> : fpprec = 16
>
> in GAMMA bigfloat-numerical for
> ((%GAMMA) ((BIGFLOAT SIMP 56) 36028797018963968 -1))
> : result = ((BIGFLOAT SIMP 69) 535045584704600947821 2)
> : fpprec = 16
>
> (%o3) 3.379891200336424b-1
>
> I have done the calculation without extra digits to bffac. Without extra digits
> the results of both calculations are correct.
>
>
Yes, I've come more or less to the same conclusion. Also, as an
experiment, I added a call to fpround before simpgamma returns the
bfloat result. This fixes the issue.
Based on this, I think something somewhere is getting confused by the
fact that gamma returns a bfloat with 69 bits instead of 56. I haven't
tracked this down yet. But the path taken is somehow different between
bfloat(gamma(3/4)/gamma(1/4)) and bfloat(gamma(3/4))/bfloat(gamma(1/4)),
because, as you show, gamma returns the same bfloat in both cases.
Ray