beside (1 & (-1.01)) (square 1) (circle 1)
Looks OK. And now:
beside (1 & (-1)) (square 1) (circle 1)
Tracking this down a bit further, I've reduced it to a smaller test case:
p :: Trail R2
p = Trail [last $ trailSegments (circle 1)] False
main = defaultMain $ beside (1 & (-1)) (strokeT p) (square 1 # rotateBy (1/8))
The problem is that the envelope for the given cubic segment returns the wrong value for v = 1 & (-1). Ultimately, the envelope ends up solving the quadratic equation
v = 1 & (-1)
-3.3306690738754696e-16 x^2 - 3.3137084989847607 x + 1.6568542494923806 = 0
The correct solution is x = 1/2, but the quadForm solver does a really poor job here and returns [-9.949077574161476e15,-0.0]. The problem seems to be the extreme tininess of the x^2 coefficient which introduces all sorts of numerical instability.
Sigh. Maybe we ought to just switch to iterative rootfinding? What do "real" graphics programs (or other programs that have to solve "real world" quadratic or cubic equations) do?
I thought I looked at quadForm before to check this, but I guess I hadn't. There is a way to avoid the problems with small values given here http://www.mpi-hd.mpg.de/astrophysik/HEA/internal/Numerical_Recipes/f5-6.pdf.
My understanding is that for cubic and quartic you are better off with some iterative method (some are given in that book which I have around here somewhere).
Ah, thanks for the link! Works like a charm. That PDF also has some stuff about solving cubics, I wonder if we should be using that too. Though I've never (yet) noticed any problems with cubForm.
Make quadForm solver more numerically stable (fixes #46)
Using the formulas given in
note #46 fix in CHANGES [ci skip]