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

 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);
!
!
!     c1 = quickmulii(c1+2,c2+2, lgefint(c1)-2,lgefint(c2)-2);
}
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);
!
!       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);
!     }
}
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,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);
!   }
!   else