intnum bug

[Paul Van Wamelen <pvanwam@unix1.sncc.lsu.edu>]

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