Improve numerical stability of ground plane distance estimation

[Ralith]

At sufficiently high elevations with respect to the ground plane, glitchy floating islands appear. This is likely due to exponential growth of the magnitude of the normal vector encoding the distance to the ground plane. We could mitigate the representation problem by throwing additional precision at the problem in the form of an explicit exponent: rather than the vector v, store the pair (v/(2^c), c). Further investigation (@MagmaMcFry?) is necessary to determine the details of computing distance/direction of the plane from this representation without running into the same issue again.

[MagmaMcFry]

Suppose you have a pair (x,c), where x is a vec4 and c is a nonnegative number, representing the unit normal x*2^c of a plane H, and we want to compute the distance d between H and a point p. We know that sinh(d) = mip(x*2^c, p) = m*2^c, where m := mip(x,p).

Case 0: m is zero. Then d = 0.

Case 1: m is positive. Compute k := c+log2(m) = log2(sinh(d)). If k is relatively small, then you can just compute d = asinh(2^k) directly. If k is relatively large (at least 26 or so), then up to double precision, sinh(d) = exp(d)/2, so log2(sinh(d)) = log2(exp(d)/2) = d*log2(e) - 1 = k, so you can compute d = (k+1)*ln(2).

Case 2: m is negative. In this case, negate m, apply the computation in case 1, then negate d.

The only possibility for numerical instability is when p is very far away from the current frame of reference (and then only when p appears to be very close to H in the local Klein model). We won't run into this case as long as p is close to the current node.

---

To compute the perpendicular direction y to H at a given point p, you don't even need c. Just project x to the tangent space at p using y := x - p*mip(x,p)/mip(p,p) , then renormalize y.

[MagmaMcFry]

As a unit test for this system, I suggest starting with the extended normal (x=(1,0,0,0), c=0) and repeatedly boosting it away from the origin using a boost by distance 10 in negative X direction. After the nth boost, the distance between the normal and the point p = (0,0,0,1) should be 10*n. Repeat this 200 times to make sure that it works even if unextended normals would overflow double-precision floats.

[patowen]

Currently, I believe the main source of this numerical precision issue is in Plane's multiplication implementation. It runs a lorentz_normalize, which is numerically unstable when the plane is far from the origin.

This isn't necessary because the plane's normal is multiplied by an isometry. I believe simply getting rid of the call to lorentz_normalize fixes this issue, but if you explore further, you will run into the maximum f64 value, in which case we will need to use the method MagmaMcFry mentioned above.