[raymath] Make every normalize function similar

[planetis-m]

Follow up to #2982 Makes normalize functions consistent, moved sqrtf call after the if statement which causes its assembly to be slightly better in both clang 17 and gcc 13. It now uses one rsqrtss instruction instead of sqrtss + divss combo https://godbolt.org/z/5K475o88n plus in clang it's branchless.

[raysan5]

@planetis-m This change modifies multiple functions and it's quite sensible, did you implement some test example to verify everything works exactly the same way? Some test cases would be really useful in this situation.

[planetis-m]

I don't have tests cases for all the functions I changed, they're quite a lot, however my justification is that the number with square root 0 can only be 0. So delaying the call to the function sqrtf after the "if value is 0" statement doesn't make any difference in the result. In fact it was already done like that in the ClampValue functions.

[planetis-m]

Relevant SO question: https://stackoverflow.com/questions/73988574/given-x0-is-it-possible-for-sqrtx-0-to-give-a-floating-point-error Maybe it's preferable to remove the check altogether?

[planetis-m]

I have tried the example in the SO answer and the condition xxyy == 0 was hit with both clang and gcc. I've also asked LLMs to try to break it with their own values for x and y but couldn't. I think it's ok to merge.

#include <stdio.h>
#include <math.h>

int main() {
    double x = 0x1p-1021;
    double y = 0x2p-1020; // or y = 0
    double xxyy = x * x + y * y;
    if (xxyy == 0) printf("caught\n");
    printf("%a %a\n", x, y);
    printf("%a %a\n", xxyy, sqrt(xxyy));

    return 0;
}

But I can switch the checks to abs(xxyy) <= EPSILON, if it's any better.

[raysan5]

@planetis-m This PR changes many functions and it should be carefully tested to verify every new implementation behaves exactly like previous one. Also note that comparing float values with == could be really dangerous due to rounding issues. Afaik, it's a not recomended practize.

[JayLCypher]

I notice there are some mismatches checking the length > 0.0f some places and != 0.0f in other places.
If I'm not mistaken, in real math a square root cannot be negative. Ergo, length > 0.0f should probably be preferred.

Regardless, most of these functions follow the same pattern of Calculate length using sqrtf and divide each vector member by length. In reality you may have V.X * Y / LEN, but since Y is 1.0f, it becomes V.X / Len. Each function could therefore be reduced to:

// Normalize provided vector
RMAPI Vector2 Vector2Normalize(const Vector2 v)
{
    const float length = sqrtf((v.x*v.x) + (v.y*v.y));
	return (length > 0.0f) ? CLITERAL(Vector2){v.x / length, v.y / length} : CLITERAL(Vector2){0};
}

Step 1: Calculate Length
Step 2: Length not 0 ? Normalize length or return 0 vector.

I may be wrong, so further observation is required.

Edit:
I notice that raymath.h is standalone and doesn't include the CLITERAL macro, so the safer bet is probably to bind result to a variable per function. Otherwise we'd need an ifndef clause on CLITERAL and redefine it in raymath.h based on if we're c++ or not.

Currently no cliteral and compound literal would look like:

// Normalize provided vector
RMAPI Vector2 Vector2Normalize(const Vector2 v)
{
    const float length = sqrtf((v.x*v.x) + (v.y*v.y));
	Vector2 result = {0};
	if (length > 0.0f) {
		result.x = v.x / length;
		result.y = v.y / length;
	}
	return result;
}

I'm not sure what your thoughts are on const-correctness and conciseness here @raysan5 , since Raylib is also for educational use, a smaller, let's call it "more clever" way of doing things is not necessarily better for teaching and learning.

[raysan5]

TODO: All affected functions should be properly tested and compared to current ones, in this case some unit tests would be required.

I will keep this issue open for some time in case anyone wants to do that work. If not, I will just close the issue because current implementation is proved to be working correctly and the benefits of this PR are dubious.

[planetis-m]

Sorry I do not know how to make comprehensive tests. Because I do not know what exactly to test for! Are we doing unit testing or integration testing? Are there any other tests files I can draw inspiration from?

I did another test though with the example raymath_vector_angle and used assert(fpclassify(ilength) == FP_NORMAL) I noticed when the length of the vector was zero, it correctly hit the if (lengthSq == 0.0f) branch (verified with print debugging) and this assert was not triggered. Then I removed the condition all together and it still works!

