SF [2159499] Full bigfloat precision for Gamma after the second call
Subject: SF [2159499] Full bigfloat precision for Gamma after the second call
From: Raymond Toy
Date: Mon, 13 Oct 2008 12:09:52 -0400
Dieter Kaiser wrote:
> As reported on Sourceforge.net I have studied the subtle bug that we do not get
> the full precision of the Gamma function in the first call. I can reproduce the
> bug with GCL 2.6.8 and CLISP 2.44.
>
> I have studied the effect in more detail and detected that we have the
> problem only for a negative argument to bffac. That is true for
> gamma(0.5b0).
>
> The reflection formula calculates:
>
> bfloat(%pi*z/sin(%pi*z))/bffac(-z,fpprec)
>
I looked into this a bit. It seems to me that the problem is actually
in the value of pi. I put a break point at *fpsin, and this is what I
see, when we change the precision from 64 to 128. The arg to sin is:
((MTIMES SIMP)
((BIGFLOAT SIMP 428)
346583711765101857447301773017885462929554634421977071896309947576827663475703202879996800763017447262173901370175446478621769728
0)
$%PI))
The bigfloat there works out to be exactly 0.5b-1, as expected:
(displa '((BIGFLOAT SIMP 428)
346583711765101857447301773017885462929554634421977071896309947576827663475703202879996800763017447262173901370175446478621769728
0))
5.0b-1
But look at the value of bigfloat%pi
((BIGFLOAT SIMP 428)
544412421367563193456719682845547953364731116819812684167733550941339415157339997517605497896021017270535788777222242612081295775
2)
0] (displa '((BIGFLOAT SIMP 428)
544412421367563193456719682845547953364731116819812684167733550941339415157339997517605497896021017270535788777222242612081295775
2))
3.1415926535897932384626433832795028841971693993751058209749445923078164062862\
089986280348256812002485034480173261900300090541701b0
Compared to the actual value of %pi:
(%i5) bfloat(%pi);
(%o5)
3.1415926535897932384626433832795028841971693993751058209749445923078164\
062862089986280348253421170679821480865132823066470938446b0
I think the problem is that when fpprec is set, fpprec1 is called. But
this only updates the bfloat versions of 1, 0, +/- 1/2. Perhaps we need
to update the bigfloat values of %pi at least, and probably also %e and
%gamma.
Ray