Re: roots function

[Karim BELABAS <Karim.Belabas@math.u-psud.fr>]

On Sat, 20 Oct 2001, Jean-Marc Sac-Epée wrote:
> I have strange problems when using roots function with some peculiar
> polynomials. With most polynomials, roots is OK, but I can show some
> examples which stop the program.

This has been corrected some time ago in the unstable version. Since it
didn't have any apparent adverse effect, I have merged the change back to
the stable CVS branch.

The patch follows, please check that it solves your problem without creating
new ones. The stable and unstable versions have diverged significantly and I
tried to minimize the patch [original diff was over 2000 lines long]

    Karim.

P.S: Applies to pari-2.1.2 or 2.1.1

Index: src/basemath/rootpol.c
===================================================================
RCS file: /home/megrez/cvsroot/pari/src/basemath/rootpol.c,v
retrieving revision 1.21.2.1
retrieving revision 1.21.2.2
diff -c -r1.21.2.1 -r1.21.2.2
*** src/basemath/rootpol.c	2001/07/03 23:24:18	1.21.2.1
--- src/basemath/rootpol.c	2001/10/25 15:04:15	1.21.2.2
***************
*** 1307,1312 ****
--- 1307,1372 ----
      if (i==0) return r;
      aux=gmul(pui,aux);
    }
+ }
+
+ static GEN
+ compute_radius(GEN* radii, GEN p, long k, double aux, double *delta)
+ {
+   long i, n = lgef(p)-3;
+   GEN rmin,rmax,p1;
+   if (k>1)
+   {
+     i=k-1; while (i>0 && !signe(radii[i])) i--;
+     rmin = pre_modulus(p,k,aux, radii[i], radii[k]);
+   }
+   else /* k=1 */
+     rmin = min_modulus(p,aux);
+   affrr(rmin, radii[k]);
+
+   if (k+1<n)
+   {
+     i=k+2; while (i<=n && !signe(radii[i])) i++;
+     rmax = pre_modulus(p,k+1,aux, radii[k+1], radii[i]);
+   }
+   else /* k+1=n */
+     rmax = max_modulus(p,aux);
+   affrr(rmax, radii[k+1]);
+
+   p1 = radii[k];
+   for (i=k-1; i>=1; i--)
+   {
+     if (!signe(radii[i]) || cmprr(radii[i], p1) > 0)
+       affrr(p1, radii[i]);
+     else
+       p1 = radii[i];
+   }
+   p1 = radii[k+1];
+   for (i=k+1; i<=n; i++)
+   {
+     if (!signe(radii[i]) || cmprr(radii[i], p1) < 0)
+       affrr(p1, radii[i]);
+     else
+       p1 = radii[i];
+   }
+   *delta = rtodbl(gmul2n(mplog(divrr(rmax,rmin)), -1));
+   if (*delta > 1.) *delta = 1.;
+   return mpsqrt(mulrr(rmin,rmax));
+ }
+
+ static GEN
+ update_radius(GEN *radii, GEN rho, double *par, double *par2)
+ {
+   GEN p1, invrho = ginv(rho);
+   long i, n = lg(radii);
+   double t, param = 0., param2 = 0.;
+   for (i=1; i<n; i++)
+   {
+     affrr(mulrr(radii[i], invrho), radii[i]);
+     p1 = ginv(subsr(1, radii[i]));
+     t = fabs(rtodbl(p1));
+     param += t; if (t > 1.) param2 += log2(t);
+   }
+   *par = param; *par2 = param2; return invrho;
  }

  /* apply the conformal mapping then split from U */
***************
*** 1314,1321 ****
  conformal_mapping(GEN p, long k, long bitprec, double aux,GEN *F,GEN *G)
  {
    long bitprec2,n=lgef(p)-3,decprec,i,ltop = avma, av;
!   GEN q,FF,GG,a,R,aux2, *gptr[3];
!   GEN rmin,rmax,rho,invrho;
    double delta,param,param2;

    bitprec2=bitprec+(long) (n*(2.*log2(2.732)+log2(1.5)))+1;
--- 1374,1381 ----
  conformal_mapping(GEN p, long k, long bitprec, double aux,GEN *F,GEN *G)
  {
    long bitprec2,n=lgef(p)-3,decprec,i,ltop = avma, av;
!   GEN q,FF,GG,a,R, *gptr[3];
!   GEN rho,invrho;
    double delta,param,param2;

    bitprec2=bitprec+(long) (n*(2.*log2(2.732)+log2(1.5)))+1;
***************
*** 1335,1395 ****
        affrr(mpsqrt(addsr(1,p1)), radius[i]);
        avma = a;
      }
-   if (k>1)
-   {
-     i=k-1; while (i>0 && !signe(radius[i])) i--;
-     rmin = pre_modulus(q,k,aux/10., radius[i], radius[k]);
-   }
-   else
-     rmin = min_modulus(q,aux/10.);
-   affrr(rmin, radius[k]);

!   if (k+1<n)
!   {
!     i=k+2; while (i<=n && !signe(radius[i])) i++;
!     rmax = pre_modulus(q,k+1,aux/10., radius[k+1], radius[i]);
!   }
!   else /* k+1=n */
!     rmax = max_modulus(q,aux/10.);
!   affrr(rmax, radius[k+1]);

-   rho = mpsqrt(mulrr(rmin,rmax));
-   delta = rtodbl(gmul2n(mplog(divrr(rmax,rmin)), -1));
-   if (delta>1.) delta=1.;
-
    bitprec2 += (long) (((double)n) * fabs(log2ir(rho)) + 1.);
-
-   invrho = ginv(rho);
    R=mygprec(invrho,bitprec2);
    q=scalepol(q,R,bitprec2);
    gptr[0]=&q; gptr[1]=&R; gptr[2]=&invrho;
    gerepilemany(av,gptr,3);

-   aux2 = radius[k];
-   for (i=k-1; i>=1; i--)
-   {
-     if (!signe(radius[i]) || cmprr(radius[i], aux2) > 0)
-       affrr(aux2, radius[i]);
-     else
-       aux2 = radius[i];
-   }
-   aux2 = radius[k+1];
-   for (i=k+1; i<=n; i++)
-   {
-     if (!signe(radius[i]) || cmprr(radius[i], aux2) < 0)
-       affrr(aux2, radius[i]);
-     else
-       aux2 = radius[i];
-   }
-   param=0.; param2=0.;
-   for (i=1; i<=n; i++)
-   {
-     double t;
-     affrr(mulrr(radius[i], invrho), radius[i]);
-     aux2 = ginv(subsr(1, radius[i]));
-     t = fabs(rtodbl(aux2));
-     param += t; if (t > 1.) param2 += log2(t);
-   }
    optimize_split(q,k,delta,bitprec2,&FF,&GG,param,param2);
    bitprec2 += n; R = ginv(R);
    FF=scalepol(FF,R,bitprec2);
--- 1395,1410 ----
        affrr(mpsqrt(addsr(1,p1)), radius[i]);
        avma = a;
      }

