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


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

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);
  
      case t_QUAD: case t_POLMOD:
!       return discsr((GEN)x[1]);
  
      case t_QFR: case t_QFI:
        p1=sqri((GEN)x[2]);
--- 2083,2089 ----
        return stoi(-4);
  
      case t_QUAD: case t_POLMOD:
!       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
--
PARI/GP Home Page: http://hasse.mathematik.-tu-muenchen.de/ntsw/pari/