diff --git a/source/blender/blenlib/BLI_math_vector.h b/source/blender/blenlib/BLI_math_vector.h
index e6e46a4d8cd..9ef2dc781de 100644
--- a/source/blender/blenlib/BLI_math_vector.h
+++ b/source/blender/blenlib/BLI_math_vector.h
@@ -285,6 +285,7 @@ void axis_sort_v3(const float axis_values[3], int r_axis_order[3]);
 /***************************** Array Functions *******************************/
 /* attempted to follow fixed length vertex functions. names could be improved*/
 double dot_vn_vn(const float *array_src_a, const float *array_src_b, const int size) ATTR_WARN_UNUSED_RESULT;
+double len_squared_vn(const float *array, const int size) ATTR_WARN_UNUSED_RESULT;
 float normalize_vn_vn(float *array_tar, const float *array_src, const int size);
 float normalize_vn(float *array_tar, const int size);
 void range_vn_i(int *array_tar, const int size, const int start);
diff --git a/source/blender/blenlib/intern/math_vector.c b/source/blender/blenlib/intern/math_vector.c
index 19865aa00bd..3fbddacaea2 100644
--- a/source/blender/blenlib/intern/math_vector.c
+++ b/source/blender/blenlib/intern/math_vector.c
@@ -665,6 +665,11 @@ void axis_sort_v3(const float axis_values[3], int r_axis_order[3])
 
 /***************************** Array Functions *******************************/
 
+MINLINE double sqr_db(double f)
+{
+	return f * f;
+}
+
 double dot_vn_vn(const float *array_src_a, const float *array_src_b, const int size)
 {
 	double d = 0.0f;
@@ -677,9 +682,20 @@ double dot_vn_vn(const float *array_src_a, const float *array_src_b, const int s
 	return d;
 }
 
+double len_squared_vn(const float *array, const int size)
+{
+	double d = 0.0f;
+	const float *array_pt = array + (size - 1);
+	int i = size;
+	while (i--) {
+		d += sqr_db((double)(*(array_pt--)));
+	}
+	return d;
+}
+
 float normalize_vn_vn(float *array_tar, const float *array_src, const int size)
 {
-	double d = dot_vn_vn(array_tar, array_src, size);
+	double d = len_squared_vn(array_src, size);
 	float d_sqrt;
 	if (d > 1.0e-35) {
 		d_sqrt = (float)sqrt(d);
diff --git a/source/blender/python/mathutils/mathutils_Vector.c b/source/blender/python/mathutils/mathutils_Vector.c
index 3dc953e22e9..96e68317e2d 100644
--- a/source/blender/python/mathutils/mathutils_Vector.c
+++ b/source/blender/python/mathutils/mathutils_Vector.c
@@ -1888,17 +1888,6 @@ static PyObject *Vector_neg(VectorObject *self)
 	return Vector_CreatePyObject_alloc(tvec, self->size, Py_TYPE(self));
 }
 
-/*------------------------vec_magnitude_nosqrt (internal) - for comparing only */
-static double vec_magnitude_nosqrt(const float *data, int size)
-{
-	/* return (double)sqrt(dot);*/
-	/* warning, line above removed because we are not using the length,
-	 * rather the comparing the sizes and for this we do not need the sqrt
-	 * for the actual length, the dot must be sqrt'd */
-	return dot_vn_vn(data, data, size);
-}
-
-
 /*------------------------tp_richcmpr
  * returns -1 exception, 0 false, 1 true */
 static PyObject *Vector_richcmpr(PyObject *objectA, PyObject *objectB, int comparison_type)
@@ -1933,15 +1922,15 @@ static PyObject *Vector_richcmpr(PyObject *objectA, PyObject *objectB, int compa
 
 	switch (comparison_type) {
 		case Py_LT:
-			lenA = vec_magnitude_nosqrt(vecA->vec, vecA->size);
-			lenB = vec_magnitude_nosqrt(vecB->vec, vecB->size);
+			lenA = len_squared_vn(vecA->vec, vecA->size);
+			lenB = len_squared_vn(vecB->vec, vecB->size);
 			if (lenA < lenB) {
 				result = 1;
 			}
 			break;
 		case Py_LE:
-			lenA = vec_magnitude_nosqrt(vecA->vec, vecA->size);
-			lenB = vec_magnitude_nosqrt(vecB->vec, vecB->size);
+			lenA = len_squared_vn(vecA->vec, vecA->size);
+			lenB = len_squared_vn(vecB->vec, vecB->size);
 			if (lenA < lenB) {
 				result = 1;
 			}
@@ -1956,15 +1945,15 @@ static PyObject *Vector_richcmpr(PyObject *objectA, PyObject *objectB, int compa
 			result = !EXPP_VectorsAreEqual(vecA->vec, vecB->vec, vecA->size, 1);
 			break;
 		case Py_GT:
-			lenA = vec_magnitude_nosqrt(vecA->vec, vecA->size);
-			lenB = vec_magnitude_nosqrt(vecB->vec, vecB->size);
+			lenA = len_squared_vn(vecA->vec, vecA->size);
+			lenB = len_squared_vn(vecB->vec, vecB->size);
 			if (lenA > lenB) {
 				result = 1;
 			}
 			break;
 		case Py_GE:
-			lenA = vec_magnitude_nosqrt(vecA->vec, vecA->size);
-			lenB = vec_magnitude_nosqrt(vecB->vec, vecB->size);
+			lenA = len_squared_vn(vecA->vec, vecA->size);
+			lenB = len_squared_vn(vecB->vec, vecB->size);
 			if (lenA > lenB) {
 				result = 1;
 			}