Pre-compute solver data for each contact: velocity bias (restitution) and tangent direction.
Must be called ONCE per frame before warm starting and solver iterations.
This follows the Box2D approach: bias and tangent are fixed for the entire solve.
One iteration of velocity-level contact resolution with accumulated impulse clamping.
reverse: if true, iterate contacts in reverse order (for alternating Gauss-Seidel).
Pre-apply cached impulses from previous frame (warm starting).
This gives the iterative solver a head start, dramatically improving convergence for stacks.