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*/
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));

```

• References:
• sqr() bug?
• From: Michael Somos <somos@grail.cba.csuohio.edu>
• Re: sqr() bug?
• From: Bill Allombert <Bill.Allombert@math.u-bordeaux.fr>
• Re: sqr() bug?
• From: Michael Somos <somos@grail.cba.csuohio.edu>