Karim BELABAS on Thu, 21 Jan 1999 15:36:41 +0100 (MET)

 Re: poldisc

```[Ilya:]
> b) I tried to do
>
> p=a*x^3+b*x^2*y+c*x*y^2+d*y^3 + e*x^2+f*x*y+g*y^2 + h*x+i*y + h
> di=poldisc(p)
>
> and it will not work with 16M stack, and does not finish in half an hour
> anyway.
>
> Should it be *that* resource-intensive?  I think I may be able to calculate
> given a sheet of paper...

PARI's internal representation for polynomials... (highly non-symetrical,
optimized for 1 or 2 variables).

It is instantaneous with default stack, *IF* you make sure the
important variables have high priority:
gp> y; p=a*x^3+b*x^2*y+c*x*y^2+d*y^3 + e*x^2+f*x*y+g*y^2 + h*x+i*y + h
^^^ (now x and y are priviledged)

> a) It looks like you cannot advice poldisc() which variable you are
> interested in.

True; easy to add, though. The patch below does what you want (see online
help).

Karim.

*** src/basemath/polarit2.c.orig	Tue Jan 19 13:51:09 1999
--- src/basemath/polarit2.c	Thu Jan 21 15:20:43 1999
***************
*** 2049,2055 ****
}

GEN
! discsr(GEN x)
{
long av=avma,tetpil,tx=typ(x),i;
GEN z,p1,p2;
--- 2062,2068 ----
}

GEN
! poldisc0(GEN x, long v)
{
long av=avma,tetpil,tx=typ(x),i;
GEN z,p1,p2;
***************
*** 2058,2064 ****
{
case t_POL:
if (gcmp0(x)) return gzero;
!       p1 = subres(x, derivpol(x));
p2 = leading_term(x); if (!gcmp1(p2)) p1 = gdiv(p1,p2);
if ((lgef(x)-3) & 2) p1 = gneg(p1);
return gerepileupto(av,p1);
--- 2071,2080 ----
{
case t_POL:
if (gcmp0(x)) return gzero;
!       if (v < 0)
!         p1 = subres(x, derivpol(x));
!       else
!         p1 = polresultant0(x, deriv(x,v), v, 0);
p2 = leading_term(x); if (!gcmp1(p2)) p1 = gdiv(p1,p2);
if ((lgef(x)-3) & 2) p1 = gneg(p1);
return gerepileupto(av,p1);
***************
*** 2067,2073 ****
return stoi(-4);

!       return discsr((GEN)x[1]);

case t_QFR: case t_QFI:
p1=sqri((GEN)x[2]);
--- 2083,2089 ----
return stoi(-4);

!       return poldisc0((GEN)x[1], v);

case t_QFR: case t_QFI:
p1=sqri((GEN)x[2]);
***************
*** 2076,2082 ****

case t_VEC: case t_COL: case t_MAT:
i=lg(x); z=cgetg(i,tx);
!       for (i--; i; i--) z[i]=(long)discsr((GEN)x[i]);
return z;
}
err(typeer,"discsr");
--- 2092,2098 ----

case t_VEC: case t_COL: case t_MAT:
i=lg(x); z=cgetg(i,tx);
!       for (i--; i; i--) z[i]=(long)poldisc0((GEN)x[i], v);
return z;
}
err(typeer,"discsr");
***************
*** 2084,2089 ****
--- 2100,2111 ----
}

GEN
+ discsr(GEN x)
+ {
+   return poldisc0(x, -1);
+ }
+
+ GEN
reduceddiscsmith(GEN pol)
{
long av=avma,tetpil,i,j,n;
*** src/headers/paridecl.h.orig	Tue Jan 19 18:25:58 1999
--- src/headers/paridecl.h	Thu Jan 21 15:17:03 1999
***************
*** 997,1002 ****
--- 997,1003 ----
GEN     nfisincl0(GEN a, GEN b,long flag);
GEN     nfisisom0(GEN a, GEN b,long flag);
GEN     nfiso(GEN a, GEN b);
+ GEN     poldisc0(GEN x, long v);
GEN     polfnf(GEN a, GEN t);
GEN     polresultant0(GEN x, GEN y,long v,long flag);
GEN     polsym(GEN x, long n);
*** src/language/init.c.orig	Fri Jan 15 13:59:44 1999
--- src/language/init.c	Thu Jan 21 15:16:44 1999
***************
*** 1495,1501 ****
{"polcompositum",25,(void*)polcompositum0,6,"GGD0,L,"},
{"polcyclo",99,(void*)cyclo,7,"LDn"},
{"poldegree",99,(void*)poldegree,7,"lGDn"},
! {"poldisc",18,(void*)discsr,7,"G"},
{"poldiscreduced",18,(void*)reduceddiscsmith,7,"G"},
{"polgalois",99,(void*)galois,6,"Gp"},
{"polinterpolate",31,(void*)polint,7,"GGDGD&"},
--- 1495,1501 ----
{"polcompositum",25,(void*)polcompositum0,6,"GGD0,L,"},
{"polcyclo",99,(void*)cyclo,7,"LDn"},
{"poldegree",99,(void*)poldegree,7,"lGDn"},
! {"poldisc",99,(void*)poldisc0,7,"GDn"},
{"poldiscreduced",18,(void*)reduceddiscsmith,7,"G"},
{"polgalois",99,(void*)galois,6,"Gp"},
{"polinterpolate",31,(void*)polint,7,"GGDGD&"},
*** src/language/helpmsg.c.orig	Tue Dec 15 16:30:18 1998
--- src/language/helpmsg.c	Thu Jan 21 15:18:10 1999
***************
*** 294,300 ****
"polcompositum(pol1,pol2,{flag=0}): vector of all possible compositums of the number fields defined by the polynomials pol1 and pol2. If (optional) flag is set (i.e non-null), output for each compositum, not only the compositum polynomial pol, but a vector [pol,al1,al2,k] where al1 (resp. al2) is a root of pol1 (resp. pol2) expressed as a polynomial modulo pol, and a small integer k such that al2+k*al1 is the chosen root of pol",
"polcyclo(n,{v=x}): n-th cyclotomic polynomial (in variable v)",
"poldegree(x,{v}): degree of the polynomial or rational function x with respect to main variable if v is omitted, with respect to v otherwise. Return -1 if x = 0, and 0 if it's a non-zero scalar",
!   "poldisc(x): discriminant of the polynomial x",
"poldiscreduced(f): vector of elementary divisors of Z[a]/f'(a)Z[a], where a is a root of the polynomial f",
"polgalois(x): Galois group of the polynomial x (see manual for group coding)",
"polinterpolate(xa,ya,{x}): polynomial interpolation at x according to data vectors xa, ya",
--- 294,300 ----
"polcompositum(pol1,pol2,{flag=0}): vector of all possible compositums of the number fields defined by the polynomials pol1 and pol2. If (optional) flag is set (i.e non-null), output for each compositum, not only the compositum polynomial pol, but a vector [pol,al1,al2,k] where al1 (resp. al2) is a root of pol1 (resp. pol2) expressed as a polynomial modulo pol, and a small integer k such that al2+k*al1 is the chosen root of pol",
"polcyclo(n,{v=x}): n-th cyclotomic polynomial (in variable v)",
"poldegree(x,{v}): degree of the polynomial or rational function x with respect to main variable if v is omitted, with respect to v otherwise. Return -1 if x = 0, and 0 if it's a non-zero scalar",
!   "poldisc(x,{v}): discriminant of the polynomial x, with respect to main variable if v is omitted, with respect to v otherwise",
"poldiscreduced(f): vector of elementary divisors of Z[a]/f'(a)Z[a], where a is a root of the polynomial f",
"polgalois(x): Galois group of the polynomial x (see manual for group coding)",
"polinterpolate(xa,ya,{x}): polynomial interpolation at x according to data vectors xa, ya",
--
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
--