Colleagues, g <- function(x) ( x^11 + 1000*x^10 + 500 *x^9 + 1 ) coeffs <- c(1, rep(0, 8), 500, 1000, 1) roots <- polyroot(coeffs) Output [1] 0.25770068+3.958197e-01i [2] -0.34615184+3.782848e-01i [3] -0.04089779-4.838134e-01i [4] 0.44124314-1.517731e-01i [5] -0.04089779+4.838134e-01i [6] -0.56201931-1.282822e-01i [7] -0.34615184-3.782848e-01i [8] 0.44124314+1.517731e-01i [9] -0.56201931+1.282822e-01i [10] 0.25770068-3.958197e-01i [11] -999.49974975+1.110223e-16i [11] -999.49974975+1.110223e-16i makes no sense since f>0 for all x Why does polyroot do this? Thomas Subia
Dear Thomas, How can f be strictly positive if the leading term is an odd power of x? For large negative x, x^11 will dominate, and f should become quite negative. This argument is also supported by a graphical plot of f. Assuming that f can indeed be negative, my suspicion is that the 11th root is still a reasonable approximation to the actual root, but the gradient in this region is high and so f(x) evaluated numerically is too sensitive to small errors (that is why f(roots[[11]]) evaluates to a seemingly huge number). Best, Lukas On Thu, Oct 2, 2025 at 6:14?PM tgs77m--- via R-help <r-help at r-project.org> wrote:> Colleagues, > > g <- function(x) ( x^11 + 1000*x^10 + 500 *x^9 + 1 ) > coeffs <- c(1, rep(0, 8), 500, 1000, 1) > roots <- polyroot(coeffs) > > Output > > [1] 0.25770068+3.958197e-01i > [2] -0.34615184+3.782848e-01i > [3] -0.04089779-4.838134e-01i > [4] 0.44124314-1.517731e-01i > [5] -0.04089779+4.838134e-01i > [6] -0.56201931-1.282822e-01i > [7] -0.34615184-3.782848e-01i > [8] 0.44124314+1.517731e-01i > [9] -0.56201931+1.282822e-01i > [10] 0.25770068-3.958197e-01i > [11] -999.49974975+1.110223e-16i > > [11] -999.49974975+1.110223e-16i makes no sense since f>0 for all x > > Why does polyroot do this? > > Thomas Subia > > ______________________________________________ > R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide > https://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. >[[alternative HTML version deleted]]
I haven't thought about this at all, but for what it's worth, sympy
gives the same answer (minus the tiny imaginary part):
from sympy import nroots
from sympy.abc import x
g = x**11 + 1000*x**10 + 500 *x**9 + 1
nroots(g)[0]
-999.499749749687
The documentation for
https://docs.sympy.org/latest/guides/solving/find-roots-polynomial.html
mentions that
nroots() can fail sometimes for polynomials that are numerically ill
conditioned, for example Wilkinson?s polynomial. Using RootOf() and
evalf() as described in Numerically Evaluate CRootOf Roots can avoid
both ill-conditioning and returning spurious complex parts because it
uses a more exact, but much slower, numerical algorithm based on
isolating intervals.
FWIW the slow/more exact method described here still gives the large
root.
from sympy import solve
g_solved = solve(g, x, dict=True);
for root in g_solved:
print(root[x].n(10))
Given the convergence of different methods, this is very unlikely to
be a bug in R or the underlying algorithm. Either there's a thinko and
-999.499... is actually a root, or this is a badly conditioned problem
for which it's extremely hard to get accurate numerical solutions ...
On 2025-10-02 12:14 p.m., tgs77m--- via R-help wrote:> Colleagues,
>
> g <- function(x) ( x^11 + 1000*x^10 + 500 *x^9 + 1 )
> coeffs <- c(1, rep(0, 8), 500, 1000, 1)
> roots <- polyroot(coeffs)
>
> Output
>
> [1] 0.25770068+3.958197e-01i
> [2] -0.34615184+3.782848e-01i
> [3] -0.04089779-4.838134e-01i
> [4] 0.44124314-1.517731e-01i
> [5] -0.04089779+4.838134e-01i
> [6] -0.56201931-1.282822e-01i
> [7] -0.34615184-3.782848e-01i
> [8] 0.44124314+1.517731e-01i
> [9] -0.56201931+1.282822e-01i
> [10] 0.25770068-3.958197e-01i
> [11] -999.49974975+1.110223e-16i
>
> [11] -999.49974975+1.110223e-16i makes no sense since f>0 for all x
>
> Why does polyroot do this?
>
> Thomas Subia
>
> ______________________________________________
> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
https://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
--
Dr. Benjamin Bolker
Professor, Mathematics & Statistics and Biology, McMaster University
Director, School of Computational Science and Engineering
> E-mail is sent at my convenience; I don't expect replies outside of
working hours.
How can an odd degree polynomial be positive for all x? Deepayan On Thu, 2 Oct, 2025, 6:14 pm tgs77m--- via R-help, <r-help at r-project.org> wrote:> Colleagues, > > g <- function(x) ( x^11 + 1000*x^10 + 500 *x^9 + 1 ) > coeffs <- c(1, rep(0, 8), 500, 1000, 1) > roots <- polyroot(coeffs) > > Output > > [1] 0.25770068+3.958197e-01i > [2] -0.34615184+3.782848e-01i > [3] -0.04089779-4.838134e-01i > [4] 0.44124314-1.517731e-01i > [5] -0.04089779+4.838134e-01i > [6] -0.56201931-1.282822e-01i > [7] -0.34615184-3.782848e-01i > [8] 0.44124314+1.517731e-01i > [9] -0.56201931+1.282822e-01i > [10] 0.25770068-3.958197e-01i > [11] -999.49974975+1.110223e-16i > > [11] -999.49974975+1.110223e-16i makes no sense since f>0 for all x > > Why does polyroot do this? > > Thomas Subia > > ______________________________________________ > R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide > https://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. >[[alternative HTML version deleted]]