| 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