I would just remove the condition if (length == 0.0f) as other libraries, like glm, don't have it either. And it also affects performance.

[planetis-m]

Sokol: https://github.com/floooh/sokol/blob/c54523c078e481d3084fa0b4630d2ce3d3e1e74f/util/sokol_gl.h#L3158
Glm: https://github.com/g-truc/glm/blob/33b4a621a697a305bc3a7610d290677b96beb181/glm/detail/func_geometric.inl#L104
vemath: https://github.com/PistonDevelopers/vecmath/blob/763fedccfd4a72048e7af6870e4f6368949fee52/src/lib.rs#L951
vmath: https://github.com/treeform/vmath/blob/7282ae1247f2f384ebeaec3826d7fa38fd0e1df1/src/vmath.nim#L1518

Not to say that I could still be wrong, but please to advance the conversation, try to show me actual breakage, it's only fair to ask as I did so much research in the topic.

[orcmid]

@planetis-m

Not to say that I could still be wrong, but please to advance the conversation, try to show me actual breakage, it's only fair to ask as I did so much research in the topic.

PPS: In reviewing more of this, I see that you did explain that you were avoiding sqrtf when the result would be 0 regardless. So I apologize about my breaking-change objection. What is needed is agreement on what result should be provided in those cases.

PS: There are differences in edge cases here and, with a quick look, it seems as if that is what the proposed changes are about. The screening of float values for exactly 0.0f is also not a great idea although I can see what is intended and how silent results have been produced (inconsistently?) when such an edge is determined. In that case, these are not breaking changes in the usual sense. This needs to be made more clear and the differences in edge behavior reviewed.

I think this could have been explained better.

I withdraw my previous comments.

[planetis-m]

Ok I took the time to write a fuzz target that you can compile with: clang++ -fsanitize=fuzzer,address -o target target.cpp

#include <math.h>
#include <stddef.h>
#include <stdint.h>
#include <cstring>
#include <iostream>

// Definition of the Vector2 struct
typedef struct Vector2 {
  float x;
  float y;
} Vector2;

// Implementation of the Vector2Normalize function
Vector2 Vector2Normalize(Vector2 v) {
  Vector2 result = { 0 };
  float length = sqrtf((v.x*v.x) + (v.y*v.y));
  if (length > 0) {
    float ilength = 1.0f/length;
    result.x = v.x*ilength;
    result.y = v.y*ilength;
  }
  return result;
}

Vector2 Vector2Normalize2(Vector2 v) {
  Vector2 result = { 0 };
  float lengthSq = (v.x*v.x) + (v.y*v.y);
  if (lengthSq == 0.0f) lengthSq = 1.0f;
  float ilength = 1.0f/sqrtf(lengthSq);
  result.x = v.x*ilength;
  result.y = v.y*ilength;
  return result;
}

// Fuzz target function
extern "C" int LLVMFuzzerTestOneInput(const uint8_t *data, size_t size) {
  if (size < sizeof(Vector2)) {
    // Not enough data to form a valid Vector2
    return 0;
  }
  // Create a Vector2 object from the input data
  Vector2 v;
  memcpy(&v, data, sizeof(Vector2));
  // // Call the function under test
  // Vector2 result = Vector2Normalize(v);
  // // Optionally, you can add some checks here to validate the result
  // // For example, you could check that the result vector has a length of 1 (normalized)
  // float result_length = sqrtf((result.x * result.x) + (result.y * result.y));
  // if (result_length < 0.99f || result_length > 1.01f) {
  //   // Report an error if the result is not properly normalized (allowing some tolerance)
  //   std::cerr << "v: (" << v.x << ", " << v.y << ")\n";
  //   std::cerr << "result: (" << result.x << ", " << result.y << ")\n";
  //   abort();
  // }
  // Call both normalization functions
  Vector2 result1 = Vector2Normalize(v);
  Vector2 result2 = Vector2Normalize2(v);

  // Check that the results are equal
  // const float tolerance = 0.001f;
  // if (fabs(result1.x - result2.x) > tolerance || fabs(result1.y - result2.y) > tolerance) {
  if ((result1.x != result2.x) || (result1.y != result2.y)) {
    // If the results differ significantly, print the vectors and abort
    std::cerr << "Mismatch:\n";
    std::cerr << "v: (" << v.x << ", " << v.y << ")\n";
    std::cerr << "result1: (" << result1.x << ", " << result1.y << ")\n";
    std::cerr << "result2: (" << result2.x << ", " << result2.y << ")\n";
    abort();
  }
  return 0; // Indicate successful execution
}

