Start simpler.

Check expected behavior of...

Nan, +0, -0, +1, -1, +Inf, -Inf, +Eps, -Eps, +1+Eps, +1-Eps, -1+Eps, -1-Eps

...on each coordinate when supplying a sensible constant for the other coordinates.

Check the outer product of the above list taken across all coordinates.

Check that you get the same answer when you call the function twice in each of these outer product cases.

This approach works for any floating point function. You will be surprised how many mistakes these tests will catch, inclusive of later maintenance on the code under test.