Karim Belabas on Fri, 28 Nov 1997 05:14:09 +0100


[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]

patch8


* This addresses the second "bug" reported by Clemens Heuberger in msg 21
(bad ordering for polred output in compatibility mode).

  In fact it's not a bug, and neither the random number generator, nor the
compatibility mode are involved at all. Only the processor and the way the
compiler chose to distribute the registers: we have a comparison between two
C doubles (quo and thickness, which happen to be equal in this example) whose
accuracy might not be the same depending whether they are in a register or
not (on my test machine, quo is and thickness is not...)

  We don't care if the inequality is reversed due to roundoff errors, except
that we may get the complex roots of certain polynomials ordered differently
on different architectures (or if we use different compilers). Since we want
PARI output to be independent of the machine on which it was produced (with
same long integer size), I'm trying to suppress that (adding a cast).

  I do not know why this did not show up in other benches: before this patch,
a simple polroots(x^4+576) produced two different outputs for me on Linux and
Solaris machines (both 32-bit).

Cheers, Karim.

=========================== patch 8 (2.0.alpha) ============================

*** src/basemath/rootpol.c.orig	Fri Nov 14 04:53:31 1997
--- src/basemath/rootpol.c	Fri Nov 28 01:58:50 1997
***************
*** 1708,1714 ****
  split_1(GEN p, long bitprec, GEN *F, GEN *G)
  {
    long bitprec2,bitprec3,i,imax,n=lgef(p)-3, polreal = isreal(p);
!   double rmax,rmin,thickness;
    GEN q,qq,newq,FF,GG,v,gr,r;
  
    r=gmax_modulus(p,0.01);
--- 1708,1714 ----
  split_1(GEN p, long bitprec, GEN *F, GEN *G)
  {
    long bitprec2,bitprec3,i,imax,n=lgef(p)-3, polreal = isreal(p);
!   double rmax,rmin,thickness,quo;
    GEN q,qq,newq,FF,GG,v,gr,r;
  
    r=gmax_modulus(p,0.01);
***************
*** 1731,1744 ****
      if (3./rmin > thickness)
      {
        rmax=max_modulus(qq,0.05);
!       if (rmax/rmin>thickness)
        {
! 	thickness=rmax/rmin;
! 	newq=qq; globalu=(GEN)v[i];
        }
      }
!     if (thickness>2.) i=imax;
!     if (polreal && (i==1 && thickness>1.5)) i=imax;
    }
    bitprec3=bitprec2+(long)((double) n*log2(3.)+1)+exponorm(newq)-exponorm(q);
    split_2(newq,bitprec3,log(thickness),&FF,&GG);
--- 1731,1744 ----
      if (3./rmin > thickness)
      {
        rmax=max_modulus(qq,0.05);
!       quo = rmax/rmin;
!       if ((float)quo > (float)thickness)
        {
! 	thickness=quo; newq=qq; globalu=(GEN)v[i];
        }
      }
!     if (thickness>2.) break;
!     if (polreal && (i==1 && thickness>1.5)) break;
    }
    bitprec3=bitprec2+(long)((double) n*log2(3.)+1)+exponorm(newq)-exponorm(q);
    split_2(newq,bitprec3,log(thickness),&FF,&GG);
--
Karim Belabas                          e-mail:
Max-Planck-Institut fuer Mathematik       karim@mpim-bonn.mpg.de
Gottfried-Claren-Str. 26               tel:
53225 Bonn (Germany)                      (00 49 228) 402-245