There's a mismatch with very small numbers:

Mismatch:
v: (7.91498e-24, 2.17773e-32)
result1: (0, 0)
result2: (7.91498e-24, 2.17773e-32)

The first function returns the zero vector while the modified one, the original vector itself.
lengthSq is set to zero with these numbers so the branch if (lengthSq == 0.0f) lengthSq = 1.0f; is executed.

Both are not unit vectors obviously, so the function's contract is broken regardless.

I also caught other cases that produce different results based on the compilation flags:

v: (1.12104e-44, 0)

Produces either itself or (0, 0) depending if -ffast-math -march=native -O3 is used.

Every vector with a nan component is also causing trouble:

v: (1, 0.0/0.0)
result1: (0, 0)
result2: (nan, nan)

If your calling code checks the results with a tolerance there's no change in output. But if it doesn't I guess there's some difference. But does it really break your code, please let me know.

[orcmid]

@planetis-m @raysan5

There's a mismatch with very small numbers:

Mismatch:
v: (7.91498e-24, 2.17773e-32)
result1: (0, 0)
result2: (7.91498e-24, 2.17773e-32)

The first function returns the zero vector while the modified one, the original vector itself. lengthSq is set to zero with these numbers so the branch if (lengthSq == 0.0f) lengthSq = 1.0f; is executed.

Both are not unit vectors obviously, so the function's contract is broken regardless.

Yes, the problem is what should always be returned is a unit vector. There are two problem cases: Underflow of the squares and overflow of the squares summed.

I don't think raymath.h should deal with the NaN or +Inf cases. Those are cases for the user to handle. Also whether or not denormalized small values are allowed.

There is mathematically only one edge case for which there is no answer. That's for (0,0). There's no way to assign a direction of a true 0-length vector.

1. So the question for the exact (0,0) cases is, what should the result be?
2. It should be possible to work around the underflow cases, but the question is whether or not it is worth it. And if not, see question 1.

I suspect these cases do not come up in correct usages of raylib.h. But having consistent responses would be valuable. I support that quest.

I do think it is a defect not to address this. And the edge-case result should be made known and part of the interface contract.

PS: It seems that returning (0,...,0) is the most popular result for the case where there is no normalization. https://stackoverflow.com/questions/722073/how-do-you-normalize-a-zero-vector

[planetis-m]

To add to what @orcmid said, I think the proper check for the function is this:

let unit_vector = normalize(v);
let tmp = length(unit_vector);
if (tmp < 0.99f || tmp > 1.01f) {
  // Handle fail case
}

which works in both cases.

Or:

let unit_vector = normalize(v);
if (unit_vector == {0}) {
  // Handle fail case
}

Which only works with the original approach. Let's please close this PR either way. It's been dragged for too long.

[JayLCypher]

If we think about this some more, we can recognize that the only situation of which a division by zero can occur in these functions are on Zero-vectors only.

By virtue of squaring, all negative values will be positive, and both positives are summed, ergo no value can cancel each other out ending up in a square root of zero. The only place this can happen is if both parts of the vector are already zero.
Thus, the only error case for Normalize is overflow in sqrtf, leading to +inf returned and division by infinity.

This means that we can by contract assert that either part of the argument vector has to be non-zero, and calling the function with a Zero-vector is illegal.
All other values will most likely become transient, meaning calling the functions with part NaN or part Inf will return NaN or Inf respectively. (This could also be removed by contract.)

This will result in a much smaller function, but will require handling of the contract (by using assert) and users will have to opt-out of the assert check for well-formed programs after verification (using -DNDEBUG flag)

Resulting code (Godbolt link):

