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