Re: zetakinit() puzzle

[Bill Allombert <allomber@math.u-bordeaux.fr>]

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);
+    
     case t_INTMOD: case t_PADIC: err(typeer,"gcos");
 
     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);
 
     case t_INTMOD: case t_PADIC: err(typeer,"gsin");
 
@@ -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);
+
     case t_INTMOD: case t_PADIC: err(typeer,"gtan");
     
     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);
 
     case t_INTMOD: case t_PADIC: err(typeer,"gcotan");
 



Cheers,
Bill.