I'm trying to create an algorithm that could tell if a point lies on a Bézier curve (or very close to it) analytically.
I know that I can check if that's the case by solving $t$ for some ${x,y}$:
$x = a_x(1-t)^3 + 3b_xt(1-t)^2 + 3c_xt^2(1-t)+d_xt^3$
$y = a_y(1-t)^3 + 3b_yt(1-t)^2 + 3c_yt^2(1-t)+d_yt^3$
and if it has solutions in range $[0,1]$, then the point belongs to a Bézier.
$a_x, a_y, ... d_x, d_y $ are Bézier's start point, controls points and end point.
I know how to do that on a piece of paper if I substitute $a_x, a_y, ... d_x, d_y $ with some values, but I have no idea how to create an algorithm for doing that with my computer (mainly because I'm not able to combine polynomial terms since they are connected to $a_x, a_y, ... d_x, d_y $). Could someone please point me to a technique that can be used in such cases or to some source code that does what I want to achieve?
You have to eliminate ("Elimination" is an important keyword) variable $t$ in order to turn your parametric representation into a cartesian implicit equation $f(x,y)=0$, because it is easy to test if $f(x,y)=0$ or, more realisticaly, $ -\varepsilon < f(x,y) < + \varepsilon$ for some small $\varepsilon$.
How to do that ? The answer, adapted to computer programs, is : by using Sylvester's resultant. I refer you to the recent answer I made here in a similar case.
In order to simplify a little the computations, let us assume that the initial point ($t=0$) is the origin (it is always possible to do that), expanding formulas for $x$ and $y$ as polynomials in $t$ :
$$\begin{cases}U_xt^3+V_xt^2+W_xt-x=0\\ U_yt^3+V_yt^2+W_yt-y=0\end{cases}$$
$$\begin{cases}U_x&:=&3b_x-3c_x+d_x\\V_x&:=&3c_x-6b_x\\ W_x&:=& 3b_x\end{cases} \ \ \text{and} \ \ \begin{cases}U_y&:=&3b_y-3c_y+d_y\\V_y&:=&3c_y-6b_y\\ W_y&:=& 3b_y\end{cases}$$
one gets the following $6 \times 6$ "determinantal" expression for $f$:
$$f(x,y)=\left|\begin{array}{rrrrrrrr} U_x&V_x&W_x&-x&0&0\\ 0&U_x&V_x&W_x&-x&0\\ 0&0&U_x&V_x&W_x&-x\\ U_y&V_y&W_y&-y&0&0 \\ 0&U_y&V_y&W_y&-y&0\\ 0&0&U_y&V_y&W_y&-y \end{array}\right|=0.$$
You can leave it under this form if you use a software like Matlab. But if you have to do the operation $\approx10^9$ times, you need to expand it once for all (using a CAS).