Karim Belabas on Mon, 30 Nov 1998 14:05:43 +0100 (MET)


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

Re: Discriminant of non-monic polynomials


[Igor:]
> Hi, I've encountered several problems with poldisc and nfdisc
> commands.  Here's the demonstration:
> ------------------------------------------------------------------------------
> ? g(pol)=subst(pol,x,2*x+1) \\ If pol is f(x) then g(pol) is f(2*x+1)
> ? poldisc(g(polcyclo(57)))
>   ***   Warning: normalizing a polynomial with 0 leading term.
>   ***   bus error: bug in GP (please report).
> 
> ? nfdisc(g(x^12-x-1))
>   ***   impossible assignment I-->S
> ? nfdisc(g(polcyclo(23)))
>   ***   the PARI stack overflows !!!
> 
>   ***   Warning: doubling stack size; new stack = 8000000.
> ? nfdisc(g(polcyclo(23)))
>   ***   the PARI stack overflows !!!
> 
>   ***   Warning: doubling stack size; new stack = 16000000.                   
> 
> ??nfdisc tells me "preferably monic".  How do I interprete this?  Does
> it mean that nfdisc will give reliable output only for monic
> polynomials?

No: we use the obvious (and highly inefficient) change of variable to reduce
to monic form (it's theoretically possible to work directly with non-monic
polynomials, but it will be a _real_ pain to adapt the code...). [then call
polred to try and lower the huge index which we just introduced].

It's usually a bad idea to use non-monic polynomials, since internally we are
going to transform them and work on a different polynomial. It's more
efficient to do the change of variable yourself and reduce the output via
polred / polredabs. Hence the warning.

                       ****** IMPORTANT ******

On the other hand, you exhibited a real (kernel) bug which I introduced in
2.0.12 (operations on rationals incorrectly assume that the stack cannot grow
more than a certain bound, which is not true anymore as soon as Karatsuba
multiplication is called, i.e operands bigger than 155 digits...). Apply the
patch at the end.

Karim.

*** src/basemath/gen1.c.orig	Fri Nov  6 16:08:14 1998
--- src/basemath/gen1.c	Wed Nov 18 17:46:38 1998
***************
*** 701,706 ****
--- 701,712 ----
  #define fix_frac_if_int(z) if (is_pm1(z[2]))\
    z = gerepileupto((long)(z+3), (GEN)z[1]);
  
+ /* assume z[1] was created last */
+ #define fix_frac_if_int_GC(z,tetpil) if (is_pm1(z[2]))\
+   z = gerepileupto((long)(z+3), (GEN)z[1]);\
+ else\
+   gerepilemanyvec((long)z, tetpil, z+1, 2);
+ 
  GEN
  fix_rfrac_if_pol(GEN x, GEN y)
  {
***************
*** 883,889 ****
  	  case t_FRAC:
              if (!signe(x)) return gzero;
              z=cgetg(3,t_FRAC);
-             (void)new_chunk((lgefint(x)+lgefint(y[1])+lgefint(y[2]))<<1);
              p1 = mppgcd(x,(GEN)y[2]);
              if (is_pm1(p1))
              {
--- 889,894 ----
***************
*** 893,902 ****
              }
              else
              {
!               x = divii(x,p1); avma = (long)z;
                z[2] = ldivii((GEN)y[2], p1);
                z[1] = lmulii((GEN)y[1], x);
!               fix_frac_if_int(z);
              }
              return z;
  
--- 898,907 ----
              }
              else
              {
!               x = divii(x,p1); tetpil = avma;
                z[2] = ldivii((GEN)y[2], p1);
                z[1] = lmulii((GEN)y[1], x);
!               fix_frac_if_int_GC(z,tetpil);
              }
              return z;
  
***************
*** 988,1002 ****
              GEN x1 = (GEN)x[1], x2 = (GEN)x[2];
              GEN y1 = (GEN)y[1], y2 = (GEN)y[2];
              z=cgetg(3,t_FRAC);
-             (void)new_chunk(lgefint(x1)+lgefint(x2)+lgefint(y1)+lgefint(y2));
              p1 = mppgcd(x1, y2);
              if (!is_pm1(p1)) { x1 = divii(x1,p1); y2 = divii(y2,p1); }
              p1 = mppgcd(x2, y1);
              if (!is_pm1(p1)) { x2 = divii(x2,p1); y1 = divii(y1,p1); }
!             avma = (long)z;
              z[2] = lmulii(x2,y2);
              z[1] = lmulii(x1,y1);
!             fix_frac_if_int(z); return z;
            }
  	  case t_FRACN: z=cgetg(3,t_FRACN);
  	    z[1]=lmulii((GEN)x[1],(GEN)y[1]);
