SF [2159499] Full bigfloat precision for Gamma after the second call

>>>>> "Stavros" == Stavros Macrakis <macrakis at alum.mit.edu> writes:

    Stavros> On Tue, Oct 14, 2008 at 2:03 PM, Raymond Toy <raymond.toy at ericsson.com>wrote:
    >> Ok, ok. :-)
    >> 

    Stavros> OK, a Maxima-level demo like this is convincing.  Thanks.

    Stavros> I think fppi is correct.  It's the computation in fppi1 that is
    >> incorrect because it uses $fpprec instead of fpprec.
    >> 

    Stavros> Try your test case again, with an addition:

    Stavros> zot(z,fpprec):= bfloat(%pi*z/sin(%pi*z))$
    Stavros> fpprec:64$
    Stavros> zot(bfloat(1/2),fpprec)$
    Stavros> fpprec:128$
    Stavros> bfloat(%pi);        <<<<<<< fppi returns without calling fppi1 for a new
    Stavros> value of %pi
    Stavros> bfloat(sin(2))      <<<<<<< add this to test; fppi calls fppi1 and gets
    Stavros> correct value
    Stavros> bfloat(%pi);        <<<<<<< now gives correct result

    Stavros> So it looks as though the problem is not the pi calculation itself (i.e.
    Stavros> fppi1), but something in the logic of fppi or fpprec1.

Warning:  Long explanation follows.

Don't you think it's a problem that the cached value in bigfloat%pi is
wrong?

Anyway, here is why I think the problem is fppi1, and, specifically,
fprt18231_.  (The tracing is cmucl-specific, but illustrates the
problem.)

zot(z,fpprec):= bfloat(%pi*z/sin(%pi*z))$
fpprect:64$

/* Trace fppi1 and fprt18231_  Print out fpprec and $fpprec on entry */
:lisp (trace fppi1 :print (list fpprec) fprt18231_ :print (list $fpprec))

/* Note that $fpprec = 64 means fpprec is 215 */
(%i3) zot(bfloat(1/2),fpprec);

  0: (FPPI1)
  0: (LIST FPPREC) = (215)
    1: (FPRT18231_)
    1: (LIST $FPPREC) = (64)
    1: FPRT18231_ returned
         (33502987467293214945479741689670280507273782561274100088049553032 26)
  0: FPPI1 returned
       ((BIGFLOAT SIMP 215)
        41356040229830605799206175877058042562270358908835127062817736159 2)
  0: (FPPI1)
  0: (LIST FPPREC) = (224)
    1: (FPRT18231_)
    1: (LIST $FPPREC) = (64)

*** This is a problem but it's harder to see.  Look below.
    1: FPRT18231_ returned
         (17153529583254126052085627745111183619724176671372339245081371152496
          26)
  0: FPPI1 returned
       ((BIGFLOAT SIMP 224)
        21174292597673270169193562049053717791882423761323585056162680913575 2)
  0: (FPPI1)
  0: (LIST FPPREC) = (438)
    1: (FPRT18231_)
    1: (LIST $FPPREC) = (64)
    1: FPRT18231_ returned
         (451619377654220682466952889210173098222272334731830641620256195774520016906666251050762585786556344451681167482041915198196893760688
          26)

*** At this point maxima wants 438 bits of pi.  But fprt18231_ sees
    $fpprec = 64 and therefore only computes 64 digits or about 215
    bits.  We really want 438 bits.  I think it was a misunderstanding
    about fpprec and $fpprec.  fpprec changes throughout the code and
    is dynamically changed to produce the desired precision.  I think
    sin is especially likely to increase fpprec so that it can do
    range reduction accurately.

  0: FPPI1 returned
       ((BIGFLOAT SIMP 438)
        557478319480384710099680955233841104245484663623488188587759156163931561121116157458028029845525521685028647707875576434771246873456
        2)

This value is incorrect, and is stored in bigfloat%pi.

In your particular example, bfloat(%pi), just reuses the incorrect
value of bigfloat%pi.  But fpprec is now 428, so fppi rounds the
stored value of bigfloat%pi back in to bigfloat%pi.  This is why the
answer is wrong.

Then bfloat(sin(2)) does something different.  It wants a 438 bit
version of pi, but the current stored version only has 428 bits.  A
new value is needed:

  0: (FPPI1)
  0: (LIST FPPREC) = (438)
    1: (FPRT18231_)
    1: (LIST $FPPREC) = (128)

*** Note that $ffprec is 128 now, which is almost the same as fpprec =
    438.  We lose some bits, but not enough for anyone to have noticed
    before.  And since fprt18231 computes 428 bits or so, the value of
    pi is close to the true value.  This value 438-bit value is saved
    in bigfloat%pi.

The final call to bfloat(%pi) gives the correct result because we only
take 428 of the 438 bits, which matches the precision (roughly) from
fprt18231_, and we get the correct value.

Does that make sense?

Ray