It won't let me add this comment to the unmodified line below, but in playing with this PR I do infrequently run into a domain error on arcsin due to floating point errors. A little hard to reproduce reliably, but it does pop up. Clipping the inside value to [-1, 1] should fix this.

/mnt/c/Users/Scott/Documents/Documents/Coding/matplotlib/lib/mpl_toolkits/mplot3d/axes3d.py:4028: RuntimeWarning: invalid value encountered in arcsin
  elev = np.arcsin( 2*( qw*qy+qz*qx)/(qw*qw+qx*qx+qy*qy+qz*qz))  # noqa E201
posx and posy should be finite values

/mnt/c/Users/Scott/Documents/Documents/Coding/matplotlib/lib/mpl_toolkits/mplot3d/axes3d.py:4028: RuntimeWarning: invalid value encountered in arcsin
  elev = np.arcsin( 2*( qw*qy+qz*qx)/(qw*qw+qx*qx+qy*qy+qz*qz))  # noqa E201
posx and posy should be finite values

I'm just so curious about the values of qw, qx, qy, and qz that cause this. I can understand the posx and posy should be finite values complaint, as a consequence of the arcsin argument getting out of range. But the invalid value encountered in arcsin is still mysterious...
Could you try for me sometime replacing the offending statement with the following:

try:
    elev = np.arcsin( 2*( qw*qy+qz*qx)/(qw*qw+qx*qx+qy*qy+qz*qz))  # noqa E201
except:
    print(repr(qw), repr(qx), repr(qy), repr(qz))
    print(repr(qw*qy), repr(qz*qx), repr(qw*qw), repr(qx*qx), repr(qy*qy), repr(qz*qz) )
    print(repr( 2*( qw*qy+qz*qx) ), repr(qw*qw+qx*qx+qy*qy+qz*qz) )
    print(repr( 2*( qw*qy+qz*qx)/(qw*qw+qx*qx+qy*qy+qz*qz) ))

(I thought the error merits an attempt at diagnosing; not that I'm opposed to the np.clip(), I just thought it would be good to understand what is going on, how this comes about.)

I just spent a few minutes trying to reproduce but couldn't - it's a tricky edge case to trigger. When I ran across it before, the value inside arcsin was only 1e-16 (or thereabouts) larger than 1, so it's just a result of numerical round-off. My guess would be a conditioning issue, where one of the values is very small relative to the others and gets rounded off to 0 in the denominator, but is big enough to still impact the numerator. I'll leave the debug lines in and let you know if it happens again, but am not concerned - this stuff end up happening fairly regularly across the codebase.

I think the np.clip got missed in your latest commit, possible to add that quickly before merge?

I think the np.clip got missed in your latest commit, possible to add that quickly before merge?

Oh sorry, of course; I put the np.clip() back in.

Thanks! Need to move the close paren to fix the test failure but the MR is looking good.

I just spent a few minutes trying to reproduce but couldn't - it's a tricky edge case to trigger.

You had mentioned before that there were floating point round-off errors on the macos github CI runners (#28823 (comment)); I think that was the original motivation to put the np.clip() in in the first place. Any chance we can trigger the edge case there once again?

(Not that I'm against the np.clip(), but I'm still surprised by the occurrence of the quotient >1. I thought the IEEE 754 standard for floating-point arithmetic requires that multiplication and addition should be correctly rounded, so I thought we would be in the clear... unless macos would not conform to IEEE 754, which would be also quite remarkable... Or I don't quite understand all of the correctly, and I should go read some more Kahan)

If I happen to catch it again I'm happy to add it as a test! I think I lost the specific test case in the other MR that was erroring, and it was dependent on the specific math operations being performed so it likely wouldn't translate to this new code.

To be clear, I don't think it's a macos-specific problem, the issue is more that there isn't perfect determinism across different platforms / python versions / etc (even with perfect determinism within a configuration). Everything is being rounded "correctly", it's just that if we're right on the edge of floating point tolerance, complier/processor/implementation/optimization differences can have different results, and values infinitesimally on the wrong side of 1 will stack up unfavorably with further operations. For a simple concrete example, some systems in calculating (a + b)*c might distribute that multiplication to a*c + b*c, and those can result in different rounding of the results. Especially if for example a >> b and so b doesn't make it into the mantissa of a + b (addition and subtraction are not associative in floating point math, which is pretty unintuitive IMO).

[...] it's just that if we're right on the edge of floating point tolerance, complier/processor/implementation/optimization differences can have different results, and values infinitesimally on the wrong side of 1 will stack up unfavorably with further operations. For a simple concrete example, some systems in calculating (a + b)*c might distribute that multiplication to a*c + b*c, and those can result in different rounding of the results.

I'm just trying to wrap my head around what the optimization difference could have been that would lead to the observed result. '(a+b)c' -> 'ac + b*c' does not quite fit the bill... My hope is that an actual example would shed light on it...