Re: zetakinit() puzzle

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

On Wed, Feb 26, 2003 at 06:01:41PM +0100, Karim BELABAS wrote:
> On Wed, 26 Feb 2003, Bill Allombert wrote:
> > 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...

I suppose gsincos take care of t_COMPLEX ? t_QUAD are more tricky to handle.

Here a new patch (that still need to change the bench).

Now, I wonder if mpsc_exact is an overkill. Should not
mpsc_exact(GEN x,long prec){return gadd(x,realzero(prec));}
be as good, and work for exact t_QUAD as well ?

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 17:54:48
@@ -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,12 @@
       gerepilemanyvec(av,tetpil,y+1,2);
       return y;
 
+    case t_INT: case t_FRAC: case t_FRACN:
+      p1=cgetr(prec); av=avma;
+      p2=mpsc_exact(x,prec); 
+      affrr(mpcos(p2),p1); avma=av;
+      return p1;
+    
     case t_INTMOD: case t_PADIC: err(typeer,"gcos");
 
     default:
@@ -1900,6 +1929,12 @@
       y[2]=lmul(p1,v);
       gerepilemanyvec(av,tetpil,y+1,2);
       return y;
+    
+    case t_INT: case t_FRAC: case t_FRACN:
+      p1=cgetr(prec); av=avma;
+      p2=mpsc_exact(x,prec);
+      affrr(mpsin(p2),p1); avma=av;
+      return p1;
 
     case t_INTMOD: case t_PADIC: err(typeer,"gsin");
 
@@ -1979,9 +2014,11 @@
   switch(typ(x))
   {
     case t_INT: case t_FRAC: case t_FRACN:
-      av=avma; p1=cgetr(prec); gaffect(x,p1); tetpil=avma;
-      mpsincos(p1,s,c); gptr[0]=s; gptr[1]=c;
-      gerepilemanysp(av,tetpil,gptr,2);
+      *s=cgetr(prec); *c=cgetr(prec); av=avma; 
+      p1=mpsc_exact(x,prec); 
+      mpsincos(p1,&ps,&pc); 
+      affrr(ps,*s); affrr(pc,*c);
+      avma=av;
       return;
 
     case t_REAL:
@@ -2098,6 +2135,11 @@
       av = avma; gsincos(x,&s,&c,prec);
       return gerepileupto(av, gdiv(s,c));
 
+    case t_INT: case t_FRAC: case t_FRACN:
+      s=cgetr(prec); av=avma; c=mpsc_exact(x,prec); 
+      affrr(mptan(c),s); avma=av;
+      return s;
+
     case t_INTMOD: case t_PADIC: err(typeer,"gtan");
     
     default:
@@ -2144,6 +2186,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:
+      s=cgetr(prec); av=avma; c=mpsc_exact(x,prec); 
+      affrr(mpcotan(c),s); avma=av;
+      return s;
 
     case t_INTMOD: case t_PADIC: err(typeer,"gcotan");
 
Cheers,
Bill.