!   rho = compute_radius(radius, q,k,aux/10., &delta);
!   invrho = update_radius(radius, rho, &param, &param2);

    bitprec2 += (long) (((double)n) * fabs(log2ir(rho)) + 1.);
    R=mygprec(invrho,bitprec2);
    q=scalepol(q,R,bitprec2);
    gptr[0]=&q; gptr[1]=&R; gptr[2]=&invrho;
    gerepilemany(av,gptr,3);

    optimize_split(q,k,delta,bitprec2,&FF,&GG,param,param2);
    bitprec2 += n; R = ginv(R);
    FF=scalepol(FF,R,bitprec2);
***************
*** 1418,1424 ****
    GEN rmin,rmax,rho,invrho;
    double kappa,aux,delta,param,param2;
    long n=lgef(p)-3,i,j,k,bitprec2;
!   GEN q,FF,GG,R,aux2;

    radius = (GEN*) cgetg(n+1, t_VEC);
    for (i=1; i<n; i++) radius[i]=realzero(3);
--- 1433,1439 ----
    GEN rmin,rmax,rho,invrho;
    double kappa,aux,delta,param,param2;
    long n=lgef(p)-3,i,j,k,bitprec2;
!   GEN q,FF,GG,R;

    radius = (GEN*) cgetg(n+1, t_VEC);
    for (i=1; i<n; i++) radius[i]=realzero(3);
***************
*** 1456,1508 ****

    if (step4==0)
    {
!     aux /= 10.;
!     if (k>1)
!     {
!       i=k-1; while (i>0 && !signe(radius[i])) i--;
!       rmin = pre_modulus(p,k,aux, radius[i], radius[k]);
!     }
!     else /* k=1 */
!       rmin = min_modulus(p,aux);
!     affrr(rmin, radius[k]);
!
!     if (k+1<n)
!     {
!       i=k+2; while (i<=n && !signe(radius[i])) i++;
!       rmax = pre_modulus(p,k+1,aux, radius[k+1], radius[i]);
!     }
!     else /* k+1=n */
!       rmax = max_modulus(p,aux);
!     affrr(rmax, radius[k+1]);
!
!     rho = mpsqrt(mulrr(rmin,rmax));
!     delta = rtodbl(gmul2n(mplog(divrr(rmax, rmin)), -1));
!
!     aux2 = radius[k];
!     for (i=k-1; i>=1; i--)
!     {
!       if (!signe(radius[i]) || cmprr(radius[i], aux2) > 0)
!         affrr(aux2, radius[i]);
!       else
!         aux2 = radius[i];
!     }
!     aux2 = radius[k+1];
!     for (i=k+1; i<=n; i++)
!     {
!       if (!signe(radius[i]) || cmprr(radius[i], aux2) < 0)
!         affrr(aux2, radius[i]);
!       else
!         aux2 = radius[i];
!     }
!     invrho = ginv(rho);
!     param = 0.;param2 = 0.;
!     for (i=1; i<=n; i++)
!     {
!       affrr(mulrr(radius[i], invrho), radius[i]);
!       aux2 = ginv(subsr(1, radius[i]));
!       param += fabs(rtodbl(aux2));
!       if (cmprs(aux2, 1) > 0) param2 += log2ir(aux2);
!     }

      bitprec2 = bitprec + (long) ((double)n * fabs(log2ir(rho)) + 1.);
      R = mygprec(invrho,bitprec2);
--- 1471,1478 ----

    if (step4==0)
    {
!     rho = compute_radius(radius, p, k, aux/10., &delta);
!     invrho = update_radius(radius, rho, &param, &param2);

      bitprec2 = bitprec + (long) ((double)n * fabs(log2ir(rho)) + 1.);
      R = mygprec(invrho,bitprec2);

-- 
Karim Belabas                    email: Karim.Belabas@math.u-psud.fr
Dep. de Mathematiques, Bat. 425
Universite Paris-Sud             Tel: (00 33) 1 69 15 57 48
F-91405 Orsay (France)           Fax: (00 33) 1 69 15 60 19
--
PARI/GP Home Page: http://www.parigp-home.de/