--- 993,1006 ----
              GEN x1 = (GEN)x[1], x2 = (GEN)x[2];
              GEN y1 = (GEN)y[1], y2 = (GEN)y[2];
              z=cgetg(3,t_FRAC);
              p1 = mppgcd(x1, y2);
              if (!is_pm1(p1)) { x1 = divii(x1,p1); y2 = divii(y2,p1); }
              p1 = mppgcd(x2, y1);
              if (!is_pm1(p1)) { x2 = divii(x2,p1); y1 = divii(y1,p1); }
!             tetpil = avma;
              z[2] = lmulii(x2,y2);
              z[1] = lmulii(x1,y1);
!             fix_frac_if_int_GC(z,tetpil); return z;
            }
  	  case t_FRACN: z=cgetg(3,t_FRACN);
  	    z[1]=lmulii((GEN)x[1],(GEN)y[1]);
***************
*** 1418,1438 ****
  	  case t_FRAC:
              if (!signe(x)) return gzero;
              z=cgetg(3,t_FRAC);
-             (void)new_chunk(lgefint(x)+lgefint(y[2])+ (lgefint(y[1])<<1));
              p1 = mppgcd(x,(GEN)y[1]);
              if (is_pm1(p1))
              {
!               avma = (long)z;
                z[2] = licopy((GEN)y[1]);
              }
              else
              {
!               x = divii(x,p1); avma = (long)z;
                z[2] = ldivii((GEN)y[1], p1);
              }
              z[1] = lmulii((GEN)y[2], x);
              fix_frac(z);
!             fix_frac_if_int(z); return z;
  		
  	  case t_FRACN:
              if (!signe(x)) return gzero;
--- 1422,1445 ----
  	  case t_FRAC:
              if (!signe(x)) return gzero;
              z=cgetg(3,t_FRAC);
              p1 = mppgcd(x,(GEN)y[1]);
              if (is_pm1(p1))
              {
!               avma = (long)z; tetpil = 0;
                z[2] = licopy((GEN)y[1]);
              }
              else
              {
!               x = divii(x,p1); tetpil = avma;
                z[2] = ldivii((GEN)y[1], p1);
              }
              z[1] = lmulii((GEN)y[2], x);
              fix_frac(z);
!             if (tetpil)
!             { fix_frac_if_int_GC(z,tetpil); }
!             else
!               fix_frac_if_int(z);
!             return z;
  		
  	  case t_FRACN:
              if (!signe(x)) return gzero;
***************
*** 1526,1548 ****
            z = cgetg(3, tx);
  	  if (tx == t_FRAC)
            {
-             (void)new_chunk((lgefint(y)+lgefint(x[2])) + (lgefint(x[1])<<1));
              p1 = mppgcd(y,(GEN)x[1]);
              if (is_pm1(p1))
              {
!               avma = (long)z;
                z[1] = licopy((GEN)x[1]);
              }
              else
              {
!               y = divii(y,p1); avma = (long)z;
                z[1] = ldivii((GEN)x[1], p1);
              }
            }
            else
              z[1] = licopy((GEN)x[1]);
            z[2] = lmulii((GEN)x[2],y);
!           fix_frac(z); return z;
  
  	  case t_REAL:
  	    l=avma; p1=cgetg(ly,ty); gaffect(x,p1);
--- 1533,1559 ----
            z = cgetg(3, tx);
  	  if (tx == t_FRAC)
            {
              p1 = mppgcd(y,(GEN)x[1]);
              if (is_pm1(p1))
              {
!               avma = (long)z; tetpil = 0;
                z[1] = licopy((GEN)x[1]);
              }
              else
              {
!               y = divii(y,p1); tetpil = avma;
                z[1] = ldivii((GEN)x[1], p1);
              }
            }
            else
+           {
+             tetpil = 0;
              z[1] = licopy((GEN)x[1]);
+           }
            z[2] = lmulii((GEN)x[2],y);
!           fix_frac(z); 
!           if (tetpil) fix_frac_if_int_GC(z,tetpil);
!           return z;
  
  	  case t_REAL:
  	    l=avma; p1=cgetg(ly,ty); gaffect(x,p1);
