Dieter Kaiser wrote:
> 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.
>
>
Some additional information. I tracked down the issue to exptbigfloat
and invertbigfloat. At one point, exptbigfloat is called to compute
gamma(1/4)^(-1), where gamma(1/4) has already been converted to a
bfloat. Then exptbigfloat calls invertbigfloat on gamma(1/4). This is
fine, I think. But when invertbigfloat calls fpquotient, there is a
problem. fpquotient is called:
(FPQUOTIENT (36028797018963968 1) (535045584704600947821 2))
The first argument is 1b0, using fpprec=56 bits. The second arg is
gamma(1/4), using 69 bits. But I'm pretty sure these fp functions
assume all args are supposed to have the same precision.
So, perhaps, the correct solution is to make simpgamma round the result
from bffac to fpprec bits before returning.
Ray