How to compute normal of surface from implicit equation for ray-marching? I want to represent ovaloids with ray-marching. I have ovaloids defined with an implicit equation. 

For instance, we can consider the equation of an ellipsoid given as: $$ \frac{x^2}{a^2} + \frac{y^2}{b^2} + \frac{z^2}{c^2} -1 = 0, $$ where $a, b, c$ are constant. 

Given this equation I need to caculate the normal for the shading of the suface. In an standard course of Curves and Surfaces I would take $f(x,y,z) = \frac{x^2}{a^2} + \frac{y^2}{b^2} + \frac{z^2}{c^2} -1 $ and compute $N(p) = \frac{\nabla f(p)}{|\nabla f(p)|}$ as the normal of the surface. 

Here I have to check that $0$ is a regular value of $f$ and that $f  $ is differentiable. However, the second check is not automate.


How is it done in real-world systems? 
Can you give an explanation of the methods used and ideally point me some implementations that exist?


Please, answer in full generality and not only considering the example given above.
The usual way of doing this with raymarching is to define your surface as a signed-distance field and use finite differences to get the gradient of the distance function at the point you’re sampling. In other words, if you have a function map(p) that returns the signed distance value at a point p, the normal at p is given by:

float epsilon = 0.001; // arbitrary — should be smaller than any surface detail in your distance function, but not so small as to get lost in float precision
float centerDistance = map(p);
float xDistance = map(p + float3(epsilon, 0, 0));
float yDistance = map(p + float3(0, epsilon, 0));
float zDistance = map(p + float3(0, 0, epsilon));
float3 normal = (float3(xDistance, yDistance, zDistance) - centerDistance) / epsilon;


Note that you can improve the accuracy by taking an extra 3 samples on the other “side” of the center position and subtracting those from the first three, then dividing by 2× the epsilon value, but if your distance function is complicated then that can be a lot of additional work for not much visible benefit.

The sample function you provide isn’t a signed-distance function as-is, and I don’t know of a straightforward way to convert it into one, but there’s an article by Íñigo Quílez (kind of a legend in this field) here with various functions you can make use of, including one for an ellipsoid. There’s also a lot of reference material available on Shadertoy—search for “raymarching” and you’ll find plenty of examples.