New plane solver

[unpairedbracket]

Opening as a draft because it still needs a docs pass, benchmarking to ensure it's actually faster/scales better than the existing algorithm and removing the old code if so.

Objective

The current plane solver is O(n^3) in the number of collision planes. It's not likely to break the bank but it'd be preferable to bring this down to O(n^2) or lower

Solution

• I figured out a GJK-like algorithm for solving the underlying quadratic programming problem, and implemented it. Steps of the algorithm are O(n), and it should terminate in at most n steps, though I've never seen the step number exceed 4. I've tried to construct sets of normals that would take 5+ steps to converge and the arrangements of planes it would take seem very specific. Blog post fully explaining the details is in progress

Testing

• Tested so far in the 3d move-and-slide example, still need to test in 2d as well
• I'm intending to benchmark the new method against the existing one for some arrangements of maybe 1 to 10 normals? to try and establish scalings empirically
• I think there are a few more micro-optimisations I could squeeze out, like imposing a specific ordering of the normals in SimplicialCone::Wedge(n1, n2) relative to the new search direction, which would save a few comparisons and dot products in the relevant branch of SimplicialCone::project_point.
• The SimplicialCone enum exists to keep the algorithm allocation-free (as opposed to keeping the current simplex points in a vec. Alternatively the first couple of iterations of the loop could be manually unrolled, which would get rid of the match and might speed things up a bit

[unpairedbracket]

Ok, I've managed to get a decent benchmark set up for this. Performance of both methods vary based on where the velocity is relative to the normal cone (the old one can early-out if all constraints are fine, the new one takes varying numbers of steps depending on whether the initial velocity is allowed, or the end result is constrained by one plane or two). Because of that I came up with a thing for sampling an endless sequence of points on the sphere evenly (based on a quasirandom sequence) to use a different velocity direction for each iteration. Uniform sampling around the whole sphere of directions might not be the most realistic distribution of directions for real-world use, but it's better than nothing.

Scaling for the new method looks better than the old one on a linear scale:

On log-log it appears that the old one is actually closer to quadratic than cubic (the brighter lines without markers are power laws, n for red and n^2 for green), and the new one is roughly linear (in case anyone reading this is unfamiliar with log-log plots, power laws with exponent α on linear plots translate to lines with gradient α on log-log plots: f(x) = x^α -> log f(x) = α log x):