Counting Distinct Real Roots of Low-Degree Polynomials

Question

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: [1], 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.

Answer 1

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[1]?(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[0])))|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[1] ?           // if p[1] is defined:
    ...                //   general case
  :                    // else:
    1 / !!+p - 1       //   return +Infinity if p[0] = 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[0])         //     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

Answer 2

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
          +         #  add a(n-2)
                    #  => 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

Answer 3

Wolfram Language (Mathematica), 25 32 31 27 bytes

+7 bytes for bugfix

-4 bytes thanks to @JungHwanMin

Length@#&@@RootIntervals@#&

Try it online!

Answer 4

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

Answer 5

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!

Answer 6

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.]...

Answer 7

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.

Answer 8

Jelly, 9 13 - 5 = 8 bytes

Monadic link expecting list of coefficients from lowest to highest.

Æ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.