Karim Belabas on Fri, 06 Nov 2015 13:36:46 +0100


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

Re: real polynomials


* John Cremona [2015-11-06 12:13]:
> I am using the library functions sturm() and sturmpart() to get the
> number of real roots of a real polynomial in the whole real line or
> just the positive or negative half-lines.  But occasionally it fails
> with
> 
>  ***   domain error in polsturm: issquarefree(pol) = 0
> 
> which I expect is only happening if the discriminant is exactly 0.

No, it can occur due to round-off errors whenever the input is inexact
and two roots are "close enough".

The documentation (??polsturm) specifies that the input must be
squarefree. (And, in 2.8.*, that it better be exact...)

> Is there a reliable way of taking the squarefreepart, or
> the radical, or a real polynomial?

Not if the input is inexact.

> It is in fact true that in my program the polynomial coefficients are
> always rational with denominator a power of 2, so will be known
> exactly.

In that case, the above problems don't occur. Do you need the number of
real roots including multiplicity ? If not, P /= gcd(P,P') is enough.

> I would like to avoid having to keep track of the
> denominators and using (scaled) integer coefficients.

Just use rational coefficients, sturm/sturmpart will cope.

If you're using 2.8.* and library mode, you should try 
ZX_sturmpart / ZX_sturm  (sturmpart itself is now deprecated).

Cheers,

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