My one-line patch for intnumgauss: what is a real zero ?

[Jacques Gélinas <jacquesg00@hotmail.com>]

The intnumgauss was well written for speed, folding the symmetric table of nodes and weights in half, and finally multiplies the shorter dot-product by two. Now this shortcut will not work with an odd number of nodes since the central node "zero" will then be used twice.

=================================== start cut

tabn(n=2) = {  my(L = pollegendre(n), R = real(polroots(L)), W);
  R = vecsort([ R[k] | k<-[1..n], R[k]>=0 ]);
  W = apply(r->2/(1-r^2)/subst(L',x,r)^2, R);  return([R~,W]);}

fleq(x,y,r=7/8) = if(x==y,1,exponent(normlp(x-y))/exponent(0.) > r);

fleq( [ intnumgauss(x=0,1,1,tabn(k)) | k<-[1..7] ],  \
                [1 + 1/1^2, 1, 1 + 1/(3/2)^2, 1,  1 + 1/(15/8)^2, 1, 1 + 1/(35/16)^2 ]) 

=================================== end paste

With the help of gp2c (a wonderful tool!)
  # echo 'M = if(!r,0,f(bpa-r));' | gp2c
   ===>   if (gequal0(r))  p1 = gen_0;  else p1 = f(gsub(bpa, r));
I patched my copy of intnum.c as shown in the appended file to use the middle zero node only once if present.

This cannot break existing code since intnumgaussinit (default) returns only even tables without 0. However, I worry about the test "gequal0".

What If I generate an odd table with 50 digits and use it later at 100 digits of precision ?
What about a half-baked floating-point zero ? Check fleq(sqrt(1e-50)).
    
Thanks for any advice.

Jacques Gélinas




diff --git a/src/language/intnum.c b/src/language/intnum.c
old mode 100644
new mode 100755
index 551924115..dc3736a1c
@@ -273,7 +293,8 @@ intnumgauss(void *E, GEN (*eval)(void*, GEN), GEN a, GEN b, GEN tab, long prec)
   {
     GEN r = gel(R,i);
     GEN P = eval(E, gadd(bpa, gmul(bma, r)));
-    GEN M = eval(E, gsub(bpa, gmul(bma, r)));
+            /* use middle node only once if present */
+    GEN M = gcmp0(r) ? gen_0 :  eval(E, gsub(bpa, gmul(bma, r)));
     S = gadd(S, gmul(gel(W,i), gadd(P,M)));
     S = gprec_wensure(S, prec2);
   }