Karim Belabas on Sat, 03 May 2008 12:46:02 +0200


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

Re: poldisc shortcut known but not used


* Phil Carmody [2008-04-16 11:15]:
> It appears that this is well known:
> 
>> There's more than one way to calculate a polynomial discriminant.
>> It  turns out that when the coefficients are ugly multivariate (but
>> exact)  expressions, one of the polresultant() algorithms is much
>> faster.   Says so right in the manual!  For a monic polynomial px of
>> degree m,  poldisc(px) and polresultant(px,deriv(px)) differ
>> formally by a factor  of (-1)^(m(m-1)/2).  Dragging out the old
>> chestnut polynomial shows  that the manual wasn't kidding about the
>> other algorithm being faster:
>>
>> ? px=x^10 + (2*B*u + (2*A + 4*B))*x^9 + (9*A*u + (9*A +
>> (9*B - 45)))*x^8 + ((-24*A - 120)*u + (-24*A + (-24*B - 120)))*x^7
>> + ((-42*B + 210)*u + (-42*A + (-84*B + 210)))*x^6 + 252*x^5 + (42*B*u
>> + (42*A + 84*B))*x^4 + (24*A*u + (24*A + (24*B - 120)))*x^3 + ((-9*A
>> - 45)*u + (-9*A + (-9*B - 45)))*x^2 + ((-2*B + 10)*u + (-2*A + (-4*B
>> + 10)))*x + 1;
>> ? disc1=poldisc(px);
>> time = 46,859 ms.
>> ? disc2=-polresultant(px,deriv(px),x,2);\\Factor of -1 is (-1)^(10*9/2)
>> time = 10,475 ms.

( You may use px' instead of deriv(px). )

> Is there any reason why this case can't be detected automatically, and
> the shortcut used internally?

Had never gotten around to understanding my own implementation
(never include algorithms you don't understand:-).

Did that yesterday.

In the svn repository, polresultant(), poldisc() [ and all functions
that call them, e.g. charpoly(t_POLMOD) ] now use a combination of
1) modular algorithms for univariate polynomials over Z/Q,
2) Ducos's subresultant for other exact polynomials (with major
improvements for "simple" polynomials).
3) and Sylvester's matrix when nothing else works

4) our previous subresultant has almost disappeared, but survives in
functions like gcd() / bezout() / bezoutres(), for the time being.
Not for long in the case of gcd(); the other two are more annoying and
I probably won't adapt 2) to this case.

Since 1) takes over the "simple base ring" case, and other "simple" cases
are now treated adequately by 2), the latter now becomes the default.
The threshold with 4) was easy to determine: after a few simple
improvements, "ugly multivariate" turned out to encompass all
multivariate polynomials, as well as the univariate ones of large degree...

polresultant(P,Q,,2) is kept for backward compatibility, but gives you
the default.

There's a new 'test-resultant' bench, exhibiting cases where the new
subresultant is clearly superior to the old one.

I'm interested in all regression cases :-).

Cheers,

    K.B.

P.S: Someday, someone may sit down and implement proper
evaluation/interpolation schemes for general multivariate polynomials,
but it's definitely not a priority... the case that matters to PARI is
Z[X,Y], and it is already treated adequately.

--
Karim Belabas, IMB (UMR 5251)  Tel: (+33) (0)5 40 00 26 17
Universite Bordeaux 1          Fax: (+33) (0)5 40 00 69 50
351, cours de la Liberation    http://www.math.u-bordeaux.fr/~belabas/
F-33405 Talence (France)       http://pari.math.u-bordeaux.fr/  [PARI/GP]
`