Vector3 normalized() does not return normalized vector for very small nonzero vectors

[wareya]

Godot version

4.0.stable

System information

Windows 10

Issue description

See title. Normalized vectors are required for a lot of operations, like axis rotation, so normalized() should always return a real normalized vector.

Extremely small nonzero vectors can be produced by normal, innocuous game code when doing complicated vector math involving cross products etc.

Steps to reproduce

Create a new Node3D scene with the following script, then run the scene and wait for the errors to start coming in:

func _ready() -> void:
    seed(85123485)

func _process(delta: float) -> void:
    var f = (Vector3(randf()-0.5, randf()-0.5, randf()-0.5)*0.00000000000000001*0.00001).normalized()
    if !f.is_normalized():
        print(f)
        rotate(f, 1.0)

Values this small CAN AND WILL pop up in real game code without doing things this silly (see "Issue description"). This is PURELY to reproduce the issue in a very small amount of code.

Minimal reproduction project

N/A

[bitsawer]

I think normalized() does always return a normalized vector unless x, y and z are really, truly zero. This kind of makes sense, I mean what would be the normal vector of a zero vector? is_normalized() uses a small epsilon to work around float imprecision, which I think is also reasonable. Fortunately, it's pretty easy to implement your own is_normalized() is you need certain precision, but in this case I think you should probably be prepared for the fact that generating a random vector can sometimes give zero or extremely small vectors that fall below is_normalized() epsilon check.

[smix8]

The Vector3 normalized() always returns a normalized Vector3 but it is still a float operation.

The Vector3 is_normalized() check therefor uses is_equal_approx() with a UNIT_EPSILON which is 0.001 by default and with PRECISE_MATH_CHECKS builds 0.00001.

[wareya]

I think normalized() does always return a normalized vector unless x, y and z are really, truly zero.

The Vector3 normalized() always returns a normalized Vector3 but it is still a float operation.

No. The values printed by the given reproduction code are sometimes zero, but not usually. Some examples:

(0.38004, 0.631367, -1.172043)
(0.065905, 0.562576, 1.065786)
(-0.560717, 0.235717, 0.862024)
(1.16271, 0.111317, -0.170676)

These were produced by the given reproduction code, and their length is very different from both one and zero.

[bitsawer]

I think this is a floating point precision issue, normalizing a very small (or large) 3d vector involves sqrt() and more division, which can produce pretty inaccurate numbers with very large or small values. Might work better with 64-bit double precision build, some of these issues are just fundamental with 32-bit floats and have to be taken into consideration when doing calculations like this. There might be some algorithms to produce more solid results, but they will probably be slower too.

[smix8]

When you are not on a double-precision build of Godot multiplying with 0.00000000000000001 is too much for normal float, remove 3-4 zeros and it works again. The editor print displays zero but it is not stored as zero and the internal operations that involve further multiplications and sqrt freak out and flow over.

[wareya]

When you are not on a double-precision build of Godot multiplying with 0.00000000000000001 is too much for normal float, remove 3-4 zeros and it works again.

This is not the issue. It is just an easy way to replicate the issue in a small amount of code. Getting very small nonzero vectors like this is common when doing complicated vector math, as I said in the opening post. I ran into this issue trying to help someone troubleshoot custom physics code that was breaking because of the cross product operation between two normalized vectors producing an extremely small but nonzero vector.

[bitsawer]

Floating point inaccuracies are just something you need to be prepared for and handle when they occur, things will start to break down when you get too large or small values. In this case it might be possible to alter the normalized() to always produce values in range that is_normalized() will accept, but it would just be a bandaid hiding the real problem. Which in this case is expecting too much accuracy from a poor 32-bit float. Like mentioned, Godot can be built to use 64-bit double precision floats which will give some more breathing room, but even they will suffer from the same fundamental problem when pushed to the limits. Just don't expect too much from floating point numbers and you'll live a happier life.

[lawnjelly]

This is partly why a lot of normalize() implementations use an epsilon. The Godot one doesn't:

	real_t lengthsq = length_squared();
	if (lengthsq == 0) {
		x = y = z = 0;
	} else {
		real_t length = Math::sqrt(lengthsq);
		x /= length;
		y /= length;
		z /= length;
	}

If you might get near zero length input, you really need to do it manually anyway, because you could get zero length output.
i.e.

	real_t lengthsq = length_squared();
	if (lengthsq < EPSILON) {
	// DO SOMETHING ELSE
	} else {
		real_t length = Math::sqrt(lengthsq);
		x /= length;
		y /= length;
		z /= length;
	}

(edited my answer a little because it can be done)
This can be done in core. We could change to:

bool normalize()
{
	real_t lengthsq = length_squared();
	if (lengthsq < EPSILON) {
	return false;
	}
	real_t length = Math::sqrt(lengthsq);
	x /= length;
	y /= length;
	z /= length;
	return true;
}

The only issue is that GDScript makes warnings when you don't use a return value I think?

Using from gdscript would then become:

