Today I found the paper "Rise and Set Time of a Satellite about an Oblate Planet" (Pedro Ramón Escobal, 1963), which presents "a closed-form solution to the satellite visibility problem". It reduces to finding the roots to this trascendental equation, which only has the eccentric anomaly $E$ as independent variable:
$$ F \triangleq a (\cos{E} - e)\mathbf{P}\cdot\mathbf{Z} + a\sqrt{1 - e^2}\sin{E} \mathbf{Q}\cdot\mathbf{Z} - G = 0 $$
where $\mathbf{P}$ and $\mathbf{Q}$ are functions of $(\Omega, \omega, i)$, $\mathbf{Z}$ is the unit vector pointing towards the geodetic zenith as a function of $E$, and $G$ is a constant related to the site geodetic coordinates and the planet ellipsoid of choice.
I implemented all the equations, you can find a rough draft on this open pull request to orbit-predictor.
However, after reflecting on this a lot I don't see the big advantage of this method. Several times in the article, Escobal says that this equation gives a "better or improved estimate", but gives no hints about how to find the nearest root of the trascendental equation. Yes, there is an approximate method, but it doesn't work for mid inclinations or altitudes above LEO, so it's not useful for the general case.
I have read several articles that propose numerical methods, and they all cite Escobal "analytical method", without really clarifying whether it works or not for their use case. It seems to me that, even though everybody cites this seminal work, nobody has really tried to reproduce or implement it.
Am I missing something? Do people just implement some sort of "brute force" method with a resolution of tens of seconds to minutes and that's it?