Re: zetakinit() puzzle

[Karim BELABAS <Karim.Belabas@math.u-psud.fr>]

On Wed, 26 Feb 2003, Bill Allombert wrote:
> On Tue, Feb 25, 2003 at 11:55:02PM +0100, Karim BELABAS wrote:
>> On Mon, 24 Feb 2003, Igor Schein wrote:
>>> Since we're talking about precision loss here, I'll mention a non-related
>>> case which I noticed some time ago:
>>>
>>> ? precision(sin(precision(2^2^7+.,38)))
>>> 9
>>>
>>> I can't judge whether it's a bug or not because I first need to understand
>>> something more basic ( see my message 2407 posted a few days ago).
>>
>> Not a bug. I don't see how you can possibly reduce this mod 2Pi without
>> suffering catastrophic cancellation.
>>
>> You actually escape a precision error by a narrow margin.
>
> What is a bug in a sense is the following:
>
> ? sin(2^22)
> %1 = 0.9751293949417070368170374003
> ? sin(2^22*1.)
> %2 = 0.9751293949417070368170374003
> ? sin(2^22+.)
> %3 = 0.9751293949417070368170274962
> The correct result being the last one.

Yes, this is quite a different situation. When the input is exact,
counter-measures can be taken.

> sin() could be smart enough to reduce mod 2Pi correctly when the
> input is exact, using sizedigit(x\3)+prec digits of Pi instead of prec.
>
> The following patch fix that for sin, cos, tan, cotan and sincos (used by
> exp(I*x)).  This is not perfect, since sometimes the result will have more
> than default(realprecision) words of precision.

This should be trivial to fix: results are t_REAL. Just affrr them to
cgetr(prec) and remove the gerepiles [ which are slightly slower anyway ].

To be extra careful, you could check for an exact 0 input and directly return
realzero(prec) / realun(prec) instead. But it shouldn't make any difference.

Of course, we have the same kind of problems with exact t_COMPLEX and
t_QUADs...

    Karim.
-- 
Karim Belabas                    Tel: (+33) (0)1 69 15 57 48
Dép. de Mathématiques, Bât. 425  Fax: (+33) (0)1 69 15 60 19
Université Paris-Sud             Email: Karim.Belabas@math.u-psud.fr
F-91405 Orsay (France)           http://www.math.u-psud.fr/~belabas/
--
PARI/GP Home Page: http://www.parigp-home.de/