Vector2 Vector2Normalize(const Vector2 v) {
  assert(v.x > 0.0f || v.y > 0.0f);
  const float length = sqrtf((v.x*v.x) + (v.y*v.y));
  Vector2 result = {
    v.x / length,
    v.y / length,
  };
  return result;
}

@planetis-m I don't know if you can fuzz this with Zero-vector being illegal, but that would be very cool!

[orcmid]

@JayLCypher @planetis-m @raysan5

Let's review the bidding here. I believe there are three factors to consider.

1. The prospect of floating-point exceptions should simply be avoided, although overflows might be possible. In any case, the treatment of such exceptions, and any mitigation of them needs to be left to the user of raymath. That's consistent with all of the raymath functions (and raylib ones too). There should simply be no divisions by 0 (or maybe too-small floats) in the functions.
2. Consider the geometric meaning of normal vectors. It happens that (0,0) has no normal. That is, there is no direction to a point away from (0,0) that can be normalized to a unit vector. One rather natural way to treat coordinates for which there is no direction from (0,0) to draw a normal toward/through is to return (0,0) as the no-normal-here result.
3. With regard to numerical analysis, it is possible, as @planetis-m managed to detect, for v.x*v.x + v.y*vy to be zero when neither of v.x and v.y are zero. It's also the case that should be considered close enough to 0 to be treated as 0.0f. Since raymath.h includes <math.h>, FLT_EPSILON might be an appropriate tolerance. Or use the EPSILON defined directly in raymath.h.

Applying these ideas leads to a flavor of the @planitis-m recommendation in raylib style:

RMAPI Vector2 Vector2Normalize(Vector2 v)
{
    Vector2 result = { 0 };
    float lengthSQ = (v.x*v.x) + (v.y*v.y);

    if (lengthSQ  >  (FLT_EPSILON)*(FLT_EPSILON))
    {
        float ilength = 1.0f/sqrtf(lengthSQ);
        result.x = v.x*ilength;
        result.y = v.y*ilength;
     }
     
    return result;
}

I favor the idea of returning exactly {0.0f, 0.0f} when the operand's length is considered too close to 0.

[JayLCypher]

Right, I made an error in reasoning for the values of the range (-1, 1), where the product of squaring actually becomes smaller. One could almost be forgiven to introduce a separate handling of this special edge case.

Regardless, an easy fix is to clamp the value to FLT_EPSILON / Raymath EPSILON / some minimal value, accepting certain inaccuracies. Doing so with the squared value using fmaxf or a ternary.

More observations:

• The 1.0f / sqrt(length) ilength variable introduces an inaccuracy of around +-0.00000005960464477539 or less.
• If either value is 0, the square root of the squared value is itself, ergo no sqrtf call is necessary and the normalized value is 1 or -1.
• We have an edge case of -infinity, causing -infinity / infinity (or -infinity * 1 / infinity) due to squaring, resulting in -nan. (Bug?)
• Introduction of tolerances will introduce inaccuracies that the original function did not have (returning 0 where original function modifies value, potentially breaking?). Additionally, it does not return -1 or 1 for the axis respectively on single-digit vectors of small numbers.

Godbolt Fuzzer (changed to C)

Man, now I see why floating-point is such a doof!

[planetis-m]

Your feedback has been invaluable! This discussion has not only enhanced my understanding, but helped grow as a programmer. Based on the conversation, it seems that either Vector2Normalize2 or Vector2Normalize4 from @JayLCypher's reply might be the functions we're interested in? Correct?

Also what about every function that inlines normalize, like Vector3Rotate, MatrixLookAt, etc which use:

    length = sqrtf(v.x*v.x + v.y*v.y + v.z*v.z);
    if (length == 0.0f) length = 1.0f;
    ilength = 1.0f/length;
    vz.x *= ilength;
    vz.y *= ilength;
    vz.z *= ilength;

Should that be changed to:

    length = sqrtf(v.x*v.x + v.y*v.y + v.z*v.z);
    if (length == 0.0f) {
        vz.x = 0;
        vz.y = 0;
        vz.z = 0;
    } else {
        vz.x /= length;
        vz.y /= length;
        vz.z /= length;
    }

PS. To ignore certain inputs when fuzzing all you need to do is check for them and return early.

[raysan5]

I'm closing this PR. If it requires another review in the future for functions alignment, feel free to open a new PR.