Karim BELABAS on Wed, 26 Feb 2003 18:01:41 +0100 (MET) |
[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]
Re: zetakinit() puzzle |
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/