var myvec = Vector3(0, 0, 0)
if myvec.normalize() == false:
	# do some error handling

# normal routine

This could be added if there was popular support, might have to have a separate function name for backward compatibility / preventing warnings. Something like try_normalize()? 🤔

Or alternatively revamp my old PR #56567 :

real_t Vector2::normalize() {
	real_t l = x * x + y * y;
	if (l != 0) {
		l = Math::sqrt(l);
		x /= l;
		y /= l;
		return l;
	}
	return 0.0;

But test against an EPSILON instead of zero. This would mean as well as returning the length (which is useful), you could test against a returned length of zero for failure.

[lawnjelly]

Thinking about this, we should probably add to the documentation that normalized() will return erroneous results if the Vector is near zero length.

[Poppyanimal]

I should note, in situations where a non-zero but really small vector is normalized, the result does lead to something closer to being normalized a lot of the time. This means calling .normalized().normalized() on a vector you want normalized can result in this error occuring less than calling it once. Basically, this can be fixed in the normalized() method itself without having to complain about how its the fault of floating point error or that we should just deal with it or that there is no solution. It should be expected behaviour that normalized() returns a vector that passes a isNormalized() check when fed a non-zero vector. Hell, it doesnt even need to have an "exact" length of 1, the point is it just needs to be normalized enough to pass normalization checks like the one that occurs inside of the rotate(axis, angle) method, and the fact that normalize() can return a non-zero vector that isnt normalized should be viewed as a bug that should be fixed and not some feature we throw tape around and ignore and leave to fester

[clayjohn]

@lawnjelly Wouldn't it be enough to add an epsilon to the length check?

something like:

real_t lengthsq = length_squared();
	if (lengthsq < CMP_EPSILON2) {
		x = y = z = 0;
	} else {
		real_t length = Math::sqrt(lengthsq);
		x /= length;
		y /= length;
		z /= length;
	}

[lawnjelly]

@lawnjelly Wouldn't it be enough to add an epsilon to the length check?

That is probably a good idea for the existing function..

however the problem remains that effectively when using this, unless you can guarantee non small lengths, you need to check for the error condition. This means you have to call e.g. length_squared() on the result, and check for 0, which is both inefficient, and messy in terms of client code side. Whereas with the alternate functions, the success or failure is returned as part of the function.

Imo the whole normalized() function isn't a good way of doing this, but we are stuck with it for backward compatibility, and we should probably add a more useful function, either returning success directly, or returning the length as in #56567 .

[kleonc]

something like:

real_t lengthsq = length_squared();
if (lengthsq < CMP_EPSILON2) {
	x = y = z = 0;
} else {
	real_t length = Math::sqrt(lengthsq);
	x /= length;
	y /= length;
	z /= length;
}

I wonder if alternatively something like this would make sense:

real_t lengthsq = length_squared();
if (lengthsq < CMP_EPSILON2) {
	if (x == 0 && y == 0 && z == 0) {
		// x = y = z = 0; // If we'd want to avoid +/- zeros.
		return; // Exact zero.
	}
	// Near-zero which could fail normalizing. Let's scale it by some power of two.
	// Is there some specific value which could guarantee the normalization down below won't fail?
	const real_t multiplier = 1024.0 * 1024.0;
	x *= multiplier;
	y *= multiplier;
	z *= multiplier;
	lengthsq = length_squared();
}
real_t length = Math::sqrt(lengthsq);
x /= length;
y /= length;
z /= length;

I mean if the input is such near-zero value then it's indeed likely to be not too meaningful because of the potentially already cumulated floating points errors, and scaling it just to not fail might enlarge the already in there error. So I'm not convinced doing something like this makes sense at all, just sharing a thought. 🙃

[lawnjelly]

So there are essentially two problems caused by small values:

1. normalization resulting in non-unit length output
2. lower resolution of the direction, caused by float error at small values

Normalization twice or pre-multiplying sound like viable ideas for dealing with (1), but don't deal with (2). And with (2) some use cases might be more sensitive to breakdown in direction resolution than others. 🤔

That suggests something like:

// breakdown in unit length
if (lengthsq < CMP_EPSILON2) {
	// breakdown in direction
	if (lengthsq < CMP_EPSILON3) {
		// zero or return zero etc
	}
	// pre-multiply
}
// normal routine

This makes things a bit slower, but performance sensitive loops using normalization should probably be hand optimized (e.g. SIMD / reciprocal sqrt etc).

The question is whether we treat the threshold for breakdown in direction as different from the breakdown in unit length, and whether they are significantly different enough to warrant this. I'll have a look in the testbed I wrote yesterday. 👍

UPDATE:
I'm admittedly not an expert with super small floating point representation, but so far the accuracy of the direction seems pretty good below the epsilon needed for unit length. This may be because the mantissa has a lot of resolution even with low numbers.

One thing we may have to be aware of is subnormals:
https://en.wikipedia.org/wiki/Subnormal_number

And making sure the epsilons work ok even when subnormals are flushed to zero, and I think making sure it works well in --fast-math may be a good idea, because we have no guarantees about the mode the CPU will be in, what flags have been used to compile Godot + libraries, what CPU it is running on etc, all of which could affect things afaik.

For instance with --fast-math, the square length calculation can result in zero, and yet the x and y still contain enough information to calculate the direction! Things are quite strange at this level and I'm always paranoid compiler optimization is not doing what we ask it to, even at -O0, so looking at the assembly is the only way to verify.

Most implementations of normalize() I've seen in the past don't use two epsilons, so we may be on thin ground here, but it may be warranted because we have the peculiar case of math checks for unit length in other functions.

What I might do is to write it as SIMD instrinsics and make sure the epsilons work there, because that's the most approximate version that we would encounter I guess.

[lawnjelly]

Ok incredibly inefficient but this is just for ensuring 32 bit calculations, with flushing denormals to zero, I think this is right:

	void normalize_SIMD()
	{
		_MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON);
		_MM_SET_DENORMALS_ZERO_MODE(_MM_DENORMALS_ZERO_ON);
		__m128 x0 = _mm_set1_ps(x);
		__m128 x1 = _mm_set1_ps(x);
		__m128 resx = _mm_mul_ps(x0, x1);
		
		__m128 y0 = _mm_set1_ps(y);
		__m128 y1 = _mm_set1_ps(y);
		__m128 resy = _mm_mul_ps(y0, y1);
		
		__m128 res = _mm_add_ps(resx, resy);
		__m128 epsilon = _mm_set1_ps(EPSILON);

		__m128 cmp_res = _mm_cmple_ps(res, epsilon);
		
		float sl;
		_mm_store_ss(&sl, res);
		uint32_t cmp;
		_mm_store_ss((float *) &cmp, cmp_res);
		
		if (!cmp)
		{
			__m128 sq = _mm_sqrt_ps(res);
			x0 = _mm_div_ps(x0, sq);
			y0 = _mm_div_ps(y0, sq);

			_mm_store_ss(&x, x0);			
			_mm_store_ss(&y, y0);			
		}
		else
		{
			out("\t.. NORMALIZATION_FAILED .. ");
		}
	}

