# Counting Distinct Real Roots of Low-Degree Polynomials

Challenge: I want to know about the real roots of polynomials. As a pure mathematician, I care about the existence of such roots, rather than their numeric values.

The challenge is to write the shortest program that takes a polynomial, of degree at most 4, and simply returns how many distinct real roots said polynomial has. Polynomials with degree 4 or less have the unique property that there exist closed forms (such as the quadratic formula), which give all their roots. You can google these forms or find some useful related information in the appendix.

Input: the coefficients of a polynomial. For simplicity, we shall only care about polynomials with integer coefficients.

You may input this as a list, get the coefficients one at a time, or use any other reasonable method. You may, for example, require polynomials of degree d to inputted as lists of length d+1.

You should specify how you convert a polynomial of degree at most 4 into a valid input.

Output: the number of distinct real roots of said polynomial. (meaning roots with multiplicity are only counted once)

You must output one of the integers 0,1,2,3,4 for valid polynomials, and trailing spaces are completely fine. (the special case of the polynomial, $$\P(x) = 0\$$, is discussed in the scoring section)

Examples: in these examples, we represent a polynomial, $$\P\$$, as a list L, where L[i] contains the coefficient of $$\x^i\$$ in $$\P\$$. (with index starting at 0)

$$\P(x) = 1\$$, input: , output: 0

$$\P(x) = 1+2x\$$, input: [1,2], output: 1

$$\P(x) = 1+x^2\$$, input: [1,0,1], output: 0

$$\P(x) = 1+x+x^2+x^3 = (x+1)(1+x^2)\$$, input: [1,1,1,1], output: 1

$$\P(x) = 2+3x+x^2 = (x+2)(x+1)\$$, input: [2,3,1], output: 2

$$\P(x) = 1-2x+x^2 = (x-1)^2\$$, input: [1,-2,1], output: 1

Scoring:

