Number of normals to a parabola from a given point

Claude's technique can be extended slightly to find the the points and normal lines given a particular point off of the parabola. Using this, I obtained an animation of some of the normals.

If $(p_x,p_y)$ are the coordinates of a point off the parabola and we'd like the line from $(p_x,p_y)$ to the point $(x,x^2)$ on the parabola to be normal to the parabola, then we need:

$$ \begin{pmatrix}p_x \\ p_y\end{pmatrix} = \begin{pmatrix}x \\ x^2\end{pmatrix} + t \begin{pmatrix}-2x \\ 1\end{pmatrix}. $$

This defines a pair of equation from which we can eliminate $x$ to obtain a cubic in $t$:

$$ 0 = 4t^3 - 4(p_y+1)t^2 + (4 p_y+1) t + p_x^2 - p_y = 0 $$

The discriminant of a cubic equation tells how many roots there are:

$$ 0 = ax^3+bx^2+cx+d\\ \Delta=18abcd-4b^{3}d+b^{2}c^{2}-4ac^{3}-27a^{2}d^{2} $$

Then $\Delta > 0$ corresponds to 3 real roots, $\Delta < 0$ to 1 real root and 2 complex roots, and $\Delta = 0$ to multiple root with all roots real.

I don't know for the case $\Delta = 0$ whether there are two distinct roots (one with multiplicity 2) or one root with multiplicity 3. Maybe both could occur.

I plotted the sign of the discriminant, with the parabola and a grid of unit size overlayed. The grey regions correspond to positive discriminant (3 real roots), white regions to negative discriminant (1 real root).

Here is the GLSL source code, for use with Fragmentarium:

#include "Progressive2D.frag"

vec3 color(vec2 p)
{
  float s = length(vec4(dFdx(p), dFdy(p)));
  if (p.y - s < p.x * p.x && p.x * p.x < p.y + s) return vec3(0.0);
  if (abs(mod(p.x + 0.5, 1.0) -0.5) < 0.5 * s) return vec3(0.5);
  if (abs(mod(p.y + 0.5, 1.0) -0.5) < 0.5 * s) return vec3(0.5);
  if ((abs(p.x) - s) * (abs(p.x) - s) < p.y &&
      p.y < (abs(p.x) + s) * (abs(p.x) + s)) return vec3(0.0);
  float a = 4.0;
  float b = -4.0 * (p.y + 1.0);
  float c = 4.0 * p.y + 1.0;
  float d = p.x * p.x - p.y;
  float discriminant
    = 18.0 * a * b * c * d
    - 4.0 * b * b * b * d
    + b * b * c * c
    - 4.0 * a * c * c * c
    - 27.0 * a * a * d * d;
  if (discriminant > 0.0) return vec3(0.7);
  if (discriminant < 0.0) return vec3(1.0);
                          return vec3(0.2);
}