I'm not convinced the FLUSH_ZERO and DENORMAL_ZERO are working though, as through my debugger the smallest values stored go down to 1.4013e-45f which corresponds to the lowest subnormal value, rather than 1.175e-38 which should be the lowest normal value.
Interestingly even at lower epsilon 1.4013e-45f there is still enough resolution to decently encode the direction. That is the error in a radian angle in 2D is of the order of 5.9604644775390625e-08, with input vector length of circa 6.0742653774578532e-19, but I'm not sure this can be relied on.

But this does suggest we could get away with 2 epsilons as described above, one threshold for a multiplier to ensure unit length, and a second for complete failure, but probably above zero to ensure consistent results on different platforms.

I would suggest maybe staying with CMP_EPSILON2 for the threshold for using a multiplier, then in theory we could go down to some factor of FLT_MIN (which by definition is meant to be the smallest value before using denormals).

As an aside, I suspect that in the recalculation of square length after applying a multiplier in @kleonc 's version, can we simplify it to:

lengthsq *= (multiplier * multiplier);

However, given we are dealing with precision issues it may not be wise. 🤔

[aaronfranke]

Since nobody has mentioned it, a high-level work-around would be to check if the vector is nearly zero before getting the normalized vector.

func _ready() -> void:
    seed(85123485)

func _process(delta: float) -> void:
    var test_vector := Vector3(randf()-0.5, randf()-0.5, randf()-0.5)*0.00000000000000001*0.00001
    var f = checked_normalized(test_vector)
    print(f.is_normalized())

func checked_normalized(vec: Vector3) -> Vector3:
    if vec.is_zero_approx():
        return Vector3.ZERO
    return vec.normalized()

Another option is we could add this into normalized(), instead of checking if the length is equal to zero we could check if it's less than a small epsilon (in which case we also don't need to return the length to solve this problem).

Also, I want to mention that you generally should not expect values like 0.00000000000000001 to work at all. If it does, great, but don't expect it to be reliable. If you often need super precise values like this, I recommend compiling Godot with double precision support enabled.

[Poppyanimal]

Not a suggestion, but I wanted to apologize for getting hot headed with my previous comment earlier this thread. I failed to maintain a non-adversarial stance and my own emotions overtook me. Thank you for considering my suggestions/comments in spite of myself at the time. I really like lawnjelly's current idea, the one in the issue just mentioned, and it would be interesting to see that after it gets some testing. This problem is one that had effected me relatively recently and I'm happy to see it be addressed with such speed and consideration as it has been. Thank you all for spending time looking and thinking about this problem as you have been. Again, sorry for my behaviour and thank you for your time

[reduz]

I am generally more on the side that this should not be fixed, and that whathever you are doing if you have an issue with this is the actual problem. An epsilon check in the normalize will probably not fix your problem anyway so, at most, could be a debug build feature that errors as if it was zero.