Suppose one wishes to test if a given algebraic surface f(x,y,z,w) = 0 in projective 3 space has singular points, that is df/dx = df/dy = df/dz = 0, and one also wishes to calculate these singular points explicitly. Is there a general algorithm for doing this, and does there exist an implementation of this algorithm in a modern programming language?
I have written my own script in MatLab which does find singular points, but it doesn't find all of them. It takes the naive approach and simply calculates the partial derivatives and attempts to solve for the system of partial derivative equations and f = 0. For a quintic polynomial which is known to have 31 singularities (the togliatti quintic), my algorithm only finds 20 singularities.
I've also looked into the programming language Magma, which has a function "SingularPoints(C) : Sch -> SetIndx" which claims to return "The singular points of the curve C which are defined over the base field of C." But beyond this, the documentation is a bit sparse, and doesn't say if this function works for curves of more than two real variables. Here's the link to the documentation:
http://magma.maths.usyd.edu.au/magma/handbook/text/1415#15982
I found the above documentation after looking at the following paper on arxiv, in the section on togliatti quintic equations:
https://arxiv.org/pdf/2206.05492 (page 109)
