Bill Allombert on Sat, 30 Oct 1999 16:54:31 +0200 (MET DST)


[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]

Re: sqr() bug and acos bug.


The "PARI philosophy" (do not forget PARI name comes from Blaise
Pascal) state that "... operations and functions should, firstly, give
as exact a result as possible ...".

Suppose X= 5+ e_1*2^10 +e_2*2^11+e_3+2^12+ a true 2-adic number (infinite developpement)
where e_i=0 or 1
and x=X+O(2^10) an approximation 

X^2=1 + 2^3 + 2^4 + O(2^11) independantly of e_1
so PARI *must* return 1 + 2^3 + 2^4 + O(2^11) for sqr(x),x^2 and x*x

But in (5+O(2^10))*(5+O(2^10)) we do not know if the 10th digits of the two factors is the same, for example
X= 5+ 2^10 +e_2*2^11+++
and 
Y= 5+ 0*2^10 +f_2^2^11+++

X^2=1 + 2^3 + 2^4 + O(2^11) 
XY=1 + 2^3 + 2^4 + 2^10 + O(2^11) 

So PARI must return 1 + 2^3 + 2^4 + O(2^10) to be safe

So this is the explanation for x*x.

Well, but I do not like the idea that a PARI function can return
different results depending of the memory address of it's
arguments. It ruins all tentatives to predicts the results.
Fortunately t_PADIC are almosts never used in library mode, mainly because there are too slow.

Secondly, such nice behaviour is not implemented for p-ading numbers
 and exponents p . e.g. : 4+O(3^5))^3=1 + 3^2 + 2*3^3 + O(3^6) and gp
 return only 1 + 3^2 + 2*3^3 + O(3^5)


If you like 2-adics, check if the sqrt function is correct (remember
there is 4 square roots for 2-adics and PARI returns only one).

> Message-Id: <199910292001.QAA13852@grail.cba.csuohio.edu>
> Message-Id: <199910291922.PAA12938@grail.cba.csuohio.edu>
____
Here a patch which corrects acos and atan, and doc of abs.
I have updated the CVS repository.
? atan(O(1))
  ***   division by zero in gdiv, gdivgs or ginv   
? \ps1
? acos(x)
  ***   bus error: bug in GP (please report).
Now we get:

? atan(O(1))
%1 = O(1)
? \ps1
   seriesprecision = 1 significant terms
? acos(x)
%2 = 1.570796326794896619231321691 + O(x)     

Bill.

Index: doc/usersch3.tex
===================================================================
RCS file: /home/megrez/cvsroot/pari/doc/usersch3.tex,v
retrieving revision 1.16
diff -u -r1.16 usersch3.tex
--- doc/usersch3.tex	1999/10/29 09:04:51	1.16
+++ doc/usersch3.tex	1999/10/30 13:25:13
@@ -794,9 +794,10 @@
 name \teb{gpi} (with no parentheses), use instead $\teb{constpi}(\var{prec})$.
 
 \subsecidx{abs}$(x)$: absolute value of $x$ (modulus if $x$ is complex).
-Polynomials, power series and rational functions are not allowed.
+Power series and rational functions are not allowed.
 Contrary to most transcendental functions, an integer is \var{not}
 converted to a real number before applying \kbd{abs}.
+If $x$ is a polynomial, returns $-x$ if the leading coefficient is a negative real number (including integer and rational) else returns $x$.
 
 \syn{gabs}{x,\var{prec}}.
 
Index: src/basemath/trans2.c
===================================================================
RCS file: /home/megrez/cvsroot/pari/src/basemath/trans2.c,v
retrieving revision 1.2
diff -u -r1.2 trans2.c
--- src/basemath/trans2.c	1999/09/23 17:50:57	1.2
+++ src/basemath/trans2.c	1999/10/30 13:25:14
@@ -130,6 +130,8 @@
 
     case t_SER:
       av=avma; if (valp(x)<0) err(negexper,"gatan");
+      if (lg(x)==2) return gcopy(x);
+      /*Now we can assume lg(x)>2*/
       p1=gdiv(derivser(x), gaddsg(1,gsqr(x)));
       y=integ(p1,varn(x)); if (valp(x)) return gerepileupto(av,y);
 
@@ -215,7 +217,7 @@
 
     case t_SER:
       if (gcmp0(x)) return gcopy(x);
-
+      /*Now, we can assume lg(x)>2*/
       av=avma; if (valp(x)<0) err(negexper,"gasin");
       p1=gdiv(derivser(x), gsqrt(gsubsg(1,gsqr(x)),prec));
       y=integ(p1,varn(x)); if (valp(x)) return gerepileupto(av,y);
@@ -317,12 +319,16 @@
 
     case t_SER: av=avma;
       if (valp(x)<0) err(negexper,"gacos");
-      p1=integ(gdiv(derivser(x), gsqrt(gsubsg(1,gsqr(x)),prec)),varn(x));
-      if (gcmp1((GEN)x[2]) && !valp(x))
+      if (lg(x)>2)
       {
-	tetpil=avma; return gerepile(av,tetpil,gneg(p1));
-      }
-      if (valp(x)) { y=mppi(prec); setexpo(y,0); }
+	p1=integ(gdiv(derivser(x), gsqrt(gsubsg(1,gsqr(x)),prec)),varn(x));
+	if (gcmp1((GEN)x[2]) && !valp(x))/*x=1+O(x^k);k>=1*/
+	{
+	  tetpil=avma; return gerepile(av,tetpil,gneg(p1));
+	}
+      } 
+      else p1=x;
+      if (lg(x)==2 || valp(x)) { y=mppi(prec); setexpo(y,0); }
       else y=gacos((GEN)x[2],prec);
       tetpil=avma; return gerepile(av,tetpil,gsub(y,p1));