Paul Van Wamelen on Mon, 14 Feb 2000 21:55:34 -0600 (CST) |
[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]
intnum bug |
Dear Pari developers, I've been looking into Ilya's intnum(x=0,1,0,3) bug. (it makes the stack overflow). There are several ways the problem could be fixed, but looking at the code (and the examples below) maybe it's the documentation that should be fixed. The intnum algorithm stops when the precision is high enough, that is if log_2(size of the integral) - log_2(latest change in approximation) gets big enough. Of course in case the integrand is zero, this will cause major problems so that the program also stops if the integrand is zero to the current precision and the latest change in approximation is less than half of the current best approximation. When both of these are exactly zero the test is not passed, therefore the bug. This "possible zero integrand" test does cause the following though: ? intnum(x=0,1,sin(x)/10^100)*10^100 - intnum(x=0,1,sin(x)) %1 = -9.469541730020726655672747007 E-15 while the algorithm is in fact capable of the correct thing: ? \p 128 realprecision = 134 significant digits (128 digits displayed) ? intnum(x=0,1,sin(x)/10^100+2) %10 = 2.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000459697694131860282599063392 ? (%-2)*10^100 %11 = 0.45969769413186028259906339255702329140370068700000000000000000000000000000000000000000000000000000000000000000000000000000000000 ? \p 28 realprecision = 28 significant digits ? % - intnum(x=0,1,sin(x)) %12 = 5.995342878624796867000000000 E-28 I therefore suggest the following fix: stop intnum when the integral is zero to the current precision (including exact zero) and forget the bit about "and the latest change in approximation is less than half of the current best approximation", but then apply the pari philosophy of "user beware" and add to the documentation that the integrand should be decent sized. The following patch does this (and fixes a minor typo in the ffinit documentation). Sincerely, Paul van Wamelen. "cvs diff -c" output: cvs server: Diffing . cvs server: Diffing Odos cvs server: Diffing config cvs server: Diffing doc Index: doc/usersch3.tex =================================================================== RCS file: /home/megrez/cvsroot/pari/doc/usersch3.tex,v retrieving revision 1.42 diff -c -r1.42 usersch3.tex *** doc/usersch3.tex 2000/02/08 20:23:56 1.42 --- doc/usersch3.tex 2000/02/14 22:52:42 *************** *** 3230,3236 **** $n$ which is irreducible over $\F_p$. For instance if \kbd{P = ffinit(3,2,y)}, you can represent elements in $\F_{3^2}$ as polmods modulo \kbd{P}. This function is rather crude and expects $p$ to be ! relatively small ($p < 2^31$). \syn{ffinit}{p,n,v}, where $v$ is a variable number. --- 3230,3236 ---- $n$ which is irreducible over $\F_p$. For instance if \kbd{P = ffinit(3,2,y)}, you can represent elements in $\F_{3^2}$ as polmods modulo \kbd{P}. This function is rather crude and expects $p$ to be ! relatively small ($p < 2^{31}$). \syn{ffinit}{p,n,v}, where $v$ is a variable number. *************** *** 5847,5853 **** this can be accomplished with a simple change of variable. Furthermore, for improper integrals, where one or both of the limits of integration are plus or minus infinity, the function must decrease sufficiently rapidly at ! infinity. This can often be accomplished through integration by parts. Note that \idx{infinity} can be represented with essentially no loss of accuracy by 1e4000. However beware of real underflow when dealing with --- 5847,5857 ---- this can be accomplished with a simple change of variable. Furthermore, for improper integrals, where one or both of the limits of integration are plus or minus infinity, the function must decrease sufficiently rapidly at ! infinity. This can often be accomplished through integration by parts. ! Finally, the function to be integrated should not be very small ! (compared to the current precision) on the entire interval. This can ! of course be accomplished by just multiplying by an appropriate ! constant. Note that \idx{infinity} can be represented with essentially no loss of accuracy by 1e4000. However beware of real underflow when dealing with cvs server: Diffing emacs cvs server: Diffing examples cvs server: Diffing misc cvs server: Diffing src cvs server: Diffing src/basemath cvs server: Diffing src/gp cvs server: Diffing src/graph cvs server: Diffing src/headers cvs server: Diffing src/kernel cvs server: Diffing src/kernel/alpha cvs server: Diffing src/kernel/hppa cvs server: Diffing src/kernel/ix86 cvs server: Diffing src/kernel/m68k cvs server: Diffing src/kernel/none cvs server: Diffing src/kernel/ppc cvs server: Diffing src/kernel/sparcv7 cvs server: Diffing src/kernel/sparcv8 cvs server: Diffing src/kernel/sparcv9 cvs server: Diffing src/language Index: src/language/sumiter.c =================================================================== RCS file: /home/megrez/cvsroot/pari/src/language/sumiter.c,v retrieving revision 1.8 diff -c -r1.8 sumiter.c *** src/language/sumiter.c 1999/10/29 09:04:52 1.8 --- src/language/sumiter.c 2000/02/14 22:52:44 *************** *** 888,894 **** { tetpil=avma; ss=polint_i(h+j-KLOC,s+j-KLOC,gzero,KLOC+1,&dss); j1=gexpo(ss); j2=gexpo(dss); lim=bit_accuracy(prec)-j-6; ! if (j1-j2 > lim || (j1 < -lim && j2<j1-1)) { pop_val(ep); if (gcmp0(gimag(ss))) ss=greal(ss); tetpil=avma; return gerepile(av,tetpil,gmulsg(sig,ss)); --- 888,894 ---- { tetpil=avma; ss=polint_i(h+j-KLOC,s+j-KLOC,gzero,KLOC+1,&dss); j1=gexpo(ss); j2=gexpo(dss); lim=bit_accuracy(prec)-j-6; ! if (j1-j2 > lim || j1 < -lim) { pop_val(ep); if (gcmp0(gimag(ss))) ss=greal(ss); tetpil=avma; return gerepile(av,tetpil,gmulsg(sig,ss)); *************** *** 935,941 **** { tetpil=avma; ss=polint_i(h+j-KLOC,s+j-KLOC,gzero,KLOC+1,&dss); j1=gexpo(ss); j2=gexpo(dss); lim=bit_accuracy(prec)-(3*j/2)-6; ! if ( j1-j2 > lim || (j1 < -lim && j2<j1-1)) { pop_val(ep); if (gcmp0(gimag(ss))) ss=greal(ss); tetpil=avma; return gerepile(av,tetpil,gmulsg(sig,ss)); --- 935,941 ---- { tetpil=avma; ss=polint_i(h+j-KLOC,s+j-KLOC,gzero,KLOC+1,&dss); j1=gexpo(ss); j2=gexpo(dss); lim=bit_accuracy(prec)-(3*j/2)-6; ! if ( j1-j2 > lim || j1 < -lim) { pop_val(ep); if (gcmp0(gimag(ss))) ss=greal(ss); tetpil=avma; return gerepile(av,tetpil,gmulsg(sig,ss)); *************** *** 993,999 **** { tetpil=avma; ss=polint_i(h+j-KLOC,s+j-KLOC,gzero,KLOC+1,&dss); j1=gexpo(ss); j2=gexpo(dss); lim=bit_accuracy(prec)-(3*j/2)-6; ! if (j1-j2 > lim || (j1 < -lim && j2 < j1-1)) { pop_val(ep); if (gcmp0(gimag(ss))) ss=greal(ss); tetpil=avma; return gerepile(av,tetpil,gmulsg(sig,ss)); --- 993,999 ---- { tetpil=avma; ss=polint_i(h+j-KLOC,s+j-KLOC,gzero,KLOC+1,&dss); j1=gexpo(ss); j2=gexpo(dss); lim=bit_accuracy(prec)-(3*j/2)-6; ! if (j1-j2 > lim || j1 < -lim) { pop_val(ep); if (gcmp0(gimag(ss))) ss=greal(ss); tetpil=avma; return gerepile(av,tetpil,gmulsg(sig,ss)); cvs server: Diffing src/modules cvs server: Diffing src/test cvs server: Diffing src/test/32 cvs server: Diffing src/test/64 cvs server: Diffing src/test/in