Do not derive them by hand, use Maxima:
p: a*x^2+b*x=y;
s: solve([ev(p,x=x1,y=y1),ev(p,x=x2,y=y2)],[a,b])[1];
a < 0,s;
m: denom(lhs(ev(p,s,ratsimp)));
e: ev(p,s) * m, ratsimp;
to print as source code use
fortran(e);
It is also a good idea to use macro to avoid typos:
bool on_par(int a, int b, int c){ // [c] on the same parabola as [a] and [b]
# define _(a,b,c) x[c]*x[a]*(x[c]-x[a])*y[b]
return _(a,b,c) - _(b,a,c) == _(a,c,b);
# undef _
}