-5 bytes if the polynomial, $$\P(x) = 0\$$, outputs a representation of infinity, such as: the infinity float in python, the Unicode symbol ∞, or a string that can be mathematically interpreted to evaluate to infinity, like "1/0". (otherwise, you don't need to handle this case)

Otherwise, shortest code wins. (however I am personally quite interested in seeing answers which don't rely on built-in root finders)

Appendix: Closed Forms

Starting with the basics:

Degree 0: $$\P(x) = a\$$, roots: none, (or all reals, if a = 0)

Degree 1: $$\P(x) = ax+b\$$, roots: $$\x=-b/a\$$

Degree 2: $$\P(x) = ax^2+bx+c\$$, roots: $$\x = \frac{-b \pm \sqrt{b^2-4ac}}{2a}\$$

Then, we have the harder ones. For these, it becomes quite verbose to state the closed forms. I've paraphrased some ways you can deduce the number of roots by considering discriminants, you may also seek to find the closed forms of the roots.

Degree 3: From Wikipedia. $$\ P(x) = ax^3+bx^2+cx+d\$$. We define the discriminant as $$\\Delta = 18abcd -4b^3d+b^2d^2-4ac^3-27a^3d^2\$$.

If $$\\Delta > 0\$$ then there are 3 distinct roots, if $$\\Delta = 0\$$, then there are two distinct real roots, otherwise then $$\ \Delta < 0\$$ and there is only one real root.

Degree 4: From Wikipedia.$$\ P(x) = ax^4+bx^3+cx^2+dx+e\$$. We define the discriminant as \begin{align} \Delta\ =\ &256 a^3 e^3 - 192 a^2 b d e^2 - 128 a^2 c^2 e^2 + 144 a^2 c d^2 e - 27 a^2 d^4 \\ &+ 144 a b^2 c e^2 - 6 a b^2 d^2 e - 80 a b c^2 d e + 18 a b c d^3 + 16 a c^4 e \\ &- 4 a c^3 d^2 - 27 b^4 e^2 + 18 b^3 c d e - 4 b^3 d^3 - 4 b^2 c^3 e + b^2 c^2 d^2 \end{align}

If $$\\Delta > 0\$$ then there are two distinct real roots. Otherwise things get more complicated.

Defining $$\P = 8ac - 3b^2\$$, $$\D = 64 a^3 e - 16 a^2 c^2 + 16 a b^2 c - 16 a^2 bd - 3 b^4\$$, then if $$\ \Delta < 0 \$$ and $$\ P < 0\$$ and $$\ D < 0 \$$ then there are 4 distinct real roots. Otherwise, if $$\ \Delta < 0\$$, there are 2 distinct real roots.

Finally, if $$\\Delta = 0\$$, we define $$\\Delta_0=c^2-3bd+12ae\$$ and $$\R= b^3+8da^2-4abc\$$.

If $$\D =0\$$, then:

• $$\ P < 0\$$ implies two distinct real roots
• $$\ P > 0\$$ implies zero distinct real roots
• $$\P = 0\$$ implies four distinct real roots

Else, if $$\\Delta_0 = 0\$$, then there are two distinct real roots.

Else, if $$\P<0\$$ and $$\D <0\$$ there are three distinct real roots.

Else, there are one real root.

• @Anush The polynomial $x^5-x+1$ is an example of a polynomial with a root which cannot be expressed with radicals. (instead, we need hypergeometric functions) Since the roots no longer have a nice expression, we can't get an easy closed form, making considering higher degrees much less simple. Feb 26, 2020 at 9:23
• @Neil How is trivial?
– user9207
Feb 26, 2020 at 17:05
• Feb 26, 2020 at 18:36
• For the record, "count the number of roots in an interval" is a task that does not require being able to write down the roots; Sturm's theorem is a bit technical to explain, but it is a simple algorithm to count the number of roots of a polynomial over an interval (or over the whole real line) (and this is used in the answers) Feb 26, 2020 at 23:31
• I'm wondering if a golfy approach here would be to evaluate the polynomial at a huge number of equally spaced points, and count the intervals where it crosses zero or nearly touches it. This would require a mesh size small enough to separate any two roots, and an epsilon where anything approaching that close to an axis must actually be multiple root touching it. I suspect something suitably small in the polynomial degree and coefficients will suffice, but don't now how to prove it.
– xnor
Feb 27, 2020 at 4:12

# JavaScript (Node.js),  285 ... 248  247 - 5 = 242 bytes

Fixed a bug and saved a couple of bytes by borrowing @Grimmy's method for the final count of the real roots

Expects the coefficients from highest to lowest.

This is an implementation of Sturm's theorem, and therefore works for any degree, which might be a bit overkill for this challenge.

p=>1/p?(h=(A,[[c,...q],[...p]]=A,n=q.length)=>n?h([(g=i=>1/p[(j=i)+n]?g(i+1,q.map(v=>p[++j]-=v*p[i],p[i]/=c)):p.slice(-n).map(v=>-v))(0),...A]):A)([p.flatMap(v=>1/p[++i]?(p.length-i)*v:[],i=k=s=0),p]).map(p=>s+=Math.sign(k*(k=p)))|s:1/!!+p-1


Try it online!

### How?

We first test whether the degree is $$\0\$$, which is an edge case:

p =>                   // p[] = polynomial coefficients
1 / p ?           // if p is defined:
...                //   general case
:                    // else:
1 / !!+p - 1       //   return +Infinity if p = 0, or 0 otherwise


For the general case, we use the recursive function $$\h\$$ to build the Sturm sequence of polynomials in reverse order $$\[p_n,...,p_1,p_0]\$$, where $$\p_0\$$ is the input and $$\p_n\$$ is of degree $$\0\$$.

The initial call to $$\h\$$ is performed with $$\[q,p]\$$ where $$\q\$$ is the derivative of $$\p\$$ computed as follows:

p.flatMap(v =>         // for each coefficient v in p[]:
1 / p[++i] ?         //   increment i; if p[i] is defined:
(p.length - i) * v //     multiply v by the exponent at this position - 1
:                    //   else:
[],                //     this is the constant term: discard it
i = k = s = 0        //   start with i = 0 (k and s are used later)
)                      // end of flatMap()


Below is the definition of $$\h\$$:

h = (                  // h is a recursive function taking:
A,                   //   A[] = current Sturm sequence, which is split into:
[[c, ...q], [...p]]  //     c = 1st coefficient of the divisor
= A,                 //     q[] = other coefficients of the divisor
//     p[] = copy of the dividend
n = q.length         //   n = number of coefficients in q[]
) =>                   //
n ?                  // if n is not equal to 0:
h([g(0), ...A])    //   preprend a new polynomial to A[], using g (see below)
:                    // else:
A                  //   end of recursion: return A[]


Where $$\g\$$ computes the remainder of $$\p/q\$$, multiplied by $$\-1\$$:

g = i =>               // g is a recursive function taking a counter i
1 / p[(j = i) + n] ? // copy i into j; if p[j + n] is defined:
g(                 //   do a recursive call:
i + 1,           //     with i + 1
q.map(v =>       //     for each coefficient v in q[]:
p[++j] -=      //       increment j; subtract from p[j]:
v * p[i],    //         v multiplied by p[i]
p[i] /= c      //       start by normalizing p[i], using c
)                //     end of map()
)                  //   end of recursive call
:                    // else:
p.slice(-n)        //   end of recursion: isolate the remainder
.map(v => -v)      //   and invert all signs


Finally, we compute the number of real roots by applying the following code to the Sturm sequence:

.map(p =>              // for each polynomial p in the sequence:
s += Math.sign(      //   add to s the sign of:
k *                //     multiply the previous leading coefficient
(k = p)         //     by the leading coefficient of p (and update k)
)                    //   NB: k and s are initialized to 0 earlier in the code
) | s                  // end of map(); yield s

• When roots fall together, suchs as in the last example ([1,-2,1]) in the question, your code counts them as separate solutions. So [1,-2,1] returns 2, but should return 1. Although +1 for implementing Sturm theorem! Feb 26, 2020 at 15:23
• @agtoever This is now fixed. Feb 26, 2020 at 16:33

# 05AB1E, 39 bytes

gŠ¦ƶ‚íλ£¬/D₂0š*₂-2£δ*˜0š2ôO+(¦¦}€нüP.±O


Try it online!

05AB1E doesn't have a built-in to get the roots of a polynomial, so we use Sturm's theorem

g                   # length of the input (degree + 1)
Š                  # push two more copies of the input
¦                 # drop the first coefficient
ƶ                # multiply each coefficient by its 1-based index
# => this is the derivative of the input
‚               # pair the input with its derivative
í              # reverse both, so the constant coeff is now last

λ£            }     # recursively compute the Sturm's sequence:
# (initially, a(n-2) = input and a(n-1) = derivative)
¬/                #  divide a(n-1) by its first element
D               #  duplicate
₂              #  push a(n-2)
0š            #  prepend a 0
*           #  multiply
₂-         #  subtract a(n-2)
2£       #  keep the first two coefficients
#  => this is the polynomial division -a(n-2) / a(n-1)
δ*                #  double-vectorized multiplication
˜               #  flatten
0š             #  prepend a 0
2ô           #  split in groups of 2
O          #  sum each group
#  => this is the remainder
(        #  negate
¦¦      #  drop the two leading zeroes

€н                  # keep only the highest-order coefficient of each polynom
üP                # pairwise product
.±              # sign of each
O             # sum

• Amazing! Could you explain how it works? Feb 26, 2020 at 14:30
• @ZacharyHunter explanation added. Feb 26, 2020 at 14:46

# Wolfram Language (Mathematica), 253231 27 bytes

+7 bytes for bugfix

-4 bytes thanks to @JungHwanMin

Length@#&@@RootIntervals@#&


Try it online!

• @JungHwanMin this was posted before as a comment.. and still fails for 2 + 3x + x^2.. I'll leave this comment up this time so nobody posts it again Feb 26, 2020 at 14:32
• Ah, right, it would be a vertical vector. This works, but is same length. Feb 26, 2020 at 14:36
• 27 bytes, for real this time Feb 26, 2020 at 14:44
• Tr fails for cubic Feb 26, 2020 at 14:51
• Ok yeah thanks good point @JungHwanMin Feb 26, 2020 at 14:52

# Jelly, 12 bytes - 5 = 7

ÆrQṠṠƑ€SµİẸ?


A monadic link accepting a list of coefficients from highest degree to lowest which yields a number or the list [inf] if all the integers given were zero. As a full program prints the number or inf.

Try it online!

### How?

ÆrQṠṠƑ€SµİẸ? - Link: list of coefficients
? - if...
Ẹ  - ...condition: any?
µ    - ...then: the monadic chain:
Ær           -   roots
Q          -   de-duplicate
Ṡ         -   sign if real; complex conjugate if complex (vectorises)
€      -   for each:
Ƒ       -     is invariant under?:
Ṡ        -       sign if real; complex conjugate if complex
S     -   sum
İ   - ...else: inverse


# Pari/GP, 8 bytes

Pari/GP has a built-in for this.

It takes an polynomial as input, e.g., 1-2*x+x^2.

polsturm


Try it online!

# Python 3, 91 bytes

lambda c:sum(f.real==f for f in set(polyroots(c)))
from numpy.polynomial.polynomial import*


Try it online!

Had to include set() to take care of "double" roots. polyroots([1, -2, 1]) gives [1., 1.]...

# R, 44 bytes

sum(!Im(unique(zapsmall(polyroot(scan())))))


Try it online!

zapsmall is necessary to deal with numerical imprecision, as otherwise e.g. 2 3 1 test case would fail.

# Jelly, 9 13 - 5 = 8 bytes

ÆrÆiṪÐḟQLµİẸ?


You can try it online! How it works:

Ær        Find the roots of the polynomial
Æi       and split each root into [real part, imaginary part]
Ðḟ   Filter the list, removing elements
Ṫ      where the imaginary part is Truthy (i.e. other than 0)
Q  Remove duplicates
L  and return the length of that.


This is the main part. Then µİẸ? checks if there is any non-zero coefficient. If not, returns infinity. If there is, just run the previous program I described.

Borrowed the final four bytes µİẸ? from Jonathan Allan's Jelly answer, increasing my count by 4 but gaining the bonus -5, for a total of -1.