Bill Allombert on Wed, 26 Feb 2003 17:44:25 +0100

 Re: zetakinit() puzzle

```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.

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.

Index: src/basemath/trans1.c
===================================================================
RCS file: /home/megrez/cvsroot/pari/src/basemath/trans1.c,v
retrieving revision 1.81
diff -u -r1.81 trans1.c
--- src/basemath/trans1.c	2003/01/15 20:46:02	1.81
+++ src/basemath/trans1.c	2003/02/26 16:41:12
@@ -1673,6 +1673,29 @@
/**                                                                **/
/********************************************************************/

+/*Transform an exact number to a real with sufficient accuracy
+ *to avoid precision loss in modulo Pi reduction*/
+
+static GEN
+mpsc_exact(GEN x, long prec)
+{
+  long t=typ(x);
+  GEN  p1;
+  long pr=prec, d;
+  switch(t)
+  {
+  case t_INT:
+    pr += lgefint(x)-2;
+    break;
+  default:
+    d=lgefint(x[1])-lgefint(x[2])+1;
+    if (d>0)
+     pr += d;
+  }
+  p1=cgetr(pr); gaffect(x,p1);
+  return p1;
+}
+
/* Reduce x0 mod Pi/2 to x in [-Pi/4, Pi/4]. Return cos(x)-1 */
static GEN
mpsc1(GEN x0, long *ptmod8)
@@ -1827,6 +1850,11 @@
gerepilemanyvec(av,tetpil,y+1,2);
return y;

+    case t_INT: case t_FRAC: case t_FRACN:
+      av=avma; p1=mpsc_exact(x,prec); tetpil=avma;
+      p1=mpcos(p1);
+      return gerepile(av,tetpil,p1);
+

default:
@@ -1900,6 +1928,11 @@
y[2]=lmul(p1,v);
gerepilemanyvec(av,tetpil,y+1,2);
return y;
+
+    case t_INT: case t_FRAC: case t_FRACN:
+      av=avma; p1=mpsc_exact(x,prec); tetpil=avma;
+      p1=mpsin(p1);
+      return gerepile(av,tetpil,p1);

@@ -1979,7 +2012,7 @@
switch(typ(x))
{
case t_INT: case t_FRAC: case t_FRACN:
-      av=avma; p1=cgetr(prec); gaffect(x,p1); tetpil=avma;
+      av=avma; p1=mpsc_exact(x,prec); tetpil=avma;
mpsincos(p1,s,c); gptr[0]=s; gptr[1]=c;
gerepilemanysp(av,tetpil,gptr,2);
return;
@@ -2098,6 +2131,11 @@
av = avma; gsincos(x,&s,&c,prec);
return gerepileupto(av, gdiv(s,c));

+    case t_INT: case t_FRAC: case t_FRACN:
+      av=avma; y=mpsc_exact(x,prec);
+      y=mptan(y);
+      return gerepileupto(av,y);
+

default:
@@ -2144,6 +2182,11 @@
case t_COMPLEX:
av = avma; gsincos(x,&s,&c,prec);
return gerepileupto(av, gdiv(c,s));
+
+    case t_INT: case t_FRAC: case t_FRACN:
+      av=avma; y=mpsc_exact(x,prec);
+      y=mpcotan(y);
+      return gerepileupto(av,y);