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.
I have no idea what is wrong.
Dieter Kaiser