***************
*** 1561,1576 ****
              {
                GEN x1 = (GEN)x[1], x2 = (GEN)x[2];
                GEN y1 = (GEN)y[1], y2 = (GEN)y[2];
-               (void)new_chunk(lgefint(x1)+lgefint(x2)+lgefint(y1)+lgefint(y2));
                p1 = mppgcd(x1, y1);
                if (!is_pm1(p1)) { x1 = divii(x1,p1); y1 = divii(y1,p1); }
                p1 = mppgcd(x2, y2);
                if (!is_pm1(p1)) { x2 = divii(x2,p1); y2 = divii(y2,p1); }
!               avma = (long)z;
                z[2] = lmulii(x2,y1);
                z[1] = lmulii(x1,y2);
                fix_frac(z);
!               fix_frac_if_int(z);
              }
              else
              {
--- 1572,1586 ----
              {
                GEN x1 = (GEN)x[1], x2 = (GEN)x[2];
                GEN y1 = (GEN)y[1], y2 = (GEN)y[2];
                p1 = mppgcd(x1, y1);
                if (!is_pm1(p1)) { x1 = divii(x1,p1); y1 = divii(y1,p1); }
                p1 = mppgcd(x2, y2);
                if (!is_pm1(p1)) { x2 = divii(x2,p1); y2 = divii(y2,p1); }
!               tetpil = avma;
                z[2] = lmulii(x2,y1);
                z[1] = lmulii(x1,y2);
                fix_frac(z);
!               fix_frac_if_int_GC(z,tetpil);
              }
              else
              {
*** src/kernel/none/mp.c.orig	Fri Nov  6 16:08:36 1998
--- src/kernel/none/mp.c	Wed Nov 18 16:34:35 1998
***************
*** 1870,1876 ****
    av=avma; a0=a+na; n0a=n0;
    while (!*a0 && n0a) { a0++; n0a--; }
  
!   if (nb > n0)
    { /* nb <= na <= n0 */
      GEN b0,c1,c2;
      long n0b;
--- 1870,1876 ----
    av=avma; a0=a+na; n0a=n0;
    while (!*a0 && n0a) { a0++; n0a--; }
  
!   if (n0a && nb > n0)
    { /* nb <= na <= n0 */
      GEN b0,c1,c2;
      long n0b;
***************
*** 1879,1892 ****
      c = quickmulii(a,b,na,nb);
      b0 = b+nb; n0b = n0;
      while (!*b0 && n0b) { b0++; n0b--; }
!     c0 = quickmulii(a0,b0, n0a,n0b);
! 
!     c2 = addiispec(a0,a, n0a,na);
!     c1 = addiispec(b0,b, n0b,nb);
! 
!     c1 = quickmulii(c1+2,c2+2, lgefint(c1)-2,lgefint(c2)-2);
!     c2 = addii(c0,c); setsigne(c2,-1);
!     c = addshiftw(c, addii(c1,c2), n0);
    }
    else
    {
--- 1879,1901 ----
      c = quickmulii(a,b,na,nb);
      b0 = b+nb; n0b = n0;
      while (!*b0 && n0b) { b0++; n0b--; }
!     if (n0b)
!     {
!       c0 = quickmulii(a0,b0, n0a,n0b);
! 
!       c2 = addiispec(a0,a, n0a,na);
!       c1 = addiispec(b0,b, n0b,nb);
!       c1 = quickmulii(c1+2,c2+2, lgefint(c1)-2,lgefint(c2)-2);
!       c2 = addiispec(c0+2, c+2, lgefint(c0)-2,lgefint(c) -2);
! 
!       c1 = subiispec(c1+2,c2+2, lgefint(c1)-2,lgefint(c2)-2);
!     }
!     else
!     {
!       c0 = gzero;
!       c1 = quickmulii(a0,b, n0a,nb);
!     }
!     c = addshiftw(c,c1, n0);
    }
    else
    {
***************
*** 1923,1934 ****
    i=(na>>1); n0=na-i; na=i;
    av=avma; a0=a+na; n0a=n0;
    while (!*a0 && n0a) { a0++; n0a--; }
- 
    c = quicksqri(a,na);
!   c0 = quicksqri(a0,n0a);
!   c1 = shifti(quickmulii(a0,a, n0a,na),1);
!   c = addshiftw(c,c1, n0);
!   c = addshiftw(c,c0, n0); i=lgefint(c);
    avma = (long)(((GEN)av) -  i);
    c0 = (GEN)avma; while (--i >= 0) c0[i]=c[i];
    return c0;
--- 1932,1948 ----
    i=(na>>1); n0=na-i; na=i;
    av=avma; a0=a+na; n0a=n0;
    while (!*a0 && n0a) { a0++; n0a--; }
    c = quicksqri(a,na);
!   if (n0a)
!   {
!     c0 = quicksqri(a0,n0a);
!     c1 = shifti(quickmulii(a0,a, n0a,na),1);
!     c = addshiftw(c,c1, n0);
!     c = addshiftw(c,c0, n0);
!   }
!   else
!     c = addshiftw(c,gzero,n0<<1);
!   i = lgefint(c);
    avma = (long)(((GEN)av) -  i);
    c0 = (GEN)avma; while (--i >= 0) c0[i]=c[i];
    return c0;

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