How to compute normal of surface from implicit equation for ray-marching?
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.