Revert "BLI: Refactor vector types & functions to use templates"

Reverted because the commit removes a lot of commits.

This reverts commit a2c1c368af48644fa8995ecbe7138cc0d7900c30.

45 files changed, 1877 insertions, 1457 deletions

diff --git a/source/blender/blenlib/BLI_array_store.h b/source/blender/blenlib/BLI_array_store.h
index 0be361d4ab9..68928f53e55 100644
--- a/source/blender/blenlib/BLI_array_store.h
+++ b/source/blender/blenlib/BLI_array_store.h
@@ -84,7 +84,7 @@ size_t BLI_array_store_calc_size_compacted_get(const BArrayStore *bs);
  */
 BArrayState *BLI_array_store_state_add(BArrayStore *bs,
                                        const void *data,
-                                       const size_t data_len,
+                                       size_t data_len,
                                        const BArrayState *state_reference);
 /**
  * Remove a state and free any unused #BChunk data.
diff --git a/source/blender/blenlib/BLI_array_utils.h b/source/blender/blenlib/BLI_array_utils.h
index 50fc2ce909b..202ae9a9786 100644
--- a/source/blender/blenlib/BLI_array_utils.h
+++ b/source/blender/blenlib/BLI_array_utils.h
@@ -52,7 +52,7 @@ void _bli_array_wrap(void *arr, uint arr_len, size_t arr_stride, int dir);
  * Access via #BLI_array_wrap
  */
 void _bli_array_permute(
-    void *arr, uint arr_len, const size_t arr_stride, const uint *order, void *arr_temp);
+    void *arr, uint arr_len, size_t arr_stride, const uint *order, void *arr_temp);
 #define BLI_array_permute(arr, arr_len, order) \
   _bli_array_permute(arr, arr_len, sizeof(*(arr)), order, NULL)
 #define BLI_array_permute_ex(arr, arr_len, order, arr_temp) \
@@ -152,7 +152,7 @@ bool _bli_array_is_zeroed(const void *arr, uint arr_len, size_t arr_stride);
  */
 bool _bli_array_iter_spiral_square(const void *arr_v,
                                    const int arr_shape[2],
-                                   const size_t elem_size,
+                                   size_t elem_size,
                                    const int center[2],
                                    bool (*test_fn)(const void *arr_item, void *user_data),
                                    void *user_data);
diff --git a/source/blender/blenlib/BLI_buffer.h b/source/blender/blenlib/BLI_buffer.h
index 7f577443cf7..e79d44fd934 100644
--- a/source/blender/blenlib/BLI_buffer.h
+++ b/source/blender/blenlib/BLI_buffer.h
@@ -74,7 +74,7 @@ enum {
 /**
  * \note Never decreases the amount of memory allocated.
  */
-void BLI_buffer_resize(BLI_Buffer *buffer, const size_t new_count);
+void BLI_buffer_resize(BLI_Buffer *buffer, size_t new_count);
 
 /**
  * Ensure size, throwing away old data, respecting #BLI_BUFFER_USE_CALLOC.
@@ -83,7 +83,7 @@ void BLI_buffer_resize(BLI_Buffer *buffer, const size_t new_count);
  * - Ignored (malloc'd).
  * - Cleared (when #BLI_BUFFER_USE_CALLOC is set).
  */
-void BLI_buffer_reinit(BLI_Buffer *buffer, const size_t new_count);
+void BLI_buffer_reinit(BLI_Buffer *buffer, size_t new_count);
 
 /**
  * Append an array of elements.
diff --git a/source/blender/blenlib/BLI_delaunay_2d.h b/source/blender/blenlib/BLI_delaunay_2d.h
index 1ee0be64cee..db0df95499f 100644
--- a/source/blender/blenlib/BLI_delaunay_2d.h
+++ b/source/blender/blenlib/BLI_delaunay_2d.h
@@ -215,9 +215,9 @@ void BLI_delaunay_2d_cdt_free(CDT_result *result);
 /* C++ Interface. */
 
 #  include "BLI_array.hh"
+#  include "BLI_double2.hh"
 #  include "BLI_math_mpq.hh"
-#  include "BLI_math_vec_mpq_types.hh"
-#  include "BLI_math_vec_types.hh"
+#  include "BLI_mpq2.hh"
 #  include "BLI_vector.hh"
 
 namespace blender::meshintersect {
diff --git a/source/blender/blenlib/BLI_double2.hh b/source/blender/blenlib/BLI_double2.hh
new file mode 100644
index 00000000000..0abff01ab2f
--- /dev/null
+++ b/source/blender/blenlib/BLI_double2.hh
@@ -0,0 +1,143 @@
+/*
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+#pragma once
+
+/** \file
+ * \ingroup bli
+ */
+
+#include "BLI_double3.hh"
+
+namespace blender {
+
+struct double2 {
+  double x, y;
+
+  double2() = default;
+
+  double2(const double *ptr) : x{ptr[0]}, y{ptr[1]}
+  {
+  }
+
+  double2(double x, double y) : x(x), y(y)
+  {
+  }
+
+  double2(const double3 &other) : x(other.x), y(other.y)
+  {
+  }
+
+  operator double *()
+  {
+    return &x;
+  }
+
+  operator const double *() const
+  {
+    return &x;
+  }
+
+  double length() const
+  {
+    return len_v2_db(*this);
+  }
+
+  friend double2 operator+(const double2 &a, const double2 &b)
+  {
+    return {a.x + b.x, a.y + b.y};
+  }
+
+  friend double2 operator-(const double2 &a, const double2 &b)
+  {
+    return {a.x - b.x, a.y - b.y};
+  }
+
+  friend double2 operator*(const double2 &a, double b)
+  {
+    return {a.x * b, a.y * b};
+  }
+
+  friend double2 operator/(const double2 &a, double b)
+  {
+    BLI_assert(b != 0.0);
+    return {a.x / b, a.y / b};
+  }
+
+  friend double2 operator*(double a, const double2 &b)
+  {
+    return b * a;
+  }
+
+  friend bool operator==(const double2 &a, const double2 &b)
+  {
+    return a.x == b.x && a.y == b.y;
+  }
+
+  friend bool operator!=(const double2 &a, const double2 &b)
+  {
+    return a.x != b.x || a.y != b.y;
+  }
+
+  friend std::ostream &operator<<(std::ostream &stream, const double2 &v)
+  {
+    stream << "(" << v.x << ", " << v.y << ")";
+    return stream;
+  }
+
+  static double dot(const double2 &a, const double2 &b)
+  {
+    return a.x * b.x + a.y * b.y;
+  }
+
+  static double2 interpolate(const double2 &a, const double2 &b, double t)
+  {
+    return a * (1 - t) + b * t;
+  }
+
+  static double2 abs(const double2 &a)
+  {
+    return double2(fabs(a.x), fabs(a.y));
+  }
+
+  static double distance(const double2 &a, const double2 &b)
+  {
+    return (a - b).length();
+  }
+
+  static double distance_squared(const double2 &a, const double2 &b)
+  {
+    double2 diff = a - b;
+    return double2::dot(diff, diff);
+  }
+
+  struct isect_result {
+    enum {
+      LINE_LINE_COLINEAR = -1,
+      LINE_LINE_NONE = 0,
+      LINE_LINE_EXACT = 1,
+      LINE_LINE_CROSS = 2,
+    } kind;
+    double lambda;
+  };
+
+  static isect_result isect_seg_seg(const double2 &v1,
+                                    const double2 &v2,
+                                    const double2 &v3,
+                                    const double2 &v4);
+};
+
+}  // namespace blender
diff --git a/source/blender/blenlib/BLI_double3.hh b/source/blender/blenlib/BLI_double3.hh
new file mode 100644
index 00000000000..ab258c9121b
--- /dev/null
+++ b/source/blender/blenlib/BLI_double3.hh
@@ -0,0 +1,246 @@
+/*
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+#pragma once
+
+/** \file
+ * \ingroup bli
+ */
+
+#include <iostream>
+
+#include "BLI_math_vector.h"
+#include "BLI_span.hh"
+
+namespace blender {
+
+struct double3 {
+  double x, y, z;
+
+  double3() = default;
+
+  double3(const double *ptr) : x{ptr[0]}, y{ptr[1]}, z{ptr[2]}
+  {
+  }
+
+  double3(const double (*ptr)[3]) : double3((const double *)ptr)
+  {
+  }
+
+  explicit double3(double value) : x(value), y(value), z(value)
+  {
+  }
+
+  explicit double3(int value) : x(value), y(value), z(value)
+  {
+  }
+
+  double3(double x, double y, double z) : x{x}, y{y}, z{z}
+  {
+  }
+
+  operator const double *() const
+  {
+    return &x;
+  }
+
+  operator double *()
+  {
+    return &x;
+  }
+
+  double normalize_and_get_length()
+  {
+    return normalize_v3_db(*this);
+  }
+
+  double3 normalized() const
+  {
+    double3 result;
+    normalize_v3_v3_db(result, *this);
+    return result;
+  }
+
+  double length() const
+  {
+    return len_v3_db(*this);
+  }
+
+  double length_squared() const
+  {
+    return len_squared_v3_db(*this);
+  }
+
+  void reflect(const double3 &normal)
+  {
+    *this = this->reflected(normal);
+  }
+
+  double3 reflected(const double3 &normal) const
+  {
+    double3 result;
+    reflect_v3_v3v3_db(result, *this, normal);
+    return result;
+  }
+
+  static double3 safe_divide(const double3 &a, const double3 &b)
+  {
+    double3 result;
+    result.x = (b.x == 0.0) ? 0.0 : a.x / b.x;
+    result.y = (b.y == 0.0) ? 0.0 : a.y / b.y;
+    result.z = (b.z == 0.0) ? 0.0 : a.z / b.z;
+    return result;
+  }
+
+  void invert()
+  {
+    x = -x;
+    y = -y;
+    z = -z;
+  }
+
+  friend double3 operator+(const double3 &a, const double3 &b)
+  {
+    return {a.x + b.x, a.y + b.y, a.z + b.z};
+  }
+
+  void operator+=(const double3 &b)
+  {
+    this->x += b.x;
+    this->y += b.y;
+    this->z += b.z;
+  }
+
+  friend double3 operator-(const double3 &a, const double3 &b)
+  {
+    return {a.x - b.x, a.y - b.y, a.z - b.z};
+  }
+
+  friend double3 operator-(const double3 &a)
+  {
+    return {-a.x, -a.y, -a.z};
+  }
+
+  void operator-=(const double3 &b)
+  {
+    this->x -= b.x;
+    this->y -= b.y;
+    this->z -= b.z;
+  }
+
+  void operator*=(const double &scalar)
+  {
+    this->x *= scalar;
+    this->y *= scalar;
+    this->z *= scalar;
+  }
+
+  void operator*=(const double3 &other)
+  {
+    this->x *= other.x;
+    this->y *= other.y;
+    this->z *= other.z;
+  }
+
+  friend double3 operator*(const double3 &a, const double3 &b)
+  {
+    return {a.x * b.x, a.y * b.y, a.z * b.z};
+  }
+
+  friend double3 operator*(const double3 &a, const double &b)
+  {
+    return {a.x * b, a.y * b, a.z * b};
+  }
+
+  friend double3 operator*(const double &a, const double3 &b)
+  {
+    return b * a;
+  }
+
+  friend double3 operator/(const double3 &a, const double &b)
+  {
+    BLI_assert(b != 0.0);
+    return {a.x / b, a.y / b, a.z / b};
+  }
+
+  friend bool operator==(const double3 &a, const double3 &b)
+  {
+    return a.x == b.x && a.y == b.y && a.z == b.z;
+  }
+
+  friend bool operator!=(const double3 &a, const double3 &b)
+  {
+    return a.x != b.x || a.y != b.y || a.z != b.z;
+  }
+
+  friend std::ostream &operator<<(std::ostream &stream, const double3 &v)
+  {
+    stream << "(" << v.x << ", " << v.y << ", " << v.z << ")";
+    return stream;
+  }
+
+  static double dot(const double3 &a, const double3 &b)
+  {
+    return a.x * b.x + a.y * b.y + a.z * b.z;
+  }
+
+  static double3 cross_high_precision(const double3 &a, const double3 &b)
+  {
+    double3 result;
+    cross_v3_v3v3_db(result, a, b);
+    return result;
+  }
+
+  static double3 project(const double3 &a, const double3 &b)
+  {
+    double3 result;
+    project_v3_v3v3_db(result, a, b);
+    return result;
+  }
+
+  static double distance(const double3 &a, const double3 &b)
+  {
+    return (a - b).length();
+  }
+
+  static double distance_squared(const double3 &a, const double3 &b)
+  {
+    double3 diff = a - b;
+    return double3::dot(diff, diff);
+  }
+
+  static double3 interpolate(const double3 &a, const double3 &b, double t)
+  {
+    return a * (1 - t) + b * t;
+  }
+
+  static double3 abs(const double3 &a)
+  {
+    return double3(fabs(a.x), fabs(a.y), fabs(a.z));
+  }
+
+  static int dominant_axis(const double3 &a)
+  {
+    double x = (a.x >= 0) ? a.x : -a.x;
+    double y = (a.y >= 0) ? a.y : -a.y;
+    double z = (a.z >= 0) ? a.z : -a.z;
+    return ((x > y) ? ((x > z) ? 0 : 2) : ((y > z) ? 1 : 2));
+  }
+
+  static double3 cross_poly(Span<double3> poly);
+};
+
+}  // namespace blender
diff --git a/source/blender/blenlib/BLI_fileops.h b/source/blender/blenlib/BLI_fileops.h
index 0022823b3de..28cb5f6d84b 100644
--- a/source/blender/blenlib/BLI_fileops.h
+++ b/source/blender/blenlib/BLI_fileops.h
@@ -155,8 +155,7 @@ double BLI_dir_free_space(const char *dir) ATTR_WARN_UNUSED_RESULT ATTR_NONNULL(
  *
  * \note can return NULL when the size is not big enough
  */
-char *BLI_current_working_dir(char *dir, const size_t maxncpy) ATTR_WARN_UNUSED_RESULT
-    ATTR_NONNULL();
+char *BLI_current_working_dir(char *dir, size_t maxncpy) ATTR_WARN_UNUSED_RESULT ATTR_NONNULL();
 eFileAttributes BLI_file_attributes(const char *path);
 
 /** \} */
diff --git a/source/blender/blenlib/BLI_float2.hh b/source/blender/blenlib/BLI_float2.hh
new file mode 100644
index 00000000000..bb4229db86e
--- /dev/null
+++ b/source/blender/blenlib/BLI_float2.hh
@@ -0,0 +1,218 @@
+/*
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+#pragma once
+
+#include "BLI_float3.hh"
+
+namespace blender {
+
+struct float2 {
+  float x, y;
+
+  float2() = default;
+
+  float2(const float *ptr) : x{ptr[0]}, y{ptr[1]}
+  {
+  }
+
+  explicit float2(float value) : x(value), y(value)
+  {
+  }
+
+  explicit float2(int value) : x(value), y(value)
+  {
+  }
+
+  float2(float x, float y) : x(x), y(y)
+  {
+  }
+
+  float2(const float3 &other) : x(other.x), y(other.y)
+  {
+  }
+
+  operator float *()
+  {
+    return &x;
+  }
+
+  operator const float *() const
+  {
+    return &x;
+  }
+
+  float length() const
+  {
+    return len_v2(*this);
+  }
+
+  float length_squared() const
+  {
+    return len_squared_v2(*this);
+  }
+
+  bool is_zero() const
+  {
+    return this->x == 0.0f && this->y == 0.0f;
+  }
+
+  float2 &operator+=(const float2 &other)
+  {
+    x += other.x;
+    y += other.y;
+    return *this;
+  }
+
+  float2 &operator-=(const float2 &other)
+  {
+    x -= other.x;
+    y -= other.y;
+    return *this;
+  }
+
+  float2 &operator*=(float factor)
+  {
+    x *= factor;
+    y *= factor;
+    return *this;
+  }
+
+  float2 &operator/=(float divisor)
+  {
+    x /= divisor;
+    y /= divisor;
+    return *this;
+  }
+
+  uint64_t hash() const
+  {
+    uint64_t x1 = *reinterpret_cast<const uint32_t *>(&x);
+    uint64_t x2 = *reinterpret_cast<const uint32_t *>(&y);
+    return (x1 * 812519) ^ (x2 * 707951);
+  }
+
+  friend float2 operator+(const float2 &a, const float2 &b)
+  {
+    return {a.x + b.x, a.y + b.y};
+  }
+
+  friend float2 operator-(const float2 &a, const float2 &b)
+  {
+    return {a.x - b.x, a.y - b.y};
+  }
+
+  friend float2 operator-(const float2 &a, const float &b)
+  {
+    return {a.x - b, a.y - b};
+  }
+
+  friend float2 operator*(const float2 &a, float b)
+  {
+    return {a.x * b, a.y * b};
+  }
+
+  friend float2 operator/(const float2 &a, float b)
+  {
+    BLI_assert(b != 0.0f);
+    return {a.x / b, a.y / b};
+  }
+
+  friend float2 operator*(float a, const float2 &b)
+  {
+    return b * a;
+  }
+
+  friend std::ostream &operator<<(std::ostream &stream, const float2 &v)
+  {
+    stream << "(" << v.x << ", " << v.y << ")";
+    return stream;
+  }
+
+  static float2 safe_divide(const float2 &a, const float b)
+  {
+    return (b != 0.0f) ? a / b : float2(0.0f);
+  }
+
+  static float2 floor(const float2 &a)
+  {
+    return float2(floorf(a.x), floorf(a.y));
+  }
+
+  /**
+   * Returns a normalized vector. The original vector is not changed.
+   */
+  float2 normalized() const
+  {
+    float2 result;
+    normalize_v2_v2(result, *this);
+    return result;
+  }
+
+  static float dot(const float2 &a, const float2 &b)
+  {
+    return a.x * b.x + a.y * b.y;
+  }
+
+  static float2 interpolate(const float2 &a, const float2 &b, float t)
+  {
+    return a * (1 - t) + b * t;
+  }
+
+  static float2 abs(const float2 &a)
+  {
+    return float2(fabsf(a.x), fabsf(a.y));
+  }
+
+  static float distance(const float2 &a, const float2 &b)
+  {
+    return (a - b).length();
+  }
+
+  static float distance_squared(const float2 &a, const float2 &b)
+  {
+    float2 diff = a - b;
+    return float2::dot(diff, diff);
+  }
+
+  struct isect_result {
+    enum {
+      LINE_LINE_COLINEAR = -1,
+      LINE_LINE_NONE = 0,
+      LINE_LINE_EXACT = 1,
+      LINE_LINE_CROSS = 2,
+    } kind;
+    float lambda;
+    float mu;
+  };
+
+  static isect_result isect_seg_seg(const float2 &v1,
+                                    const float2 &v2,
+                                    const float2 &v3,
+                                    const float2 &v4);
+
+  friend bool operator==(const float2 &a, const float2 &b)
+  {
+    return a.x == b.x && a.y == b.y;
+  }
+
+  friend bool operator!=(const float2 &a, const float2 &b)
+  {
+    return !(a == b);
+  }
+};
+
+}  // namespace blender
diff --git a/source/blender/blenlib/BLI_float3.hh b/source/blender/blenlib/BLI_float3.hh
new file mode 100644
index 00000000000..765f524fb31
--- /dev/null
+++ b/source/blender/blenlib/BLI_float3.hh
@@ -0,0 +1,320 @@
+/*
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+#pragma once
+
+#include <iostream>
+
+#include "BLI_math_vector.h"
+
+namespace blender {
+
+struct float3 {
+  float x, y, z;
+
+  float3() = default;
+
+  float3(const float *ptr) : x{ptr[0]}, y{ptr[1]}, z{ptr[2]}
+  {
+  }
+
+  float3(const float (*ptr)[3]) : float3(static_cast<const float *>(ptr[0]))
+  {
+  }
+
+  explicit float3(float value) : x(value), y(value), z(value)
+  {
+  }
+
+  explicit float3(int value) : x(value), y(value), z(value)
+  {
+  }
+
+  float3(float x, float y, float z) : x{x}, y{y}, z{z}
+  {
+  }
+
+  operator const float *() const
+  {
+    return &x;
+  }
+
+  operator float *()
+  {
+    return &x;
+  }
+
+  friend float3 operator+(const float3 &a, const float3 &b)
+  {
+    return {a.x + b.x, a.y + b.y, a.z + b.z};
+  }
+
+  friend float3 operator+(const float3 &a, const float &b)
+  {
+    return {a.x + b, a.y + b, a.z + b};
+  }
+
+  float3 &operator+=(const float3 &b)
+  {
+    this->x += b.x;
+    this->y += b.y;
+    this->z += b.z;
+    return *this;
+  }
+
+  friend float3 operator-(const float3 &a, const float3 &b)
+  {
+    return {a.x - b.x, a.y - b.y, a.z - b.z};
+  }
+
+  friend float3 operator-(const float3 &a)
+  {
+    return {-a.x, -a.y, -a.z};
+  }
+
+  friend float3 operator-(const float3 &a, const float &b)
+  {
+    return {a.x - b, a.y - b, a.z - b};
+  }
+
+  float3 &operator-=(const float3 &b)
+  {
+    this->x -= b.x;
+    this->y -= b.y;
+    this->z -= b.z;
+    return *this;
+  }
+
+  float3 &operator*=(float scalar)
+  {
+    this->x *= scalar;
+    this->y *= scalar;
+    this->z *= scalar;
+    return *this;
+  }
+
+  float3 &operator*=(const float3 &other)
+  {
+    this->x *= other.x;
+    this->y *= other.y;
+    this->z *= other.z;
+    return *this;
+  }
+
+  friend float3 operator*(const float3 &a, const float3 &b)
+  {
+    return {a.x * b.x, a.y * b.y, a.z * b.z};
+  }
+
+  friend float3 operator*(const float3 &a, float b)
+  {
+    return {a.x * b, a.y * b, a.z * b};
+  }
+
+  friend float3 operator*(float a, const float3 &b)
+  {
+    return b * a;
+  }
+
+  friend float3 operator/(const float3 &a, float b)
+  {
+    BLI_assert(b != 0.0f);
+    return {a.x / b, a.y / b, a.z / b};
+  }
+
+  friend std::ostream &operator<<(std::ostream &stream, const float3 &v)
+  {
+    stream << "(" << v.x << ", " << v.y << ", " << v.z << ")";
+    return stream;
+  }
+
+  friend bool operator==(const float3 &a, const float3 &b)
+  {
+    return a.x == b.x && a.y == b.y && a.z == b.z;
+  }
+
+  friend bool operator!=(const float3 &a, const float3 &b)
+  {
+    return !(a == b);
+  }
+
+  float normalize_and_get_length()
+  {
+    return normalize_v3(*this);
+  }
+
+  /**
+   * Normalizes the vector in place.
+   */
+  void normalize()
+  {
+    normalize_v3(*this);
+  }
+
+  /**
+   * Returns a normalized vector. The original vector is not changed.
+   */
+  float3 normalized() const
+  {
+    float3 result;
+    normalize_v3_v3(result, *this);
+    return result;
+  }
+
+  float length() const
+  {
+    return len_v3(*this);
+  }
+
+  float length_squared() const
+  {
+    return len_squared_v3(*this);
+  }
+
+  bool is_zero() const
+  {
+    return this->x == 0.0f && this->y == 0.0f && this->z == 0.0f;
+  }
+
+  void reflect(const float3 &normal)
+  {
+    *this = this->reflected(normal);
+  }
+
+  float3 reflected(const float3 &normal) const
+  {
+    float3 result;
+    reflect_v3_v3v3(result, *this, normal);
+    return result;
+  }
+
+  static float3 refract(const float3 &incident, const float3 &normal, const float eta)
+  {
+    float3 result;
+    float k = 1.0f - eta * eta * (1.0f - dot(normal, incident) * dot(normal, incident));
+    if (k < 0.0f) {
+      result = float3(0.0f);
+    }
+    else {
+      result = eta * incident - (eta * dot(normal, incident) + sqrt(k)) * normal;
+    }
+    return result;
+  }
+
+  static float3 faceforward(const float3 &vector, const float3 &incident, const float3 &reference)
+  {
+    return dot(reference, incident) < 0.0f ? vector : -vector;
+  }
+
+  static float3 safe_divide(const float3 &a, const float3 &b)
+  {
+    float3 result;
+    result.x = (b.x == 0.0f) ? 0.0f : a.x / b.x;
+    result.y = (b.y == 0.0f) ? 0.0f : a.y / b.y;
+    result.z = (b.z == 0.0f) ? 0.0f : a.z / b.z;
+    return result;
+  }
+
+  static float3 min(const float3 &a, const float3 &b)
+  {
+    return {a.x < b.x ? a.x : b.x, a.y < b.y ? a.y : b.y, a.z < b.z ? a.z : b.z};
+  }
+
+  static float3 max(const float3 &a, const float3 &b)
+  {
+    return {a.x > b.x ? a.x : b.x, a.y > b.y ? a.y : b.y, a.z > b.z ? a.z : b.z};
+  }
+
+  static void min_max(const float3 &vector, float3 &min, float3 &max)
+  {
+    min = float3::min(vector, min);
+    max = float3::max(vector, max);
+  }
+
+  static float3 safe_divide(const float3 &a, const float b)
+  {
+    return (b != 0.0f) ? a / b : float3(0.0f);
+  }
+
+  static float3 floor(const float3 &a)
+  {
+    return float3(floorf(a.x), floorf(a.y), floorf(a.z));
+  }
+
+  void invert()
+  {
+    x = -x;
+    y = -y;
+    z = -z;
+  }
+
+  uint64_t hash() const
+  {
+    uint64_t x1 = *reinterpret_cast<const uint32_t *>(&x);
+    uint64_t x2 = *reinterpret_cast<const uint32_t *>(&y);
+    uint64_t x3 = *reinterpret_cast<const uint32_t *>(&z);
+    return (x1 * 435109) ^ (x2 * 380867) ^ (x3 * 1059217);
+  }
+
+  static float dot(const float3 &a, const float3 &b)
+  {
+    return a.x * b.x + a.y * b.y + a.z * b.z;
+  }
+
+  static float3 cross_high_precision(const float3 &a, const float3 &b)
+  {
+    float3 result;
+    cross_v3_v3v3_hi_prec(result, a, b);
+    return result;
+  }
+
+  static float3 cross(const float3 &a, const float3 &b)
+  {
+    float3 result;
+    cross_v3_v3v3(result, a, b);
+    return result;
+  }
+
+  static float3 project(const float3 &a, const float3 &b)
+  {
+    float3 result;
+    project_v3_v3v3(result, a, b);
+    return result;
+  }
+
+  static float distance(const float3 &a, const float3 &b)
+  {
+    return (a - b).length();
+  }
+
+  static float distance_squared(const float3 &a, const float3 &b)
+  {
+    float3 diff = a - b;
+    return float3::dot(diff, diff);
+  }
+
+  static float3 interpolate(const float3 &a, const float3 &b, float t)
+  {
+    return a * (1 - t) + b * t;
+  }
+
+  static float3 abs(const float3 &a)
+  {
+    return float3(fabsf(a.x), fabsf(a.y), fabsf(a.z));
+  }
+};
+
+}  // namespace blender
diff --git a/source/blender/blenlib/BLI_float4.hh b/source/blender/blenlib/BLI_float4.hh
new file mode 100644
index 00000000000..5b487f6d029
--- /dev/null
+++ b/source/blender/blenlib/BLI_float4.hh
@@ -0,0 +1,138 @@
+/*
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+#pragma once
+
+namespace blender {
+
+struct float4 {
+  float x, y, z, w;
+
+  float4() = default;
+
+  float4(const float *ptr) : x{ptr[0]}, y{ptr[1]}, z{ptr[2]}, w{ptr[3]}
+  {
+  }
+
+  explicit float4(float value) : x(value), y(value), z(value), w(value)
+  {
+  }
+
+  explicit float4(int value) : x(value), y(value), z(value), w(value)
+  {
+  }
+
+  float4(float x, float y, float z, float w) : x(x), y(y), z(z), w(w)
+  {
+  }
+
+  operator float *()
+  {
+    return &x;
+  }
+
+  friend float4 operator+(const float4 &a, const float &b)
+  {
+    return {a.x + b, a.y + b, a.z + b, a.w + b};
+  }
+
+  operator const float *() const
+  {
+    return &x;
+  }
+
+  float4 &operator+=(const float4 &other)
+  {
+    x += other.x;
+    y += other.y;
+    z += other.z;
+    w += other.w;
+    return *this;
+  }
+
+  friend float4 operator-(const float4 &a, const float4 &b)
+  {
+    return {a.x - b.x, a.y - b.y, a.z - b.z, a.w - b.w};
+  }
+
+  friend float4 operator-(const float4 &a, const float &b)
+  {
+    return {a.x - b, a.y - b, a.z - b, a.w - b};
+  }
+
+  friend float4 operator+(const float4 &a, const float4 &b)
+  {
+    return {a.x + b.x, a.y + b.y, a.z + b.z, a.w + b.w};
+  }
+
+  friend float4 operator/(const float4 &a, float f)
+  {
+    BLI_assert(f != 0.0f);
+    return a * (1.0f / f);
+  }
+
+  float4 &operator*=(float factor)
+  {
+    x *= factor;
+    y *= factor;
+    z *= factor;
+    w *= factor;
+    return *this;
+  }
+
+  friend float4 operator*(const float4 &a, float b)
+  {
+    return {a.x * b, a.y * b, a.z * b, a.w * b};
+  }
+
+  friend float4 operator*(float a, const float4 &b)
+  {
+    return b * a;
+  }
+
+  float length() const
+  {
+    return len_v4(*this);
+  }
+
+  static float distance(const float4 &a, const float4 &b)
+  {
+    return (a - b).length();
+  }
+
+  static float4 safe_divide(const float4 &a, const float b)
+  {
+    return (b != 0.0f) ? a / b : float4(0.0f);
+  }
+
+  static float4 interpolate(const float4 &a, const float4 &b, float t)
+  {
+    return a * (1 - t) + b * t;
+  }
+
+  static float4 floor(const float4 &a)
+  {
+    return float4(floorf(a.x), floorf(a.y), floorf(a.z), floorf(a.w));
+  }
+
+  static float4 normalize(const float4 &a)
+  {
+    const float t = len_v4(a);
+    return (t != 0.0f) ? a / t : float4(0.0f);
+  }
+};
+
+}  // namespace blender
diff --git a/source/blender/blenlib/BLI_float4x4.hh b/source/blender/blenlib/BLI_float4x4.hh
index 81c969d02d0..b7f839f4ddf 100644
--- a/source/blender/blenlib/BLI_float4x4.hh
+++ b/source/blender/blenlib/BLI_float4x4.hh
@@ -16,9 +16,8 @@
 
 #pragma once
 
+#include "BLI_float3.hh"
 #include "BLI_math_matrix.h"
-#include "BLI_math_vec_types.hh"
-#include "BLI_math_vector.h"
 
 namespace blender {
 
@@ -64,7 +63,7 @@ struct float4x4 {
      * Without the negation, the result would be a so called improper rotation. That means it
      * contains a reflection. Such an improper rotation matrix could not be converted to another
      * representation of a rotation such as euler angles. */
-    const float3 cross = -math::cross(forward, up);
+    const float3 cross = -float3::cross(forward, up);
 
     float4x4 matrix;
     matrix.values[0][0] = forward.x;
diff --git a/source/blender/blenlib/BLI_gsqueue.h b/source/blender/blenlib/BLI_gsqueue.h
index 8b32c09b56b..a35c743c80b 100644
--- a/source/blender/blenlib/BLI_gsqueue.h
+++ b/source/blender/blenlib/BLI_gsqueue.h
@@ -31,7 +31,7 @@ extern "C" {
 
 typedef struct _GSQueue GSQueue;
 
-GSQueue *BLI_gsqueue_new(const size_t elem_size);
+GSQueue *BLI_gsqueue_new(size_t elem_size);
 /**
  * Returns true if the queue is empty, false otherwise.
  */
diff --git a/source/blender/blenlib/BLI_listbase.h b/source/blender/blenlib/BLI_listbase.h
index 64852b95ae4..f73d1f22502 100644
--- a/source/blender/blenlib/BLI_listbase.h
+++ b/source/blender/blenlib/BLI_listbase.h
@@ -85,7 +85,7 @@ void *BLI_findptr(const struct ListBase *listbase,
  */
 void *BLI_listbase_bytes_find(const ListBase *listbase,
                               const void *bytes,
-                              const size_t bytes_size,
+                              size_t bytes_size,
                               int offset) ATTR_WARN_UNUSED_RESULT ATTR_NONNULL(1, 2);
 /**
  * Find the first item in the list that matches the given string, or the given index as fallback.
@@ -96,7 +96,7 @@ void *BLI_listbase_bytes_find(const ListBase *listbase,
  */
 void *BLI_listbase_string_or_index_find(const struct ListBase *listbase,
                                         const char *string,
-                                        const size_t string_offset,
+                                        size_t string_offset,
                                         int index) ATTR_WARN_UNUSED_RESULT ATTR_NONNULL(1);
 
 /* Find backwards. */
@@ -133,7 +133,7 @@ void *BLI_rfindptr(const struct ListBase *listbase,
  */
 void *BLI_listbase_bytes_rfind(const ListBase *listbase,
                                const void *bytes,
-                               const size_t bytes_size,
+                               size_t bytes_size,
                                int offset) ATTR_WARN_UNUSED_RESULT ATTR_NONNULL(1, 2);
 
 /**
diff --git a/source/blender/blenlib/BLI_math_boolean.hh b/source/blender/blenlib/BLI_math_boolean.hh
index 8cf93c82dec..20fd00b2aa4 100644
--- a/source/blender/blenlib/BLI_math_boolean.hh
+++ b/source/blender/blenlib/BLI_math_boolean.hh
@@ -21,11 +21,13 @@
  * \brief Math vector functions needed specifically for mesh intersect and boolean.
  */
 
-#include "BLI_math_vec_types.hh"
+#include "BLI_double2.hh"
+#include "BLI_double3.hh"
 
 #ifdef WITH_GMP
 #  include "BLI_math_mpq.hh"
-#  include "BLI_math_vec_mpq_types.hh"
+#  include "BLI_mpq2.hh"
+#  include "BLI_mpq3.hh"
 #endif
 
 namespace blender {
diff --git a/source/blender/blenlib/BLI_math_vec_mpq_types.hh b/source/blender/blenlib/BLI_math_vec_mpq_types.hh
deleted file mode 100644
index 36eb0cac83c..00000000000
--- a/source/blender/blenlib/BLI_math_vec_mpq_types.hh
+++ /dev/null
@@ -1,91 +0,0 @@
-/*
- * This program is free software; you can redistribute it and/or
- * modify it under the terms of the GNU General Public License
- * as published by the Free Software Foundation; either version 2
- * of the License, or (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program; if not, write to the Free Software Foundation,
- * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
- */
-
-#pragma once
-
-/** \file
- * \ingroup bli
- */
-
-#ifdef WITH_GMP
-
-#  include "BLI_math_mpq.hh"
-#  include "BLI_math_vec_types.hh"
-
-namespace blender {
-
-using mpq2 = vec_base<mpq_class, 2>;
-using mpq3 = vec_base<mpq_class, 3>;
-
-namespace math {
-
-uint64_t hash_mpq_class(const mpq_class &value);
-
-template<> inline uint64_t vector_hash(const mpq2 &vec)
-{
-  return hash_mpq_class(vec.x) ^ (hash_mpq_class(vec.y) * 33);
-}
-
-template<> inline uint64_t vector_hash(const mpq3 &vec)
-{
-  return hash_mpq_class(vec.x) ^ (hash_mpq_class(vec.y) * 33) ^ (hash_mpq_class(vec.z) * 33 * 37);
-}
-
-/**
- * Cannot do this exactly in rational arithmetic!
- * Approximate by going in and out of doubles.
- */
-template<> inline mpq_class length(const mpq2 &a)
-{
-  return mpq_class(sqrt(length_squared(a).get_d()));
-}
-
-/**
- * Cannot do this exactly in rational arithmetic!
- * Approximate by going in and out of doubles.
- */
-template<> inline mpq_class length(const mpq3 &a)
-{
-  return mpq_class(sqrt(length_squared(a).get_d()));
-}
-
-/**
- * The buffer avoids allocating a temporary variable.
- */
-inline mpq_class distance_squared_with_buffer(const mpq3 &a, const mpq3 &b, mpq3 &buffer)
-{
-  buffer = a;
-  buffer -= b;
-  return dot(buffer, buffer);
-}
-
-/**
- * The buffer avoids allocating a temporary variable.
- */
-inline mpq_class dot_with_buffer(const mpq3 &a, const mpq3 &b, mpq3 &buffer)
-{
-  buffer = a;
-  buffer *= b;
-  buffer.x += buffer.y;
-  buffer.x += buffer.z;
-  return buffer.x;
-}
-
-}  // namespace math
-
-}  // namespace blender
-
-#endif /* WITH_GMP */
diff --git a/source/blender/blenlib/BLI_math_vec_types.hh b/source/blender/blenlib/BLI_math_vec_types.hh
deleted file mode 100644
index 52aacd294e4..00000000000
--- a/source/blender/blenlib/BLI_math_vec_types.hh
+++ /dev/null
@@ -1,566 +0,0 @@
-/*
- * This program is free software; you can redistribute it and/or
- * modify it under the terms of the GNU General Public License
- * as published by the Free Software Foundation; either version 2
- * of the License, or (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program; if not, write to the Free Software Foundation,
- * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
- *
- * Copyright 2022, Blender Foundation.
- */
-
-#pragma once
-
-/** \file
- * \ingroup bli
- */
-
-#include <array>
-#include <cmath>
-#include <iostream>
-#include <type_traits>
-
-#include "BLI_math_vector.hh"
-#include "BLI_utildefines.h"
-
-namespace blender {
-
-/* clang-format off */
-template<typename T>
-using as_uint_type = std::conditional_t<sizeof(T) == sizeof(uint8_t), uint8_t,
-                     std::conditional_t<sizeof(T) == sizeof(uint16_t), uint16_t,
-                     std::conditional_t<sizeof(T) == sizeof(uint32_t), uint32_t,
-                     std::conditional_t<sizeof(T) == sizeof(uint64_t), uint64_t, void>>>>;
-/* clang-format on */
-
-template<typename T, int Size> struct vec_struct_base {
-  std::array<T, Size> values;
-};
-
-template<typename T> struct vec_struct_base<T, 2> {
-  T x, y;
-};
-
-template<typename T> struct vec_struct_base<T, 3> {
-  T x, y, z;
-};
-
-template<typename T> struct vec_struct_base<T, 4> {
-  T x, y, z, w;
-};
-
-template<typename T, int Size> struct vec_base : public vec_struct_base<T, Size> {
-
-  static constexpr int type_length = Size;
-
-  using base_type = T;
-  using uint_type = vec_base<as_uint_type<T>, Size>;
-
-  vec_base() = default;
-
-  explicit vec_base(uint value)
-  {
-    for (int i = 0; i < Size; i++) {
-      (*this)[i] = static_cast<T>(value);
-    }
-  }
-
-  explicit vec_base(int value)
-  {
-    for (int i = 0; i < Size; i++) {
-      (*this)[i] = static_cast<T>(value);
-    }
-  }
-
-  explicit vec_base(float value)
-  {
-    for (int i = 0; i < Size; i++) {
-      (*this)[i] = static_cast<T>(value);
-    }
-  }
-
-  explicit vec_base(double value)
-  {
-    for (int i = 0; i < Size; i++) {
-      (*this)[i] = static_cast<T>(value);
-    }
-  }
-
-/* Workaround issue with template BLI_ENABLE_IF((Size == 2)) not working. */
-#define BLI_ENABLE_IF_VEC(_size, _test) int S = _size, BLI_ENABLE_IF((S _test))
-
-  template<BLI_ENABLE_IF_VEC(Size, == 2)> vec_base(T _x, T _y)
-  {
-    (*this)[0] = _x;
-    (*this)[1] = _y;
-  }
-
-  template<BLI_ENABLE_IF_VEC(Size, == 3)> vec_base(T _x, T _y, T _z)
-  {
-    (*this)[0] = _x;
-    (*this)[1] = _y;
-    (*this)[2] = _z;
-  }
-
-  template<BLI_ENABLE_IF_VEC(Size, == 4)> vec_base(T _x, T _y, T _z, T _w)
-  {
-    (*this)[0] = _x;
-    (*this)[1] = _y;
-    (*this)[2] = _z;
-    (*this)[3] = _w;
-  }
-
-  /** Mixed scalar-vector constructors. */
-
-  template<typename U, BLI_ENABLE_IF_VEC(Size, == 3)>
-  constexpr vec_base(const vec_base<U, 2> &xy, T z)
-      : vec_base(static_cast<T>(xy.x), static_cast<T>(xy.y), z)
-  {
-  }
-
-  template<typename U, BLI_ENABLE_IF_VEC(Size, == 3)>
-  constexpr vec_base(T x, const vec_base<U, 2> &yz)
-      : vec_base(x, static_cast<T>(yz.x), static_cast<T>(yz.y))
-  {
-  }
-
-  template<typename U, BLI_ENABLE_IF_VEC(Size, == 4)>
-  vec_base(vec_base<U, 3> xyz, T w)
-      : vec_base(
-            static_cast<T>(xyz.x), static_cast<T>(xyz.y), static_cast<T>(xyz.z), static_cast<T>(w))
-  {
-  }
-
-  template<typename U, BLI_ENABLE_IF_VEC(Size, == 4)>
-  vec_base(T x, vec_base<U, 3> yzw)
-      : vec_base(
-            static_cast<T>(x), static_cast<T>(yzw.x), static_cast<T>(yzw.y), static_cast<T>(yzw.z))
-  {
-  }
-
-  template<typename U, typename V, BLI_ENABLE_IF_VEC(Size, == 4)>
-  vec_base(vec_base<U, 2> xy, vec_base<V, 2> zw)
-      : vec_base(
-            static_cast<T>(xy.x), static_cast<T>(xy.y), static_cast<T>(zw.x), static_cast<T>(zw.y))
-  {
-  }
-
-  template<typename U, BLI_ENABLE_IF_VEC(Size, == 4)>
-  vec_base(vec_base<U, 2> xy, T z, T w)
-      : vec_base(static_cast<T>(xy.x), static_cast<T>(xy.y), static_cast<T>(z), static_cast<T>(w))
-  {
-  }
-
-  template<typename U, BLI_ENABLE_IF_VEC(Size, == 4)>
-  vec_base(T x, vec_base<U, 2> yz, T w)
-      : vec_base(static_cast<T>(x), static_cast<T>(yz.x), static_cast<T>(yz.y), static_cast<T>(w))
-  {
-  }
-
-  template<typename U, BLI_ENABLE_IF_VEC(Size, == 4)>
-  vec_base(T x, T y, vec_base<U, 2> zw)
-      : vec_base(static_cast<T>(x), static_cast<T>(y), static_cast<T>(zw.x), static_cast<T>(zw.y))
-  {
-  }
-
-  /** Masking. */
-
-  template<typename U, int OtherSize, BLI_ENABLE_IF(OtherSize > Size)>
-  explicit vec_base(const vec_base<U, OtherSize> &other)
-  {
-    for (int i = 0; i < Size; i++) {
-      (*this)[i] = static_cast<T>(other[i]);
-    }
-  }
-
-#undef BLI_ENABLE_IF_VEC
-
-  /** Conversion from pointers (from C-style vectors). */
-
-  vec_base(const T *ptr)
-  {
-    for (int i = 0; i < Size; i++) {
-      (*this)[i] = ptr[i];
-    }
-  }
-
-  vec_base(const T (*ptr)[Size]) : vec_base(static_cast<const T *>(ptr[0]))
-  {
-  }
-
-  /** Conversion from other vector types. */
-
-  template<typename U> explicit vec_base(const vec_base<U, Size> &vec)
-  {
-    for (int i = 0; i < Size; i++) {
-      (*this)[i] = static_cast<T>(vec[i]);
-    }
-  }
-
-  /** C-style pointer dereference. */
-
-  operator const T *() const
-  {
-    return reinterpret_cast<const T *>(this);
-  }
-
-  operator T *()
-  {
-    return reinterpret_cast<T *>(this);
-  }
-
-  /** Array access. */
-
-  const T &operator[](int index) const
-  {
-    BLI_assert(index >= 0);
-    BLI_assert(index < Size);
-    return reinterpret_cast<const T *>(this)[index];
-  }
-
-  T &operator[](int index)
-  {
-    BLI_assert(index >= 0);
-    BLI_assert(index < Size);
-    return reinterpret_cast<T *>(this)[index];
-  }
-
-  /** Internal Operators Macro. */
-
-#define BLI_INT_OP(_T) template<typename U = _T, BLI_ENABLE_IF((std::is_integral_v<U>))>
-
-#define BLI_VEC_OP_IMPL(_result, _i, _op) \
-  vec_base _result; \
-  for (int _i = 0; _i < Size; _i++) { \
-    _op; \
-  } \
-  return _result;
-
-#define BLI_VEC_OP_IMPL_SELF(_i, _op) \
-  for (int _i = 0; _i < Size; _i++) { \
-    _op; \
-  } \
-  return *this;
-
-  /** Arithmetic operators. */
-
-  friend vec_base operator+(const vec_base &a, const vec_base &b)
-  {
-    BLI_VEC_OP_IMPL(ret, i, ret[i] = a[i] + b[i]);
-  }
-
-  friend vec_base operator+(const vec_base &a, const T &b)
-  {
-    BLI_VEC_OP_IMPL(ret, i, ret[i] = a[i] + b);
-  }
-
-  friend vec_base operator+(const T &a, const vec_base &b)
-  {
-    return b + a;
-  }
-
-  vec_base &operator+=(const vec_base &b)
-  {
-    BLI_VEC_OP_IMPL_SELF(i, (*this)[i] += b[i]);
-  }
-
-  vec_base &operator+=(const T &b)
-  {
-    BLI_VEC_OP_IMPL_SELF(i, (*this)[i] += b);
-  }
-
-  friend vec_base operator-(const vec_base &a)
-  {
-    BLI_VEC_OP_IMPL(ret, i, ret[i] = -a[i]);
-  }
-
-  friend vec_base operator-(const vec_base &a, const vec_base &b)
-  {
-    BLI_VEC_OP_IMPL(ret, i, ret[i] = a[i] - b[i]);
-  }
-
-  friend vec_base operator-(const vec_base &a, const T &b)
-  {
-    BLI_VEC_OP_IMPL(ret, i, ret[i] = a[i] - b);
-  }
-
-  friend vec_base operator-(const T &a, const vec_base &b)
-  {
-    BLI_VEC_OP_IMPL(ret, i, ret[i] = a - b[i]);
-  }
-
-  vec_base &operator-=(const vec_base &b)
-  {
-    BLI_VEC_OP_IMPL_SELF(i, (*this)[i] -= b[i]);
-  }
-
-  vec_base &operator-=(const T &b)
-  {
-    BLI_VEC_OP_IMPL_SELF(i, (*this)[i] -= b);
-  }
-
-  friend vec_base operator*(const vec_base &a, const vec_base &b)
-  {
-    BLI_VEC_OP_IMPL(ret, i, ret[i] = a[i] * b[i]);
-  }
-
-  friend vec_base operator*(const vec_base &a, T b)
-  {
-    BLI_VEC_OP_IMPL(ret, i, ret[i] = a[i] * b);
-  }
-
-  friend vec_base operator*(T a, const vec_base &b)
-  {
-    return b * a;
-  }
-
-  vec_base &operator*=(T b)
-  {
-    BLI_VEC_OP_IMPL_SELF(i, (*this)[i] *= b);
-  }
-
-  vec_base &operator*=(const vec_base &b)
-  {
-    BLI_VEC_OP_IMPL_SELF(i, (*this)[i] *= b[i]);
-  }
-
-  friend vec_base operator/(const vec_base &a, const vec_base &b)
-  {
-    BLI_assert(!math::is_any_zero(b));
-    BLI_VEC_OP_IMPL(ret, i, ret[i] = a[i] / b[i]);
-  }
-
-  friend vec_base operator/(const vec_base &a, T b)
-  {
-    BLI_assert(b != T(0));
-    BLI_VEC_OP_IMPL(ret, i, ret[i] = a[i] / b);
-  }
-
-  friend vec_base operator/(T a, const vec_base &b)
-  {
-    BLI_assert(!math::is_any_zero(b));
-    BLI_VEC_OP_IMPL(ret, i, ret[i] = a / b[i]);
-  }
-
-  vec_base &operator/=(T b)
-  {
-    BLI_assert(b != T(0));
-    BLI_VEC_OP_IMPL_SELF(i, (*this)[i] /= b);
-  }
-
-  vec_base &operator/=(const vec_base &b)
-  {
-    BLI_assert(!math::is_any_zero(b));
-    BLI_VEC_OP_IMPL_SELF(i, (*this)[i] /= b[i]);
-  }
-
-  /** Binary operators. */
-
-  BLI_INT_OP(T) friend vec_base operator&(const vec_base &a, const vec_base &b)
-  {
-    BLI_VEC_OP_IMPL(ret, i, ret[i] = a[i] & b[i]);
-  }
-
-  BLI_INT_OP(T) friend vec_base operator&(const vec_base &a, T b)
-  {
-    BLI_VEC_OP_IMPL(ret, i, ret[i] = a[i] & b);
-  }
-
-  BLI_INT_OP(T) friend vec_base operator&(T a, const vec_base &b)
-  {
-    return b & a;
-  }
-
-  BLI_INT_OP(T) vec_base &operator&=(T b)
-  {
-    BLI_VEC_OP_IMPL_SELF(i, (*this)[i] &= b);
-  }
-
-  BLI_INT_OP(T) vec_base &operator&=(const vec_base &b)
-  {
-    BLI_VEC_OP_IMPL_SELF(i, (*this)[i] &= b[i]);
-  }
-
-  BLI_INT_OP(T) friend vec_base operator|(const vec_base &a, const vec_base &b)
-  {
-    BLI_VEC_OP_IMPL(ret, i, ret[i] = a[i] | b[i]);
-  }
-
-  BLI_INT_OP(T) friend vec_base operator|(const vec_base &a, T b)
-  {
-    BLI_VEC_OP_IMPL(ret, i, ret[i] = a[i] | b);
-  }
-
-  BLI_INT_OP(T) friend vec_base operator|(T a, const vec_base &b)
-  {
-    return b | a;
-  }
-
-  BLI_INT_OP(T) vec_base &operator|=(T b)
-  {
-    BLI_VEC_OP_IMPL_SELF(i, (*this)[i] |= b);
-  }
-
-  BLI_INT_OP(T) vec_base &operator|=(const vec_base &b)
-  {
-    BLI_VEC_OP_IMPL_SELF(i, (*this)[i] |= b[i]);
-  }
-
-  BLI_INT_OP(T) friend vec_base operator^(const vec_base &a, const vec_base &b)
-  {
-    BLI_VEC_OP_IMPL(ret, i, ret[i] = a[i] ^ b[i]);
-  }
-
-  BLI_INT_OP(T) friend vec_base operator^(const vec_base &a, T b)
-  {
-    BLI_VEC_OP_IMPL(ret, i, ret[i] = a[i] ^ b);
-  }
-
-  BLI_INT_OP(T) friend vec_base operator^(T a, const vec_base &b)
-  {
-    return b ^ a;
-  }
-
-  BLI_INT_OP(T) vec_base &operator^=(T b)
-  {
-    BLI_VEC_OP_IMPL_SELF(i, (*this)[i] ^= b);
-  }
-
-  BLI_INT_OP(T) vec_base &operator^=(const vec_base &b)
-  {
-    BLI_VEC_OP_IMPL_SELF(i, (*this)[i] ^= b[i]);
-  }
-
-  BLI_INT_OP(T) friend vec_base operator~(const vec_base &a)
-  {
-    BLI_VEC_OP_IMPL(ret, i, ret[i] = ~a[i]);
-  }
-
-  /** Bit-shift operators. */
-
-  BLI_INT_OP(T) friend vec_base operator<<(const vec_base &a, const vec_base &b)
-  {
-    BLI_VEC_OP_IMPL(ret, i, ret[i] = a[i] << b[i]);
-  }
-
-  BLI_INT_OP(T) friend vec_base operator<<(const vec_base &a, T b)
-  {
-    BLI_VEC_OP_IMPL(ret, i, ret[i] = a[i] << b);
-  }
-
-  BLI_INT_OP(T) vec_base &operator<<=(T b)
-  {
-    BLI_VEC_OP_IMPL_SELF(i, (*this)[i] <<= b);
-  }
-
-  BLI_INT_OP(T) vec_base &operator<<=(const vec_base &b)
-  {
-    BLI_VEC_OP_IMPL_SELF(i, (*this)[i] <<= b[i]);
-  }
-
-  BLI_INT_OP(T) friend vec_base operator>>(const vec_base &a, const vec_base &b)
-  {
-    BLI_VEC_OP_IMPL(ret, i, ret[i] = a[i] >> b[i]);
-  }
-
-  BLI_INT_OP(T) friend vec_base operator>>(const vec_base &a, T b)
-  {
-    BLI_VEC_OP_IMPL(ret, i, ret[i] = a[i] >> b);
-  }
-
-  BLI_INT_OP(T) vec_base &operator>>=(T b)
-  {
-    BLI_VEC_OP_IMPL_SELF(i, (*this)[i] >>= b);
-  }
-
-  BLI_INT_OP(T) vec_base &operator>>=(const vec_base &b)
-  {
-    BLI_VEC_OP_IMPL_SELF(i, (*this)[i] >>= b[i]);
-  }
-
-  /** Modulo operators. */
-
-  BLI_INT_OP(T) friend vec_base operator%(const vec_base &a, const vec_base &b)
-  {
-    BLI_assert(!math::is_any_zero(b));
-    BLI_VEC_OP_IMPL(ret, i, ret[i] = a[i] % b[i]);
-  }
-
-  BLI_INT_OP(T) friend vec_base operator%(const vec_base &a, T b)
-  {
-    BLI_assert(b != 0);
-    BLI_VEC_OP_IMPL(ret, i, ret[i] = a[i] % b);
-  }
-
-  BLI_INT_OP(T) friend vec_base operator%(T a, const vec_base &b)
-  {
-    BLI_assert(!math::is_any_zero(b));
-    BLI_VEC_OP_IMPL(ret, i, ret[i] = a % b[i]);
-  }
-
-#undef BLI_INT_OP
-#undef BLI_VEC_OP_IMPL
-#undef BLI_VEC_OP_IMPL_SELF
-
-  /** Compare. */
-
-  friend bool operator==(const vec_base &a, const vec_base &b)
-  {
-    for (int i = 0; i < Size; i++) {
-      if (a[i] != b[i]) {
-        return false;
-      }
-    }
-    return true;
-  }
-
-  friend bool operator!=(const vec_base &a, const vec_base &b)
-  {
-    return !(a == b);
-  }
-
-  /** Misc. */
-
-  uint64_t hash() const
-  {
-    return math::vector_hash(*this);
-  }
-
-  friend std::ostream &operator<<(std::ostream &stream, const vec_base &v)
-  {
-    stream << "(";
-    for (int i = 0; i < Size; i++) {
-      stream << v[i];
-      if (i != Size - 1) {
-        stream << ", ";
-      }
-    }
-    stream << ")";
-    return stream;
-  }
-};
-
-using int2 = vec_base<int32_t, 2>;
-using int3 = vec_base<int32_t, 3>;
-using int4 = vec_base<int32_t, 4>;
-
-using uint2 = vec_base<uint32_t, 2>;
-using uint3 = vec_base<uint32_t, 3>;
-using uint4 = vec_base<uint32_t, 4>;
-
-using float2 = vec_base<float, 2>;
-using float3 = vec_base<float, 3>;
-using float4 = vec_base<float, 4>;
-
-using double2 = vec_base<double, 2>;
-using double3 = vec_base<double, 3>;
-using double4 = vec_base<double, 4>;
-
-}  // namespace blender
diff --git a/source/blender/blenlib/BLI_math_vector.hh b/source/blender/blenlib/BLI_math_vector.hh
deleted file mode 100644
index e7d765df842..00000000000
--- a/source/blender/blenlib/BLI_math_vector.hh
+++ /dev/null
@@ -1,399 +0,0 @@
-/*
- * This program is free software; you can redistribute it and/or
- * modify it under the terms of the GNU General Public License
- * as published by the Free Software Foundation; either version 2
- * of the License, or (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program; if not, write to the Free Software Foundation,
- * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
- *
- * Copyright 2022, Blender Foundation.
- */
-
-#pragma once
-
-/** \file
- * \ingroup bli
- */
-
-#include <cmath>
-#include <type_traits>
-
-#include "BLI_math_base_safe.h"
-#include "BLI_math_vector.h"
-#include "BLI_span.hh"
-#include "BLI_utildefines.h"
-
-#ifdef WITH_GMP
-#  include "BLI_math_mpq.hh"
-#endif
-
-namespace blender::math {
-
-#ifndef NDEBUG
-#  define BLI_ASSERT_UNIT(v) \
-    { \
-      const float _test_unit = length_squared(v); \
-      BLI_assert(!(std::abs(_test_unit - 1.0f) >= BLI_ASSERT_UNIT_EPSILON) || \
-                 !(std::abs(_test_unit) >= BLI_ASSERT_UNIT_EPSILON)); \
-    } \
-    (void)0
-#else
-#  define BLI_ASSERT_UNIT(v) (void)(v)
-#endif
-
-#define bT typename T::base_type
-
-#ifdef WITH_GMP
-#  define BLI_ENABLE_IF_FLT_VEC(T) \
-    BLI_ENABLE_IF((std::is_floating_point_v<typename T::base_type> || \
-                   std::is_same_v<typename T::base_type, mpq_class>))
-#else
-#  define BLI_ENABLE_IF_FLT_VEC(T) BLI_ENABLE_IF((std::is_floating_point_v<typename T::base_type>))
-#endif
-
-#define BLI_ENABLE_IF_INT_VEC(T) BLI_ENABLE_IF((std::is_integral_v<typename T::base_type>))
-
-template<typename T> inline bool is_zero(const T &a)
-{
-  for (int i = 0; i < T::type_length; i++) {
-    if (a[i] != bT(0)) {
-      return false;
-    }
-  }
-  return true;
-}
-
-template<typename T> inline bool is_any_zero(const T &a)
-{
-  for (int i = 0; i < T::type_length; i++) {
-    if (a[i] == bT(0)) {
-      return true;
-    }
-  }
-  return false;
-}
-
-template<typename T> inline T abs(const T &a)
-{
-  T result;
-  for (int i = 0; i < T::type_length; i++) {
-    result[i] = a[i] >= 0 ? a[i] : -a[i];
-  }
-  return result;
-}
-
-template<typename T> inline T min(const T &a, const T &b)
-{
-  T result;
-  for (int i = 0; i < T::type_length; i++) {
-    result[i] = a[i] < b[i] ? a[i] : b[i];
-  }
-  return result;
-}
-
-template<typename T> inline T max(const T &a, const T &b)
-{
-  T result;
-  for (int i = 0; i < T::type_length; i++) {
-    result[i] = a[i] > b[i] ? a[i] : b[i];
-  }
-  return result;
-}
-
-template<typename T> inline T clamp(const T &a, const T &min_v, const T &max_v)
-{
-  T result = a;
-  for (int i = 0; i < T::type_length; i++) {
-    CLAMP(result[i], min_v[i], max_v[i]);
-  }
-  return result;
-}
-
-template<typename T> inline T clamp(const T &a, const bT &min_v, const bT &max_v)
-{
-  T result = a;
-  for (int i = 0; i < T::type_length; i++) {
-    CLAMP(result[i], min_v, max_v);
-  }
-  return result;
-}
-
-template<typename T, BLI_ENABLE_IF_FLT_VEC(T)> inline T mod(const T &a, const T &b)
-{
-  T result;
-  for (int i = 0; i < T::type_length; i++) {
-    BLI_assert(b[i] != 0);
-    result[i] = std::fmod(a[i], b[i]);
-  }
-  return result;
-}
-
-template<typename T, BLI_ENABLE_IF_FLT_VEC(T)> inline T mod(const T &a, bT b)
-{
-  BLI_assert(b != 0);
-  T result;
-  for (int i = 0; i < T::type_length; i++) {
-    result[i] = std::fmod(a[i], b);
-  }
-  return result;
-}
-
-template<typename T, BLI_ENABLE_IF_FLT_VEC(T)> inline T safe_mod(const T &a, const T &b)
-{
-  T result;
-  for (int i = 0; i < T::type_length; i++) {
-    result[i] = (b[i] != 0) ? std::fmod(a[i], b[i]) : 0;
-  }
-  return result;
-}
-
-template<typename T, BLI_ENABLE_IF_FLT_VEC(T)> inline T safe_mod(const T &a, bT b)
-{
-  if (b == 0) {
-    return T(0.0f);
-  }
-  T result;
-  for (int i = 0; i < T::type_length; i++) {
-    result[i] = std::fmod(a[i], b);
-  }
-  return result;
-}
-
-template<typename T> inline void min_max(const T &vector, T &min_vec, T &max_vec)
-{
-  min_vec = min(vector, min_vec);
-  max_vec = max(vector, max_vec);
-}
-
-template<typename T, BLI_ENABLE_IF_FLT_VEC(T)> inline T safe_divide(const T &a, const T &b)
-{
-  T result;
-  for (int i = 0; i < T::type_length; i++) {
-    result[i] = (b[i] == 0) ? 0 : a[i] / b[i];
-  }
-  return result;
-}
-
-template<typename T, BLI_ENABLE_IF_FLT_VEC(T)> inline T safe_divide(const T &a, const bT b)
-{
-  return (b != 0) ? a / b : T(0.0f);
-}
-
-template<typename T, BLI_ENABLE_IF_FLT_VEC(T)> inline T floor(const T &a)
-{
-  T result;
-  for (int i = 0; i < T::type_length; i++) {
-    result[i] = std::floor(a[i]);
-  }
-  return result;
-}
-
-template<typename T, BLI_ENABLE_IF_FLT_VEC(T)> inline T ceil(const T &a)
-{
-  T result;
-  for (int i = 0; i < T::type_length; i++) {
-    result[i] = std::ceil(a[i]);
-  }
-  return result;
-}
-
-template<typename T, BLI_ENABLE_IF_FLT_VEC(T)> inline T fract(const T &a)
-{
-  T result;
-  for (int i = 0; i < T::type_length; i++) {
-    result[i] = a[i] - std::floor(a[i]);
-  }
-  return result;
-}
-
-template<typename T, BLI_ENABLE_IF_FLT_VEC(T)> inline bT dot(const T &a, const T &b)
-{
-  bT result = a[0] * b[0];
-  for (int i = 1; i < T::type_length; i++) {
-    result += a[i] * b[i];
-  }
-  return result;
-}
-
-template<typename T> inline bT length_manhattan(const T &a)
-{
-  bT result = std::abs(a[0]);
-  for (int i = 1; i < T::type_length; i++) {
-    result += std::abs(a[i]);
-  }
-  return result;
-}
-
-template<typename T, BLI_ENABLE_IF_FLT_VEC(T)> inline bT length_squared(const T &a)
-{
-  return dot(a, a);
-}
-
-template<typename T, BLI_ENABLE_IF_FLT_VEC(T)> inline bT length(const T &a)
-{
-  return std::sqrt(length_squared(a));
-}
-
-template<typename T, BLI_ENABLE_IF_FLT_VEC(T)> inline bT distance_manhattan(const T &a, const T &b)
-{
-  return length_manhattan(a - b);
-}
-
-template<typename T, BLI_ENABLE_IF_FLT_VEC(T)> inline bT distance_squared(const T &a, const T &b)
-{
-  return length_squared(a - b);
-}
-
-template<typename T, BLI_ENABLE_IF_FLT_VEC(T)> inline bT distance(const T &a, const T &b)
-{
-  return length(a - b);
-}
-
-template<typename T> uint64_t vector_hash(const T &vec)
-{
-  BLI_STATIC_ASSERT(T::type_length <= 4, "Longer types need to implement vector_hash themself.");
-  const typename T::uint_type &uvec = *reinterpret_cast<const typename T::uint_type *>(&vec);
-  uint64_t result;
-  result = uvec[0] * uint64_t(435109);
-  if constexpr (T::type_length > 1) {
-    result ^= uvec[1] * uint64_t(380867);
-  }
-  if constexpr (T::type_length > 2) {
-    result ^= uvec[2] * uint64_t(1059217);
-  }
-  if constexpr (T::type_length > 3) {
-    result ^= uvec[3] * uint64_t(2002613);
-  }
-  return result;
-}
-
-template<typename T, BLI_ENABLE_IF_FLT_VEC(T)> inline T reflect(const T &incident, const T &normal)
-{
-  BLI_ASSERT_UNIT(normal);
-  return incident - 2.0 * dot(normal, incident) * normal;
-}
-
-template<typename T, BLI_ENABLE_IF_FLT_VEC(T)>
-inline T refract(const T &incident, const T &normal, const bT eta)
-{
-  float dot_ni = dot(normal, incident);
-  float k = 1.0f - eta * eta * (1.0f - dot_ni * dot_ni);
-  if (k < 0.0f) {
-    return T(0.0f);
-  }
-  return eta * incident - (eta * dot_ni + sqrt(k)) * normal;
-}
-
-template<typename T, BLI_ENABLE_IF_FLT_VEC(T)> inline T project(const T &p, const T &v_proj)
-{
-  if (UNLIKELY(is_zero(v_proj))) {
-    return T(0.0f);
-  }
-  return v_proj * (dot(p, v_proj) / dot(v_proj, v_proj));
-}
-
-template<typename T, BLI_ENABLE_IF_FLT_VEC(T)>
-inline T normalize_and_get_length(const T &v, bT &out_length)
-{
-  out_length = length_squared(v);
-  /* A larger value causes normalize errors in a scaled down models with camera extreme close. */
-  constexpr bT threshold = std::is_same_v<bT, double> ? 1.0e-70 : 1.0e-35f;
-  if (out_length > threshold) {
-    out_length = sqrt(out_length);
-    return v / out_length;
-  }
-  /* Either the vector is small or one of it's values contained `nan`. */
-  out_length = 0.0;
-  return T(0.0);
-}
-
-template<typename T, BLI_ENABLE_IF_FLT_VEC(T)> inline T normalize(const T &v)
-{
-  bT len;
-  return normalize_and_get_length(v, len);
-}
-
-template<typename T, BLI_ENABLE_IF_FLT_VEC(T), BLI_ENABLE_IF((T::type_length == 3))>
-inline T cross(const T &a, const T &b)
-{
-  return {a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x};
-}
-
-template<typename T,
-         BLI_ENABLE_IF((std::is_same_v<bT, float>)),
-         BLI_ENABLE_IF((T::type_length == 3))>
-inline T cross_high_precision(const T &a, const T &b)
-{
-  return {(float)((double)a.y * b.z - (double)a.z * b.y),
-          (float)((double)a.z * b.x - (double)a.x * b.z),
-          (float)((double)a.x * b.y - (double)a.y * b.x)};
-}
-
-template<typename T, BLI_ENABLE_IF_FLT_VEC(T), BLI_ENABLE_IF((T::type_length == 3))>
-inline T cross_poly(Span<T> poly)
-{
-  /* Newell's Method. */
-  int nv = static_cast<int>(poly.size());
-  if (nv < 3) {
-    return T(0, 0, 0);
-  }
-  const T *v_prev = &poly[nv - 1];
-  const T *v_curr = &poly[0];
-  T n(0, 0, 0);
-  for (int i = 0; i < nv;) {
-    n[0] = n[0] + ((*v_prev)[1] - (*v_curr)[1]) * ((*v_prev)[2] + (*v_curr)[2]);
-    n[1] = n[1] + ((*v_prev)[2] - (*v_curr)[2]) * ((*v_prev)[0] + (*v_curr)[0]);
-    n[2] = n[2] + ((*v_prev)[0] - (*v_curr)[0]) * ((*v_prev)[1] + (*v_curr)[1]);
-    v_prev = v_curr;
-    ++i;
-    if (i < nv) {
-      v_curr = &poly[i];
-    }
-  }
-  return n;
-}
-
-template<typename T, BLI_ENABLE_IF_FLT_VEC(T)> inline T interpolate(const T &a, const T &b, bT t)
-{
-  return a * (1 - t) + b * t;
-}
-
-template<typename T, BLI_ENABLE_IF_FLT_VEC(T)>
-inline T faceforward(const T &vector, const T &incident, const T &reference)
-{
-  return (dot(reference, incident) < 0) ? vector : -vector;
-}
-
-template<typename T> inline int dominant_axis(const T &a)
-{
-  T b = abs(a);
-  return ((b.x > b.y) ? ((b.x > b.z) ? 0 : 2) : ((b.y > b.z) ? 1 : 2));
-}
-
-/** Intersections. */
-
-template<typename T> struct isect_result {
-  enum {
-    LINE_LINE_COLINEAR = -1,
-    LINE_LINE_NONE = 0,
-    LINE_LINE_EXACT = 1,
-    LINE_LINE_CROSS = 2,
-  } kind;
-  bT lambda;
-};
-
-template<typename T, BLI_ENABLE_IF_FLT_VEC(T)>
-isect_result<T> isect_seg_seg(const T &v1, const T &v2, const T &v3, const T &v4);
-
-#undef BLI_ENABLE_IF_FLT_VEC
-#undef BLI_ENABLE_IF_INT_VEC
-#undef bT
-
-}  // namespace blender::math
diff --git a/source/blender/blenlib/BLI_memarena.h b/source/blender/blenlib/BLI_memarena.h
index 4ac4712bc8c..bcfe2efc5e5 100644
--- a/source/blender/blenlib/BLI_memarena.h
+++ b/source/blender/blenlib/BLI_memarena.h
@@ -38,13 +38,13 @@ extern "C" {
 struct MemArena;
 typedef struct MemArena MemArena;
 
-struct MemArena *BLI_memarena_new(const size_t bufsize,
+struct MemArena *BLI_memarena_new(size_t bufsize,
                                   const char *name) ATTR_WARN_UNUSED_RESULT ATTR_RETURNS_NONNULL
     ATTR_NONNULL(2) ATTR_MALLOC;
 void BLI_memarena_free(struct MemArena *ma) ATTR_NONNULL(1);
 void BLI_memarena_use_malloc(struct MemArena *ma) ATTR_NONNULL(1);
 void BLI_memarena_use_calloc(struct MemArena *ma) ATTR_NONNULL(1);
-void BLI_memarena_use_align(struct MemArena *ma, const size_t align) ATTR_NONNULL(1);
+void BLI_memarena_use_align(struct MemArena *ma, size_t align) ATTR_NONNULL(1);
 void *BLI_memarena_alloc(struct MemArena *ma, size_t size) ATTR_WARN_UNUSED_RESULT
     ATTR_NONNULL(1) ATTR_MALLOC ATTR_ALLOC_SIZE(2);
 void *BLI_memarena_calloc(struct MemArena *ma, size_t size) ATTR_WARN_UNUSED_RESULT
diff --git a/source/blender/blenlib/BLI_memory_utils.h b/source/blender/blenlib/BLI_memory_utils.h
index 79e25e26040..09d8f646905 100644
--- a/source/blender/blenlib/BLI_memory_utils.h
+++ b/source/blender/blenlib/BLI_memory_utils.h
@@ -29,7 +29,7 @@ extern "C" {
 
 /* it may be defined already */
 #ifndef __BLI_UTILDEFINES_H__
-bool BLI_memory_is_zero(const void *arr, const size_t size);
+bool BLI_memory_is_zero(const void *arr, size_t size);
 #endif
 
 #ifdef __cplusplus
diff --git a/source/blender/blenlib/BLI_memory_utils.hh b/source/blender/blenlib/BLI_memory_utils.hh
index 9a5be79b61e..37691017c12 100644
--- a/source/blender/blenlib/BLI_memory_utils.hh
+++ b/source/blender/blenlib/BLI_memory_utils.hh
@@ -557,4 +557,13 @@ Container &move_assign_container(Container &dst, Container &&src) noexcept(
   return dst;
 }
 
+/**
+ * Utility macro that wraps `std::enable_if` to make it a bit easier to use and less verbose for
+ * SFINAE in common cases.
+ *
+ * \note Often one has to invoke this macro with double parenthesis. That's because the condition
+ * often contains a comma and angle brackets are not recognized as parenthesis by the preprocessor.
+ */
+#define BLI_ENABLE_IF(condition) typename std::enable_if_t<condition> * = nullptr
+
 }  // namespace blender
diff --git a/source/blender/blenlib/BLI_mesh_intersect.hh b/source/blender/blenlib/BLI_mesh_intersect.hh
index 0ebee6f16a8..71a8abb822f 100644
--- a/source/blender/blenlib/BLI_mesh_intersect.hh
+++ b/source/blender/blenlib/BLI_mesh_intersect.hh
@@ -28,11 +28,12 @@
 #  include <iostream>
 
 #  include "BLI_array.hh"
+#  include "BLI_double3.hh"
+#  include "BLI_float3.hh"
 #  include "BLI_index_range.hh"
 #  include "BLI_map.hh"
 #  include "BLI_math_mpq.hh"
-#  include "BLI_math_vec_mpq_types.hh"
-#  include "BLI_math_vec_types.hh"
+#  include "BLI_mpq3.hh"
 #  include "BLI_span.hh"
 #  include "BLI_utility_mixins.hh"
 #  include "BLI_vector.hh"
diff --git a/source/blender/blenlib/BLI_mpq2.hh b/source/blender/blenlib/BLI_mpq2.hh
new file mode 100644
index 00000000000..18bc8821d9c
--- /dev/null
+++ b/source/blender/blenlib/BLI_mpq2.hh
@@ -0,0 +1,184 @@
+/*
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+#pragma once
+
+/** \file
+ * \ingroup bli
+ */
+
+#ifdef WITH_GMP
+
+#  include "BLI_math_mpq.hh"
+#  include "BLI_mpq3.hh"
+
+namespace blender {
+
+struct mpq2 {
+  mpq_class x, y;
+
+  mpq2() = default;
+
+  mpq2(const mpq_class *ptr) : x{ptr[0]}, y{ptr[1]}
+  {
+  }
+
+  mpq2(mpq_class x, mpq_class y) : x(x), y(y)
+  {
+  }
+
+  mpq2(const mpq2 &other) : x(other.x), y(other.y)
+  {
+  }
+
+  mpq2(mpq2 &&other) noexcept : x(std::move(other.x)), y(std::move(other.y))
+  {
+  }
+
+  ~mpq2() = default;
+
+  mpq2 &operator=(const mpq2 &other)
+  {
+    if (this != &other) {
+      x = other.x;
+      y = other.y;
+    }
+    return *this;
+  }
+
+  mpq2 &operator=(mpq2 &&other) noexcept
+  {
+    x = std::move(other.x);
+    y = std::move(other.y);
+    return *this;
+  }
+
+  mpq2(const mpq3 &other) : x(other.x), y(other.y)
+  {
+  }
+
+  operator mpq_class *()
+  {
+    return &x;
+  }
+
+  operator const mpq_class *() const
+  {
+    return &x;
+  }
+
+  /**
+   * Cannot do this exactly in rational arithmetic!
+   * Approximate by going in and out of doubles.
+   */
+  mpq_class length() const
+  {
+    mpq_class lsquared = dot(*this, *this);
+    return mpq_class(sqrt(lsquared.get_d()));
+  }
+
+  friend mpq2 operator+(const mpq2 &a, const mpq2 &b)
+  {
+    return {a.x + b.x, a.y + b.y};
+  }
+
+  friend mpq2 operator-(const mpq2 &a, const mpq2 &b)
+  {
+    return {a.x - b.x, a.y - b.y};
+  }
+
+  friend mpq2 operator*(const mpq2 &a, mpq_class b)
+  {
+    return {a.x * b, a.y * b};
+  }
+
+  friend mpq2 operator/(const mpq2 &a, mpq_class b)
+  {
+    BLI_assert(b != 0);
+    return {a.x / b, a.y / b};
+  }
+
+  friend mpq2 operator*(mpq_class a, const mpq2 &b)
+  {
+    return b * a;
+  }
+
+  friend bool operator==(const mpq2 &a, const mpq2 &b)
+  {
+    return a.x == b.x && a.y == b.y;
+  }
+
+  friend bool operator!=(const mpq2 &a, const mpq2 &b)
+  {
+    return a.x != b.x || a.y != b.y;
+  }
+
+  friend std::ostream &operator<<(std::ostream &stream, const mpq2 &v)
+  {
+    stream << "(" << v.x << ", " << v.y << ")";
+    return stream;
+  }
+
+  static mpq_class dot(const mpq2 &a, const mpq2 &b)
+  {
+    return a.x * b.x + a.y * b.y;
+  }
+
+  static mpq2 interpolate(const mpq2 &a, const mpq2 &b, mpq_class t)
+  {
+    return a * (1 - t) + b * t;
+  }
+
+  static mpq2 abs(const mpq2 &a)
+  {
+    mpq_class abs_x = (a.x >= 0) ? a.x : -a.x;
+    mpq_class abs_y = (a.y >= 0) ? a.y : -a.y;
+    return mpq2(abs_x, abs_y);
+  }
+
+  static mpq_class distance(const mpq2 &a, const mpq2 &b)
+  {
+    return (a - b).length();
+  }
+
+  static mpq_class distance_squared(const mpq2 &a, const mpq2 &b)
+  {
+    mpq2 diff = a - b;
+    return dot(diff, diff);
+  }
+
+  struct isect_result {
+    enum {
+      LINE_LINE_COLINEAR = -1,
+      LINE_LINE_NONE = 0,
+      LINE_LINE_EXACT = 1,
+      LINE_LINE_CROSS = 2,
+    } kind;
+    mpq_class lambda;
+  };
+
+  static isect_result isect_seg_seg(const mpq2 &v1,
+                                    const mpq2 &v2,
+                                    const mpq2 &v3,
+                                    const mpq2 &v4);
+
+  /** There is a sensible use for hashing on exact arithmetic types. */
+  uint64_t hash() const;
+};
+
+}  // namespace blender
+
+#endif /* WITH_GMP */
diff --git a/source/blender/blenlib/BLI_mpq3.hh b/source/blender/blenlib/BLI_mpq3.hh
new file mode 100644
index 00000000000..b9eda2ad7e1
--- /dev/null
+++ b/source/blender/blenlib/BLI_mpq3.hh
@@ -0,0 +1,297 @@
+/*
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+#pragma once
+
+/** \file
+ * \ingroup bli
+ */
+
+#ifdef WITH_GMP
+
+#  include <iostream>
+
+#  include "BLI_math.h"
+#  include "BLI_math_mpq.hh"
+#  include "BLI_span.hh"
+
+namespace blender {
+
+struct mpq3 {
+  mpq_class x, y, z;
+
+  mpq3() = default;
+
+  mpq3(const mpq_class *ptr) : x{ptr[0]}, y{ptr[1]}, z{ptr[2]}
+  {
+  }
+
+  mpq3(const mpq_class (*ptr)[3]) : mpq3((const mpq_class *)ptr)
+  {
+  }
+
+  explicit mpq3(mpq_class value) : x(value), y(value), z(value)
+  {
+  }
+
+  explicit mpq3(int value) : x(value), y(value), z(value)
+  {
+  }
+
+  mpq3(mpq_class x, mpq_class y, mpq_class z) : x{x}, y{y}, z{z}
+  {
+  }
+
+  operator const mpq_class *() const
+  {
+    return &x;
+  }
+
+  operator mpq_class *()
+  {
+    return &x;
+  }
+
+  /* Cannot do this exactly in rational arithmetic!
+   * Approximate by going in and out of doubles.
+   */
+  mpq_class normalize_and_get_length()
+  {
+    double dv[3] = {x.get_d(), y.get_d(), z.get_d()};
+    double len = normalize_v3_db(dv);
+    this->x = mpq_class(dv[0]);
+    this->y = mpq_class(dv[1]);
+    this->z = mpq_class(dv[2]);
+    return len;
+  }
+
+  mpq3 normalized() const
+  {
+    double dv[3] = {x.get_d(), y.get_d(), z.get_d()};
+    double dr[3];
+    normalize_v3_v3_db(dr, dv);
+    return mpq3(mpq_class(dr[0]), mpq_class(dr[1]), mpq_class(dr[2]));
+  }
+
+  /* Cannot do this exactly in rational arithmetic!
+   * Approximate by going in and out of double.
+   */
+  mpq_class length() const
+  {
+    mpq_class lsquared = this->length_squared();
+    double dsquared = lsquared.get_d();
+    double d = sqrt(dsquared);
+    return mpq_class(d);
+  }
+
+  mpq_class length_squared() const
+  {
+    return x * x + y * y + z * z;
+  }
+
+  void reflect(const mpq3 &normal)
+  {
+    *this = this->reflected(normal);
+  }
+
+  mpq3 reflected(const mpq3 &normal) const
+  {
+    mpq3 result;
+    const mpq_class dot2 = 2 * dot(*this, normal);
+    result.x = this->x - (dot2 * normal.x);
+    result.y = this->y - (dot2 * normal.y);
+    result.z = this->z - (dot2 * normal.z);
+    return result;
+  }
+
+  static mpq3 safe_divide(const mpq3 &a, const mpq3 &b)
+  {
+    mpq3 result;
+    result.x = (b.x == 0) ? mpq_class(0) : a.x / b.x;
+    result.y = (b.y == 0) ? mpq_class(0) : a.y / b.y;
+    result.z = (b.z == 0) ? mpq_class(0) : a.z / b.z;
+    return result;
+  }
+
+  void invert()
+  {
+    x = -x;
+    y = -y;
+    z = -z;
+  }
+
+  friend mpq3 operator+(const mpq3 &a, const mpq3 &b)
+  {
+    return mpq3(a.x + b.x, a.y + b.y, a.z + b.z);
+  }
+
+  void operator+=(const mpq3 &b)
+  {
+    this->x += b.x;
+    this->y += b.y;
+    this->z += b.z;
+  }
+
+  friend mpq3 operator-(const mpq3 &a, const mpq3 &b)
+  {
+    return mpq3(a.x - b.x, a.y - b.y, a.z - b.z);
+  }
+
+  friend mpq3 operator-(const mpq3 &a)
+  {
+    return mpq3(-a.x, -a.y, -a.z);
+  }
+
+  void operator-=(const mpq3 &b)
+  {
+    this->x -= b.x;
+    this->y -= b.y;
+    this->z -= b.z;
+  }
+
+  void operator*=(mpq_class scalar)
+  {
+    this->x *= scalar;
+    this->y *= scalar;
+    this->z *= scalar;
+  }
+
+  void operator*=(const mpq3 &other)
+  {
+    this->x *= other.x;
+    this->y *= other.y;
+    this->z *= other.z;
+  }
+
+  friend mpq3 operator*(const mpq3 &a, const mpq3 &b)
+  {
+    return {a.x * b.x, a.y * b.y, a.z * b.z};
+  }
+
+  friend mpq3 operator*(const mpq3 &a, const mpq_class &b)
+  {
+    return mpq3(a.x * b, a.y * b, a.z * b);
+  }
+
+  friend mpq3 operator*(const mpq_class &a, const mpq3 &b)
+  {
+    return mpq3(a * b.x, a * b.y, a * b.z);
+  }
+
+  friend mpq3 operator/(const mpq3 &a, const mpq_class &b)
+  {
+    BLI_assert(b != 0);
+    return mpq3(a.x / b, a.y / b, a.z / b);
+  }
+
+  friend bool operator==(const mpq3 &a, const mpq3 &b)
+  {
+    return a.x == b.x && a.y == b.y && a.z == b.z;
+  }
+
+  friend bool operator!=(const mpq3 &a, const mpq3 &b)
+  {
+    return a.x != b.x || a.y != b.y || a.z != b.z;
+  }
+
+  friend std::ostream &operator<<(std::ostream &stream, const mpq3 &v)
+  {
+    stream << "(" << v.x << ", " << v.y << ", " << v.z << ")";
+    return stream;
+  }
+
+  static mpq_class dot(const mpq3 &a, const mpq3 &b)
+  {
+    return a.x * b.x + a.y * b.y + a.z * b.z;
+  }
+
+  static mpq_class dot_with_buffer(const mpq3 &a, const mpq3 &b, mpq3 &buffer)
+  {
+    buffer = a;
+    buffer *= b;
+    buffer.x += buffer.y;
+    buffer.x += buffer.z;
+    return buffer.x;
+  }
+
+  static mpq3 cross(const mpq3 &a, const mpq3 &b)
+  {
+    return mpq3(a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0]);
+  }
+
+  static mpq3 cross_high_precision(const mpq3 &a, const mpq3 &b)
+  {
+    return cross(a, b);
+  }
+
+  static mpq3 project(const mpq3 &a, const mpq3 &b)
+  {
+    const mpq_class mul = mpq3::dot(a, b) / mpq3::dot(b, b);
+    return mpq3(mul * b[0], mul * b[1], mul * b[2]);
+  }
+
+  static mpq_class distance(const mpq3 &a, const mpq3 &b)
+  {
+    mpq3 diff(a.x - b.x, a.y - b.y, a.z - b.z);
+    return diff.length();
+  }
+
+  static mpq_class distance_squared(const mpq3 &a, const mpq3 &b)
+  {
+    mpq3 diff(a.x - b.x, a.y - b.y, a.z - b.z);
+    return mpq3::dot(diff, diff);
+  }
+
+  static mpq_class distance_squared_with_buffer(const mpq3 &a, const mpq3 &b, mpq3 &buffer)
+  {
+    buffer = a;
+    buffer -= b;
+    return mpq3::dot(buffer, buffer);
+  }
+
+  static mpq3 interpolate(const mpq3 &a, const mpq3 &b, mpq_class t)
+  {
+    mpq_class s = 1 - t;
+    return mpq3(a.x * s + b.x * t, a.y * s + b.y * t, a.z * s + b.z * t);
+  }
+
+  static mpq3 abs(const mpq3 &a)
+  {
+    mpq_class abs_x = (a.x >= 0) ? a.x : -a.x;
+    mpq_class abs_y = (a.y >= 0) ? a.y : -a.y;
+    mpq_class abs_z = (a.z >= 0) ? a.z : -a.z;
+    return mpq3(abs_x, abs_y, abs_z);
+  }
+
+  static int dominant_axis(const mpq3 &a)
+  {
+    mpq_class x = (a.x >= 0) ? a.x : -a.x;
+    mpq_class y = (a.y >= 0) ? a.y : -a.y;
+    mpq_class z = (a.z >= 0) ? a.z : -a.z;
+    return ((x > y) ? ((x > z) ? 0 : 2) : ((y > z) ? 1 : 2));
+  }
+
+  static mpq3 cross_poly(Span<mpq3> poly);
+
+  /** There is a sensible use for hashing on exact arithmetic types. */
+  uint64_t hash() const;
+};
+
+uint64_t hash_mpq_class(const mpq_class &value);
+
+}  // namespace blender
+
+#endif /* WITH_GMP */
diff --git a/source/blender/blenlib/BLI_noise.hh b/source/blender/blenlib/BLI_noise.hh
index 297c65c250a..4f68ef17ca2 100644
--- a/source/blender/blenlib/BLI_noise.hh
+++ b/source/blender/blenlib/BLI_noise.hh
@@ -16,7 +16,9 @@
 
 #pragma once
 
-#include "BLI_math_vec_types.hh"
+#include "BLI_float2.hh"
+#include "BLI_float3.hh"
+#include "BLI_float4.hh"
 
 namespace blender::noise {
 
diff --git a/source/blender/blenlib/BLI_path_util.h b/source/blender/blenlib/BLI_path_util.h
index 16f479cb3b8..658cc0c3825 100644
--- a/source/blender/blenlib/BLI_path_util.h
+++ b/source/blender/blenlib/BLI_path_util.h
@@ -75,16 +75,15 @@ bool BLI_make_existing_file(const char *name);
  * - Doesn't use CWD, or deal with relative paths.
  * - Only fill's in \a dir and \a file when they are non NULL.
  */
-void BLI_split_dirfile(
-    const char *string, char *dir, char *file, const size_t dirlen, const size_t filelen);
+void BLI_split_dirfile(const char *string, char *dir, char *file, size_t dirlen, size_t filelen);
 /**
  * Copies the parent directory part of string into `dir`, max length `dirlen`.
  */
-void BLI_split_dir_part(const char *string, char *dir, const size_t dirlen);
+void BLI_split_dir_part(const char *string, char *dir, size_t dirlen);
 /**
  * Copies the leaf filename part of string into `file`, max length `filelen`.
  */
-void BLI_split_file_part(const char *string, char *file, const size_t filelen);
+void BLI_split_file_part(const char *string, char *file, size_t filelen);
 /**
  * Returns a pointer to the last extension (e.g. the position of the last period).
  * Returns NULL if there is no extension.
@@ -94,7 +93,7 @@ const char *BLI_path_extension(const char *filepath) ATTR_NONNULL();
 /**
  * Append a filename to a dir, ensuring slash separates.
  */
-void BLI_path_append(char *__restrict dst, const size_t maxlen, const char *__restrict file)
+void BLI_path_append(char *__restrict dst, size_t maxlen, const char *__restrict file)
     ATTR_NONNULL();
 /**
  * Simple appending of filename to dir, does not check for valid path!
@@ -104,7 +103,7 @@ void BLI_path_append(char *__restrict dst, const size_t maxlen, const char *__re
  * that de-duplicates separators and can handle an arbitrary number of paths.
  */
 void BLI_join_dirfile(char *__restrict dst,
-                      const size_t maxlen,
+                      size_t maxlen,
                       const char *__restrict dir,
                       const char *__restrict file) ATTR_NONNULL();
 /**
@@ -114,7 +113,7 @@ void BLI_join_dirfile(char *__restrict dst,
  * \note If you want a trailing slash, add `SEP_STR` as the last path argument,
  * duplicate slashes will be cleaned up.
  */
-size_t BLI_path_join(char *__restrict dst, const size_t dst_len, const char *path_first, ...)
+size_t BLI_path_join(char *__restrict dst, size_t dst_len, const char *path_first, ...)
     ATTR_NONNULL(1, 3) ATTR_SENTINEL(0);
 /**
  * Like Python's `os.path.basename()`
@@ -164,12 +163,12 @@ void BLI_path_slash_rstrip(char *string) ATTR_NONNULL();
 void BLI_path_slash_native(char *path) ATTR_NONNULL();
 
 #ifdef _WIN32
-bool BLI_path_program_extensions_add_win32(char *name, const size_t maxlen);
+bool BLI_path_program_extensions_add_win32(char *name, size_t maxlen);
 #endif
 /**
  * Search for a binary (executable)
  */
-bool BLI_path_program_search(char *fullname, const size_t maxlen, const char *name);
+bool BLI_path_program_search(char *fullname, size_t maxlen, const char *name);
 
 /**
  * \return true when `str` end with `ext` (case insensitive).
@@ -353,7 +352,7 @@ bool BLI_path_is_abs_from_cwd(const char *path) ATTR_NONNULL();
  * This is _not_ something Blender's internal paths support, instead they use the "//" prefix.
  * In most cases #BLI_path_abs should be used instead.
  */
-bool BLI_path_abs_from_cwd(char *path, const size_t maxlen) ATTR_NONNULL();
+bool BLI_path_abs_from_cwd(char *path, size_t maxlen) ATTR_NONNULL();
 /**
  * Replaces `file` with a relative version (prefixed by "//") such that #BLI_path_abs, given
  * the same `relfile`, will convert it back to its original value.
diff --git a/source/blender/blenlib/BLI_rand.hh b/source/blender/blenlib/BLI_rand.hh
index 667d6df8996..cc9e9b374d7 100644
--- a/source/blender/blenlib/BLI_rand.hh
+++ b/source/blender/blenlib/BLI_rand.hh
@@ -20,8 +20,9 @@
 
 #pragma once
 
+#include "BLI_float2.hh"
+#include "BLI_float3.hh"
 #include "BLI_math.h"
-#include "BLI_math_vec_types.hh"
 #include "BLI_span.hh"
 #include "BLI_utildefines.h"
 
diff --git a/source/blender/blenlib/BLI_stack.h b/source/blender/blenlib/BLI_stack.h
index 653f5f61c9e..eb4e69a42d4 100644
--- a/source/blender/blenlib/BLI_stack.h
+++ b/source/blender/blenlib/BLI_stack.h
@@ -28,13 +28,13 @@ extern "C" {
 
 typedef struct BLI_Stack BLI_Stack;
 
-BLI_Stack *BLI_stack_new_ex(const size_t elem_size,
+BLI_Stack *BLI_stack_new_ex(size_t elem_size,
                             const char *description,
-                            const size_t chunk_size) ATTR_WARN_UNUSED_RESULT ATTR_NONNULL();
+                            size_t chunk_size) ATTR_WARN_UNUSED_RESULT ATTR_NONNULL();
 /**
  * Create a new homogeneous stack with elements of 'elem_size' bytes.
  */
-BLI_Stack *BLI_stack_new(const size_t elem_size, const char *description) ATTR_WARN_UNUSED_RESULT
+BLI_Stack *BLI_stack_new(size_t elem_size, const char *description) ATTR_WARN_UNUSED_RESULT
     ATTR_NONNULL();
 
 /**
diff --git a/source/blender/blenlib/BLI_string.h b/source/blender/blenlib/BLI_string.h
index a82e97914db..8177545911c 100644
--- a/source/blender/blenlib/BLI_string.h
+++ b/source/blender/blenlib/BLI_string.h
@@ -42,8 +42,7 @@ extern "C" {
  * \param len: The number of bytes to duplicate
  * \retval Returns the duplicated string
  */
-char *BLI_strdupn(const char *str, const size_t len) ATTR_MALLOC ATTR_WARN_UNUSED_RESULT
-    ATTR_NONNULL();
+char *BLI_strdupn(const char *str, size_t len) ATTR_MALLOC ATTR_WARN_UNUSED_RESULT ATTR_NONNULL();
 
 /**
  * Duplicates the cstring \a str into a newly mallocN'd
@@ -74,8 +73,7 @@ char *BLI_strdupcat(const char *__restrict str1,
  * the size of dst)
  * \retval Returns dst
  */
-char *BLI_strncpy(char *__restrict dst, const char *__restrict src, const size_t maxncpy)
-    ATTR_NONNULL();
+char *BLI_strncpy(char *__restrict dst, const char *__restrict src, size_t maxncpy) ATTR_NONNULL();
 
 /**
  * Like BLI_strncpy but ensures dst is always padded by given char,
@@ -107,7 +105,7 @@ char *BLI_strncpy_ensure_pad(char *__restrict dst,
  */
 size_t BLI_strncpy_rlen(char *__restrict dst,
                         const char *__restrict src,
-                        const size_t maxncpy) ATTR_WARN_UNUSED_RESULT ATTR_NONNULL();
+                        size_t maxncpy) ATTR_WARN_UNUSED_RESULT ATTR_NONNULL();
 
 size_t BLI_strcpy_rlen(char *__restrict dst, const char *__restrict src) ATTR_WARN_UNUSED_RESULT
     ATTR_NONNULL();
@@ -186,7 +184,7 @@ void BLI_str_replace_char(char *str, char src, char dst) ATTR_NONNULL();
  * \note Larger tables should use a hash table.
  */
 bool BLI_str_replace_table_exact(char *string,
-                                 const size_t string_len,
+                                 size_t string_len,
                                  const char *replace_table[][2],
                                  int replace_table_len);
 
@@ -235,7 +233,7 @@ char *BLI_sprintfN(const char *__restrict format, ...) ATTR_WARN_UNUSED_RESULT
  *
  * \note This is used for creating animation paths in blend files.
  */
-size_t BLI_str_escape(char *__restrict dst, const char *__restrict src, const size_t dst_maxncpy)
+size_t BLI_str_escape(char *__restrict dst, const char *__restrict src, size_t dst_maxncpy)
     ATTR_NONNULL();
 /**
  * This roughly matches C and Python's string escaping with double quotes - `"`.
@@ -251,9 +249,9 @@ size_t BLI_str_escape(char *__restrict dst, const char *__restrict src, const si
  */
 size_t BLI_str_unescape_ex(char *__restrict dst,
                            const char *__restrict src,
-                           const size_t src_maxncpy,
+                           size_t src_maxncpy,
                            /* Additional arguments. */
-                           const size_t dst_maxncpy,
+                           size_t dst_maxncpy,
                            bool *r_is_complete) ATTR_NONNULL();
 /**
  * See #BLI_str_unescape_ex doc-string.
@@ -265,7 +263,7 @@ size_t BLI_str_unescape_ex(char *__restrict dst,
  *
  * \note This is used for parsing animation paths in blend files (runs often).
  */
-size_t BLI_str_unescape(char *__restrict dst, const char *__restrict src, const size_t src_maxncpy)
+size_t BLI_str_unescape(char *__restrict dst, const char *__restrict src, size_t src_maxncpy)
     ATTR_NONNULL();
 
 /**
@@ -359,10 +357,10 @@ int BLI_strcmp_ignore_pad(const char *str1, const char *str2, char pad) ATTR_WAR
 /**
  * Determine the length of a fixed-size string.
  */
-size_t BLI_strnlen(const char *str, const size_t maxlen) ATTR_WARN_UNUSED_RESULT ATTR_NONNULL();
+size_t BLI_strnlen(const char *str, size_t maxlen) ATTR_WARN_UNUSED_RESULT ATTR_NONNULL();
 
-void BLI_str_tolower_ascii(char *str, const size_t len) ATTR_NONNULL();
-void BLI_str_toupper_ascii(char *str, const size_t len) ATTR_NONNULL();
+void BLI_str_tolower_ascii(char *str, size_t len) ATTR_NONNULL();
+void BLI_str_toupper_ascii(char *str, size_t len) ATTR_NONNULL();
 /**
  * Strip white-space from end of the string.
  */
@@ -479,7 +477,7 @@ bool BLI_string_all_words_matched(const char *name,
  * \return The number of words found in \a str
  */
 int BLI_string_find_split_words(const char *str,
-                                const size_t len,
+                                size_t len,
                                 char delim,
                                 int r_words[][2],
                                 int words_max) ATTR_WARN_UNUSED_RESULT ATTR_NONNULL();
diff --git a/source/blender/blenlib/BLI_string_utf8.h b/source/blender/blenlib/BLI_string_utf8.h
index 82622d442fb..108a2f5fc7d 100644
--- a/source/blender/blenlib/BLI_string_utf8.h
+++ b/source/blender/blenlib/BLI_string_utf8.h
@@ -110,14 +110,12 @@ size_t BLI_str_utf8_from_unicode_len(unsigned int c) ATTR_WARN_UNUSED_RESULT;
  *
  * \return number of bytes written.
  */
-size_t BLI_str_utf8_from_unicode(unsigned int c, char *outbuf, const size_t outbuf_len)
-    ATTR_NONNULL(2);
+size_t BLI_str_utf8_from_unicode(unsigned int c, char *outbuf, size_t outbuf_len) ATTR_NONNULL(2);
 size_t BLI_str_utf8_as_utf32(char32_t *__restrict dst_w,
                              const char *__restrict src_c,
-                             const size_t maxncpy) ATTR_NONNULL(1, 2);
-size_t BLI_str_utf32_as_utf8(char *__restrict dst,
-                             const char32_t *__restrict src,
-                             const size_t maxncpy) ATTR_NONNULL(1, 2);
+                             size_t maxncpy) ATTR_NONNULL(1, 2);
+size_t BLI_str_utf32_as_utf8(char *__restrict dst, const char32_t *__restrict src, size_t maxncpy)
+    ATTR_NONNULL(1, 2);
 /**
  * \return The UTF-32 len in UTF-8.
  */
@@ -162,21 +160,20 @@ size_t BLI_wstrlen_utf8(const wchar_t *src) ATTR_NONNULL(1) ATTR_WARN_UNUSED_RES
 size_t BLI_strlen_utf8_ex(const char *strc, size_t *r_len_bytes)
     ATTR_NONNULL(1, 2) ATTR_WARN_UNUSED_RESULT;
 size_t BLI_strlen_utf8(const char *strc) ATTR_NONNULL(1) ATTR_WARN_UNUSED_RESULT;
-size_t BLI_strnlen_utf8_ex(const char *strc, const size_t maxlen, size_t *r_len_bytes)
+size_t BLI_strnlen_utf8_ex(const char *strc, size_t maxlen, size_t *r_len_bytes)
     ATTR_NONNULL(1, 3);
 /**
  * \param strc: the string to measure the length.
  * \param maxlen: the string length (in bytes)
  * \return the unicode length (not in bytes!)
  */
-size_t BLI_strnlen_utf8(const char *strc, const size_t maxlen)
-    ATTR_NONNULL(1) ATTR_WARN_UNUSED_RESULT;
+size_t BLI_strnlen_utf8(const char *strc, size_t maxlen) ATTR_NONNULL(1) ATTR_WARN_UNUSED_RESULT;
 size_t BLI_strncpy_wchar_as_utf8(char *__restrict dst,
                                  const wchar_t *__restrict src,
-                                 const size_t maxncpy) ATTR_NONNULL(1, 2);
+                                 size_t maxncpy) ATTR_NONNULL(1, 2);
 size_t BLI_strncpy_wchar_from_utf8(wchar_t *__restrict dst,
                                    const char *__restrict src,
-                                   const size_t maxncpy) ATTR_NONNULL(1, 2);
+                                   size_t maxncpy) ATTR_NONNULL(1, 2);
 
 /**
  * Count columns that character/string occupies (based on `wcwidth.co`).
diff --git a/source/blender/blenlib/BLI_string_utils.h b/source/blender/blenlib/BLI_string_utils.h
index fd3918ff217..818bfe8182b 100644
--- a/source/blender/blenlib/BLI_string_utils.h
+++ b/source/blender/blenlib/BLI_string_utils.h
@@ -57,11 +57,11 @@ bool BLI_string_is_decimal(const char *string) ATTR_NONNULL();
  * Based on `BLI_split_dirfile()` / `os.path.splitext()`,
  * `"a.b.c"` -> (`"a.b"`, `".c"`).
  */
-void BLI_string_split_suffix(const char *string, char *r_body, char *r_suf, const size_t str_len);
+void BLI_string_split_suffix(const char *string, char *r_body, char *r_suf, size_t str_len);
 /**
  * `"a.b.c"` -> (`"a."`, `"b.c"`).
  */
-void BLI_string_split_prefix(const char *string, char *r_pre, char *r_body, const size_t str_len);
+void BLI_string_split_prefix(const char *string, char *r_pre, char *r_body, size_t str_len);
 
 /**
  * Join strings, return newly allocated string.
@@ -127,7 +127,7 @@ char *BLI_string_join_array_by_sep_char_with_tableN(char sep,
 size_t BLI_string_flip_side_name(char *r_name,
                                  const char *from_name,
                                  bool strip_number,
-                                 const size_t name_len);
+                                 size_t name_len);
 
 /**
  * Ensures name is unique (according to criteria specified by caller in unique_check callback),
diff --git a/source/blender/blenlib/BLI_timecode.h b/source/blender/blenlib/BLI_timecode.h
index f0349e289ac..1cd18dc86ab 100644
--- a/source/blender/blenlib/BLI_timecode.h
+++ b/source/blender/blenlib/BLI_timecode.h
@@ -42,7 +42,7 @@ extern "C" {
  * \return length of \a str
  */
 size_t BLI_timecode_string_from_time(char *str,
-                                     const size_t maxncpy,
+                                     size_t maxncpy,
                                      int brevity_level,
                                      float time_seconds,
                                      double fps,
@@ -56,7 +56,7 @@ size_t BLI_timecode_string_from_time(char *str,
  * \param time_seconds: time total time in seconds
  * \return length of \a str
  */
-size_t BLI_timecode_string_from_time_simple(char *str, const size_t maxncpy, double time_seconds)
+size_t BLI_timecode_string_from_time_simple(char *str, size_t maxncpy, double time_seconds)
     ATTR_NONNULL();
 
 /**
@@ -72,7 +72,7 @@ size_t BLI_timecode_string_from_time_simple(char *str, const size_t maxncpy, dou
  * \note in some cases this is used to print non-seconds values.
  */
 size_t BLI_timecode_string_from_time_seconds(char *str,
-                                             const size_t maxncpy,
+                                             size_t maxncpy,
                                              int brevity_level,
                                              float time_seconds) ATTR_NONNULL();
 
diff --git a/source/blender/blenlib/BLI_utildefines.h b/source/blender/blenlib/BLI_utildefines.h
index 9fe092fe525..35d4158de59 100644
--- a/source/blender/blenlib/BLI_utildefines.h
+++ b/source/blender/blenlib/BLI_utildefines.h
@@ -638,7 +638,7 @@ extern "C" {
 /**
  * Check if memory is zeroed, as with `memset(arr, 0, arr_size)`.
  */
-extern bool BLI_memory_is_zero(const void *arr, const size_t arr_size);
+extern bool BLI_memory_is_zero(const void *arr, size_t arr_size);
 #endif
 
 #define MEMCMP_STRUCT_AFTER_IS_ZERO(struct_var, member) \
@@ -840,15 +840,6 @@ extern bool BLI_memory_is_zero(const void *arr, const size_t arr_size);
 /** No-op for expressions we don't want to instantiate, but must remain valid. */
 #define EXPR_NOP(expr) (void)(0 ? ((void)(expr), 1) : 0)
 
-/**
- * Utility macro that wraps `std::enable_if` to make it a bit easier to use and less verbose for
- * SFINAE in common cases.
- *
- * \note Often one has to invoke this macro with double parenthesis. That's because the condition
- * often contains a comma and angle brackets are not recognized as parenthesis by the preprocessor.
- */
-#define BLI_ENABLE_IF(condition) typename std::enable_if_t<(condition)> * = nullptr
-
 /** \} */
 
 #ifdef __cplusplus
diff --git a/source/blender/blenlib/CMakeLists.txt b/source/blender/blenlib/CMakeLists.txt
index 90c6760019a..3958fd8e2d2 100644
--- a/source/blender/blenlib/CMakeLists.txt
+++ b/source/blender/blenlib/CMakeLists.txt
@@ -192,6 +192,8 @@ set(SRC
   BLI_dlrbTree.h
   BLI_dot_export.hh
   BLI_dot_export_attribute_enums.hh
+  BLI_double2.hh
+  BLI_double3.hh
   BLI_dynlib.h
   BLI_dynstr.h
   BLI_easing.h
@@ -205,6 +207,9 @@ set(SRC
   BLI_fileops.hh
   BLI_fileops_types.h
   BLI_filereader.h
+  BLI_float2.hh
+  BLI_float3.hh
+  BLI_float4.hh
   BLI_float4x4.hh
   BLI_fnmatch.h
   BLI_function_ref.hh
@@ -253,8 +258,6 @@ set(SRC
   BLI_math_statistics.h
   BLI_math_time.h
   BLI_math_vector.h
-  BLI_math_vec_types.hh
-  BLI_math_vec_mpq_types.hh
   BLI_memarena.h
   BLI_memblock.h
   BLI_memiter.h
@@ -264,6 +267,8 @@ set(SRC
   BLI_mesh_boolean.hh
   BLI_mesh_intersect.hh
   BLI_mmap.h
+  BLI_mpq2.hh
+  BLI_mpq3.hh
   BLI_multi_value_map.hh
   BLI_noise.h
   BLI_noise.hh
@@ -439,7 +444,6 @@ if(WITH_GTESTS)
     tests/BLI_math_rotation_test.cc
     tests/BLI_math_solvers_test.cc
     tests/BLI_math_time_test.cc
-    tests/BLI_math_vec_types_test.cc
     tests/BLI_math_vector_test.cc
     tests/BLI_memiter_test.cc
     tests/BLI_memory_utils_test.cc
diff --git a/source/blender/blenlib/intern/BLI_mempool_private.h b/source/blender/blenlib/intern/BLI_mempool_private.h
index 90569d87c41..03b0b11297b 100644
--- a/source/blender/blenlib/intern/BLI_mempool_private.h
+++ b/source/blender/blenlib/intern/BLI_mempool_private.h
@@ -54,8 +54,9 @@ typedef struct ParallelMempoolTaskData {
  *
  * See #BLI_task_parallel_mempool implementation for detailed usage example.
  */
-ParallelMempoolTaskData *mempool_iter_threadsafe_create(BLI_mempool *pool, const size_t num_iter)
-    ATTR_WARN_UNUSED_RESULT ATTR_NONNULL();
+ParallelMempoolTaskData *mempool_iter_threadsafe_create(BLI_mempool *pool,
+                                                        size_t num_iter) ATTR_WARN_UNUSED_RESULT
+    ATTR_NONNULL();
 void mempool_iter_threadsafe_destroy(ParallelMempoolTaskData *iter_arr) ATTR_NONNULL();
 
 /**
diff --git a/source/blender/blenlib/intern/delaunay_2d.cc b/source/blender/blenlib/intern/delaunay_2d.cc
index 842e6cb6135..53e881a9fc7 100644
--- a/source/blender/blenlib/intern/delaunay_2d.cc
+++ b/source/blender/blenlib/intern/delaunay_2d.cc
@@ -25,10 +25,11 @@
 #include <sstream>
 
 #include "BLI_array.hh"
+#include "BLI_double2.hh"
 #include "BLI_linklist.h"
 #include "BLI_math_boolean.hh"
 #include "BLI_math_mpq.hh"
-#include "BLI_math_vec_mpq_types.hh"
+#include "BLI_mpq2.hh"
 #include "BLI_set.hh"
 #include "BLI_task.hh"
 #include "BLI_vector.hh"
@@ -37,8 +38,6 @@
 
 namespace blender::meshintersect {
 
-using namespace blender::math;
-
 /* Throughout this file, template argument T will be an
  * arithmetic-like type, like float, double, or mpq_class. */
 
@@ -789,11 +788,11 @@ bool in_line<mpq_class>(const FatCo<mpq_class> &a,
   }
   vec2<mpq_class> exact_ab = b.exact - a.exact;
   vec2<mpq_class> exact_ac = c.exact - a.exact;
-  if (dot(exact_ab, exact_ac) < 0) {
+  if (vec2<mpq_class>::dot(exact_ab, exact_ac) < 0) {
     return false;
   }
   vec2<mpq_class> exact_bc = c.exact - b.exact;
-  return dot(exact_bc, exact_ac) >= 0;
+  return vec2<mpq_class>::dot(exact_bc, exact_ac) >= 0;
 }
 #endif
 
@@ -802,11 +801,11 @@ bool in_line<double>(const FatCo<double> &a, const FatCo<double> &b, const FatCo
 {
   vec2<double> ab = b.approx - a.approx;
   vec2<double> ac = c.approx - a.approx;
-  if (dot(ab, ac) < 0) {
+  if (vec2<double>::dot(ab, ac) < 0) {
     return false;
   }
   vec2<double> bc = c.approx - b.approx;
-  return dot(bc, ac) >= 0;
+  return vec2<double>::dot(bc, ac) >= 0;
 }
 
 template<> CDTVert<double>::CDTVert(const vec2<double> &pt)
@@ -1082,7 +1081,7 @@ template<typename T> CDTEdge<T> *CDTArrangement<T>::split_edge(SymEdge<T> *se, T
   SymEdge<T> *sesymprev = prev(sesym);
   SymEdge<T> *sesymprevsym = sym(sesymprev);
   SymEdge<T> *senext = se->next;
-  CDTVert<T> *v = this->add_vert(interpolate(*a, *b, lambda));
+  CDTVert<T> *v = this->add_vert(vec2<T>::interpolate(*a, *b, lambda));
   CDTEdge<T> *e = this->add_edge(v, se->next->vert, se->face, sesym->face);
   sesym->vert = v;
   SymEdge<T> *newse = &e->symedges[0];
@@ -1705,16 +1704,16 @@ void fill_crossdata_for_intersect(const FatCo<T> &curco,
   BLI_assert(se_vcva->vert == vc && se_vcva->next->vert == va);
   BLI_assert(se_vcvb->vert == vc && se_vcvb->next->vert == vb);
   UNUSED_VARS_NDEBUG(vc);
-  auto isect = isect_seg_seg<vec2<T>>(va->co.exact, vb->co.exact, curco.exact, v2->co.exact);
+  auto isect = vec2<T>::isect_seg_seg(va->co.exact, vb->co.exact, curco.exact, v2->co.exact);
   T &lambda = isect.lambda;
   switch (isect.kind) {
-    case isect_result<vec2<T>>::LINE_LINE_CROSS: {
+    case vec2<T>::isect_result::LINE_LINE_CROSS: {
 #ifdef WITH_GMP
       if (!std::is_same<T, mpq_class>::value) {
 #else
       if (true) {
 #endif
-        double len_ab = distance(va->co.approx, vb->co.approx);
+        double len_ab = vec2<double>::distance(va->co.approx, vb->co.approx);
         if (lambda * len_ab <= epsilon) {
           fill_crossdata_for_through_vert(va, se_vcva, cd, cd_next);
         }
@@ -1736,7 +1735,7 @@ void fill_crossdata_for_intersect(const FatCo<T> &curco,
       }
       break;
     }
-    case isect_result<vec2<T>>::LINE_LINE_EXACT: {
+    case vec2<T>::isect_result::LINE_LINE_EXACT: {
       if (lambda == 0) {
         fill_crossdata_for_through_vert(va, se_vcva, cd, cd_next);
       }
@@ -1751,7 +1750,7 @@ void fill_crossdata_for_intersect(const FatCo<T> &curco,
       }
       break;
     }
-    case isect_result<vec2<T>>::LINE_LINE_NONE: {
+    case vec2<T>::isect_result::LINE_LINE_NONE: {
 #ifdef WITH_GMP
       if (std::is_same<T, mpq_class>::value) {
         BLI_assert(false);
@@ -1767,9 +1766,9 @@ void fill_crossdata_for_intersect(const FatCo<T> &curco,
       }
       break;
     }
-    case isect_result<vec2<T>>::LINE_LINE_COLINEAR: {
-      if (distance_squared(va->co.approx, v2->co.approx) <=
-          distance_squared(vb->co.approx, v2->co.approx)) {
+    case vec2<T>::isect_result::LINE_LINE_COLINEAR: {
+      if (vec2<double>::distance_squared(va->co.approx, v2->co.approx) <=
+          vec2<double>::distance_squared(vb->co.approx, v2->co.approx)) {
         fill_crossdata_for_through_vert(va, se_vcva, cd, cd_next);
       }
       else {
@@ -1846,7 +1845,7 @@ void get_next_crossing_from_edge(CrossData<T> *cd,
 {
   CDTVert<T> *va = cd->in->vert;
   CDTVert<T> *vb = cd->in->next->vert;
-  vec2<T> curco = interpolate(va->co.exact, vb->co.exact, cd->lambda);
+  vec2<T> curco = vec2<T>::interpolate(va->co.exact, vb->co.exact, cd->lambda);
   FatCo<T> fat_curco(curco);
   SymEdge<T> *se_ac = sym(cd->in)->next;
   CDTVert<T> *vc = se_ac->next->vert;
@@ -2387,7 +2386,7 @@ template<typename T> void remove_non_constraint_edges_leave_valid_bmesh(CDT_stat
       dissolvable_edges[i].e = e;
       const vec2<double> &co1 = e->symedges[0].vert->co.approx;
       const vec2<double> &co2 = e->symedges[1].vert->co.approx;
-      dissolvable_edges[i].len_squared = distance_squared(co1, co2);
+      dissolvable_edges[i].len_squared = vec2<double>::distance_squared(co1, co2);
       i++;
     }
   }
@@ -2570,18 +2569,18 @@ template<typename T> void detect_holes(CDT_state<T> *cdt_state)
           if (e->symedges[0].face->visit_index == e->symedges[1].face->visit_index) {
             continue; /* Don't count hits on edges between faces in same region. */
           }
-          auto isect = isect_seg_seg<vec2<T>>(ray_end.exact,
+          auto isect = vec2<T>::isect_seg_seg(ray_end.exact,
                                               mid.exact,
                                               e->symedges[0].vert->co.exact,
                                               e->symedges[1].vert->co.exact);
           switch (isect.kind) {
-            case isect_result<vec2<T>>::LINE_LINE_CROSS: {
+            case vec2<T>::isect_result::LINE_LINE_CROSS: {
               hits++;
               break;
             }
-            case isect_result<vec2<T>>::LINE_LINE_EXACT:
-            case isect_result<vec2<T>>::LINE_LINE_NONE:
-            case isect_result<vec2<T>>::LINE_LINE_COLINEAR:
+            case vec2<T>::isect_result::LINE_LINE_EXACT:
+            case vec2<T>::isect_result::LINE_LINE_NONE:
+            case vec2<T>::isect_result::LINE_LINE_COLINEAR:
               break;
           }
         }
diff --git a/source/blender/blenlib/intern/math_boolean.cc b/source/blender/blenlib/intern/math_boolean.cc
index 0bae3c23f79..c16755868aa 100644
--- a/source/blender/blenlib/intern/math_boolean.cc
+++ b/source/blender/blenlib/intern/math_boolean.cc
@@ -18,10 +18,15 @@
  * \ingroup bli
  */
 
+#include "BLI_double2.hh"
+#include "BLI_double3.hh"
+#include "BLI_float2.hh"
+#include "BLI_float3.hh"
 #include "BLI_hash.hh"
 #include "BLI_math_boolean.hh"
 #include "BLI_math_mpq.hh"
-#include "BLI_math_vec_types.hh"
+#include "BLI_mpq2.hh"
+#include "BLI_mpq3.hh"
 #include "BLI_span.hh"
 #include "BLI_utildefines.h"
 
diff --git a/source/blender/blenlib/intern/math_vec.cc b/source/blender/blenlib/intern/math_vec.cc
index 6fab6c9a383..223c0e273f0 100644
--- a/source/blender/blenlib/intern/math_vec.cc
+++ b/source/blender/blenlib/intern/math_vec.cc
@@ -18,83 +18,89 @@
  * \ingroup bli
  */
 
+#include "BLI_double2.hh"
+#include "BLI_double3.hh"
+#include "BLI_float2.hh"
+#include "BLI_float3.hh"
 #include "BLI_hash.hh"
-#include "BLI_math_vec_mpq_types.hh"
-#include "BLI_math_vector.hh"
+#include "BLI_math_mpq.hh"
+#include "BLI_mpq2.hh"
+#include "BLI_mpq3.hh"
 #include "BLI_span.hh"
 #include "BLI_utildefines.h"
 
-namespace blender::math {
+namespace blender {
 
-template<>
-isect_result<float2> isect_seg_seg(const float2 &v1,
-                                   const float2 &v2,
-                                   const float2 &v3,
-                                   const float2 &v4)
+float2::isect_result float2::isect_seg_seg(const float2 &v1,
+                                           const float2 &v2,
+                                           const float2 &v3,
+                                           const float2 &v4)
 {
-  isect_result<float2> ans;
+  float2::isect_result ans;
   float div = (v2[0] - v1[0]) * (v4[1] - v3[1]) - (v2[1] - v1[1]) * (v4[0] - v3[0]);
   if (div == 0.0f) {
     ans.lambda = 0.0f;
-    ans.kind = isect_result<float2>::LINE_LINE_COLINEAR;
+    ans.mu = 0.0f;
+    ans.kind = float2::isect_result::LINE_LINE_COLINEAR;
   }
   else {
     ans.lambda = ((v1[1] - v3[1]) * (v4[0] - v3[0]) - (v1[0] - v3[0]) * (v4[1] - v3[1])) / div;
-    float mu = ((v1[1] - v3[1]) * (v2[0] - v1[0]) - (v1[0] - v3[0]) * (v2[1] - v1[1])) / div;
-    if (ans.lambda >= 0.0f && ans.lambda <= 1.0f && mu >= 0.0f && mu <= 1.0f) {
-      if (ans.lambda == 0.0f || ans.lambda == 1.0f || mu == 0.0f || mu == 1.0f) {
-        ans.kind = isect_result<float2>::LINE_LINE_EXACT;
+    ans.mu = ((v1[1] - v3[1]) * (v2[0] - v1[0]) - (v1[0] - v3[0]) * (v2[1] - v1[1])) / div;
+    if (ans.lambda >= 0.0f && ans.lambda <= 1.0f && ans.mu >= 0.0f && ans.mu <= 1.0f) {
+      if (ans.lambda == 0.0f || ans.lambda == 1.0f || ans.mu == 0.0f || ans.mu == 1.0f) {
+        ans.kind = float2::isect_result::LINE_LINE_EXACT;
       }
       else {
-        ans.kind = isect_result<float2>::LINE_LINE_CROSS;
+        ans.kind = float2::isect_result::LINE_LINE_CROSS;
       }
     }
     else {
-      ans.kind = isect_result<float2>::LINE_LINE_NONE;
+      ans.kind = float2::isect_result::LINE_LINE_NONE;
     }
   }
   return ans;
 }
 
-template<>
-isect_result<double2> isect_seg_seg(const double2 &v1,
-                                    const double2 &v2,
-                                    const double2 &v3,
-                                    const double2 &v4)
+double2::isect_result double2::isect_seg_seg(const double2 &v1,
+                                             const double2 &v2,
+                                             const double2 &v3,
+                                             const double2 &v4)
 {
-  isect_result<double2> ans;
+  double2::isect_result ans;
   double div = (v2[0] - v1[0]) * (v4[1] - v3[1]) - (v2[1] - v1[1]) * (v4[0] - v3[0]);
   if (div == 0.0) {
     ans.lambda = 0.0;
-    ans.kind = isect_result<double2>::LINE_LINE_COLINEAR;
+    ans.kind = double2::isect_result::LINE_LINE_COLINEAR;
   }
   else {
     ans.lambda = ((v1[1] - v3[1]) * (v4[0] - v3[0]) - (v1[0] - v3[0]) * (v4[1] - v3[1])) / div;
     double mu = ((v1[1] - v3[1]) * (v2[0] - v1[0]) - (v1[0] - v3[0]) * (v2[1] - v1[1])) / div;
     if (ans.lambda >= 0.0 && ans.lambda <= 1.0 && mu >= 0.0 && mu <= 1.0) {
       if (ans.lambda == 0.0 || ans.lambda == 1.0 || mu == 0.0 || mu == 1.0) {
-        ans.kind = isect_result<double2>::LINE_LINE_EXACT;
+        ans.kind = double2::isect_result::LINE_LINE_EXACT;
       }
       else {
-        ans.kind = isect_result<double2>::LINE_LINE_CROSS;
+        ans.kind = double2::isect_result::LINE_LINE_CROSS;
       }
     }
     else {
-      ans.kind = isect_result<double2>::LINE_LINE_NONE;
+      ans.kind = double2::isect_result::LINE_LINE_NONE;
     }
   }
   return ans;
 }
 
 #ifdef WITH_GMP
-template<>
-isect_result<mpq2> isect_seg_seg(const mpq2 &v1, const mpq2 &v2, const mpq2 &v3, const mpq2 &v4)
+mpq2::isect_result mpq2::isect_seg_seg(const mpq2 &v1,
+                                       const mpq2 &v2,
+                                       const mpq2 &v3,
+                                       const mpq2 &v4)
 {
-  isect_result<mpq2> ans;
+  mpq2::isect_result ans;
   mpq_class div = (v2[0] - v1[0]) * (v4[1] - v3[1]) - (v2[1] - v1[1]) * (v4[0] - v3[0]);
   if (div == 0.0) {
     ans.lambda = 0.0;
-    ans.kind = isect_result<mpq2>::LINE_LINE_COLINEAR;
+    ans.kind = mpq2::isect_result::LINE_LINE_COLINEAR;
   }
   else {
     ans.lambda = ((v1[1] - v3[1]) * (v4[0] - v3[0]) - (v1[0] - v3[0]) * (v4[1] - v3[1])) / div;
@@ -103,21 +109,66 @@ isect_result<mpq2> isect_seg_seg(const mpq2 &v1, const mpq2 &v2, const mpq2 &v3,
     if (ans.lambda >= 0 && ans.lambda <= 1 &&
         ((div > 0 && mudiv >= 0 && mudiv <= div) || (div < 0 && mudiv <= 0 && mudiv >= div))) {
       if (ans.lambda == 0 || ans.lambda == 1 || mudiv == 0 || mudiv == div) {
-        ans.kind = isect_result<mpq2>::LINE_LINE_EXACT;
+        ans.kind = mpq2::isect_result::LINE_LINE_EXACT;
       }
       else {
-        ans.kind = isect_result<mpq2>::LINE_LINE_CROSS;
+        ans.kind = mpq2::isect_result::LINE_LINE_CROSS;
       }
     }
     else {
-      ans.kind = isect_result<mpq2>::LINE_LINE_NONE;
+      ans.kind = mpq2::isect_result::LINE_LINE_NONE;
     }
   }
   return ans;
 }
 #endif
 
+double3 double3::cross_poly(Span<double3> poly)
+{
+  /* Newell's Method. */
+  int nv = static_cast<int>(poly.size());
+  if (nv < 3) {
+    return double3(0, 0, 0);
+  }
+  const double3 *v_prev = &poly[nv - 1];
+  const double3 *v_curr = &poly[0];
+  double3 n(0, 0, 0);
+  for (int i = 0; i < nv;) {
+    n[0] = n[0] + ((*v_prev)[1] - (*v_curr)[1]) * ((*v_prev)[2] + (*v_curr)[2]);
+    n[1] = n[1] + ((*v_prev)[2] - (*v_curr)[2]) * ((*v_prev)[0] + (*v_curr)[0]);
+    n[2] = n[2] + ((*v_prev)[0] - (*v_curr)[0]) * ((*v_prev)[1] + (*v_curr)[1]);
+    v_prev = v_curr;
+    ++i;
+    if (i < nv) {
+      v_curr = &poly[i];
+    }
+  }
+  return n;
+}
+
 #ifdef WITH_GMP
+mpq3 mpq3::cross_poly(Span<mpq3> poly)
+{
+  /* Newell's Method. */
+  int nv = static_cast<int>(poly.size());
+  if (nv < 3) {
+    return mpq3(0);
+  }
+  const mpq3 *v_prev = &poly[nv - 1];
+  const mpq3 *v_curr = &poly[0];
+  mpq3 n(0);
+  for (int i = 0; i < nv;) {
+    n[0] = n[0] + ((*v_prev)[1] - (*v_curr)[1]) * ((*v_prev)[2] + (*v_curr)[2]);
+    n[1] = n[1] + ((*v_prev)[2] - (*v_curr)[2]) * ((*v_prev)[0] + (*v_curr)[0]);
+    n[2] = n[2] + ((*v_prev)[0] - (*v_curr)[0]) * ((*v_prev)[1] + (*v_curr)[1]);
+    v_prev = v_curr;
+    ++i;
+    if (i < nv) {
+      v_curr = &poly[i];
+    }
+  }
+  return n;
+}
 
 uint64_t hash_mpq_class(const mpq_class &value)
 {
@@ -125,6 +176,20 @@ uint64_t hash_mpq_class(const mpq_class &value)
   return get_default_hash(static_cast<float>(value.get_d()));
 }
 
+uint64_t mpq2::hash() const
+{
+  uint64_t hashx = hash_mpq_class(this->x);
+  uint64_t hashy = hash_mpq_class(this->y);
+  return hashx ^ (hashy * 33);
+}
+
+uint64_t mpq3::hash() const
+{
+  uint64_t hashx = hash_mpq_class(this->x);
+  uint64_t hashy = hash_mpq_class(this->y);
+  uint64_t hashz = hash_mpq_class(this->z);
+  return hashx ^ (hashy * 33) ^ (hashz * 33 * 37);
+}
 #endif
 
-}  // namespace blender::math
+}  // namespace blender
diff --git a/source/blender/blenlib/intern/mesh_boolean.cc b/source/blender/blenlib/intern/mesh_boolean.cc
index a3eae1896d3..ce4db0c6b9d 100644
--- a/source/blender/blenlib/intern/mesh_boolean.cc
+++ b/source/blender/blenlib/intern/mesh_boolean.cc
@@ -28,6 +28,8 @@
 #  include "BLI_array.hh"
 #  include "BLI_assert.h"
 #  include "BLI_delaunay_2d.h"
+#  include "BLI_double3.hh"
+#  include "BLI_float3.hh"
 #  include "BLI_hash.hh"
 #  include "BLI_kdopbvh.h"
 #  include "BLI_map.hh"
@@ -35,9 +37,8 @@
 #  include "BLI_math_boolean.hh"
 #  include "BLI_math_geom.h"
 #  include "BLI_math_mpq.hh"
-#  include "BLI_math_vec_mpq_types.hh"
-#  include "BLI_math_vec_types.hh"
 #  include "BLI_mesh_intersect.hh"
+#  include "BLI_mpq3.hh"
 #  include "BLI_set.hh"
 #  include "BLI_span.hh"
 #  include "BLI_stack.hh"
@@ -1632,13 +1633,13 @@ static Edge find_good_sorting_edge(const Vert *testp,
   ordinate[axis_next] = -abscissa[axis];
   ordinate[axis_next_next] = 0;
   /* By construction, dot(abscissa, ordinate) == 0, so they are perpendicular. */
-  mpq3 normal = math::cross(abscissa, ordinate);
+  mpq3 normal = mpq3::cross(abscissa, ordinate);
   if (dbg_level > 0) {
     std::cout << "abscissa = " << abscissa << "\n";
     std::cout << "ordinate = " << ordinate << "\n";
     std::cout << "normal = " << normal << "\n";
   }
-  mpq_class nlen2 = math::length_squared(normal);
+  mpq_class nlen2 = normal.length_squared();
   mpq_class max_abs_slope = -1;
   Edge esort;
   const Vector<Edge> &edges = tmtopo.vert_edges(closestp);
@@ -1647,12 +1648,12 @@ static Edge find_good_sorting_edge(const Vert *testp,
     const mpq3 &co_other = v_other->co_exact;
     mpq3 evec = co_other - co_closest;
     /* Get projection of evec onto plane of abscissa and ordinate. */
-    mpq3 proj_evec = evec - (math::dot(evec, normal) / nlen2) * normal;
+    mpq3 proj_evec = evec - (mpq3::dot(evec, normal) / nlen2) * normal;
     /* The projection calculations along the abscissa and ordinate should
      * be scaled by 1/abscissa and 1/ordinate respectively,
      * but we can skip: it won't affect which `evec` has the maximum slope. */
-    mpq_class evec_a = math::dot(proj_evec, abscissa);
-    mpq_class evec_o = math::dot(proj_evec, ordinate);
+    mpq_class evec_a = mpq3::dot(proj_evec, abscissa);
+    mpq_class evec_o = mpq3::dot(proj_evec, ordinate);
     if (dbg_level > 0) {
       std::cout << "e = " << e << "\n";
       std::cout << "v_other = " << v_other << "\n";
@@ -1790,8 +1791,8 @@ static mpq_class closest_on_tri_to_point(const mpq3 &p,
   ap = p;
   ap -= a;
 
-  mpq_class d1 = math::dot_with_buffer(ab, ap, m);
-  mpq_class d2 = math::dot_with_buffer(ac, ap, m);
+  mpq_class d1 = mpq3::dot_with_buffer(ab, ap, m);
+  mpq_class d2 = mpq3::dot_with_buffer(ac, ap, m);
   if (d1 <= 0 && d2 <= 0) {
     /* Barycentric coordinates (1,0,0). */
     *r_edge = -1;
@@ -1799,13 +1800,13 @@ static mpq_class closest_on_tri_to_point(const mpq3 &p,
     if (dbg_level > 0) {
       std::cout << "  answer = a\n";
     }
-    return math::distance_squared_with_buffer(p, a, m);
+    return mpq3::distance_squared_with_buffer(p, a, m);
   }
   /* Check if p in vertex region outside b. */
   bp = p;
   bp -= b;
-  mpq_class d3 = math::dot_with_buffer(ab, bp, m);
-  mpq_class d4 = math::dot_with_buffer(ac, bp, m);
+  mpq_class d3 = mpq3::dot_with_buffer(ab, bp, m);
+  mpq_class d4 = mpq3::dot_with_buffer(ac, bp, m);
   if (d3 >= 0 && d4 <= d3) {
     /* Barycentric coordinates (0,1,0). */
     *r_edge = -1;
@@ -1813,7 +1814,7 @@ static mpq_class closest_on_tri_to_point(const mpq3 &p,
     if (dbg_level > 0) {
       std::cout << "  answer = b\n";
     }
-    return math::distance_squared_with_buffer(p, b, m);
+    return mpq3::distance_squared_with_buffer(p, b, m);
   }
   /* Check if p in region of ab. */
   mpq_class vc = d1 * d4 - d3 * d2;
@@ -1828,13 +1829,13 @@ static mpq_class closest_on_tri_to_point(const mpq3 &p,
     if (dbg_level > 0) {
       std::cout << "  answer = on ab at " << r << "\n";
     }
-    return math::distance_squared_with_buffer(p, r, m);
+    return mpq3::distance_squared_with_buffer(p, r, m);
   }
   /* Check if p in vertex region outside c. */
   cp = p;
   cp -= c;
-  mpq_class d5 = math::dot_with_buffer(ab, cp, m);
-  mpq_class d6 = math::dot_with_buffer(ac, cp, m);
+  mpq_class d5 = mpq3::dot_with_buffer(ab, cp, m);
+  mpq_class d6 = mpq3::dot_with_buffer(ac, cp, m);
   if (d6 >= 0 && d5 <= d6) {
     /* Barycentric coordinates (0,0,1). */
     *r_edge = -1;
@@ -1842,7 +1843,7 @@ static mpq_class closest_on_tri_to_point(const mpq3 &p,
     if (dbg_level > 0) {
       std::cout << "  answer = c\n";
     }
-    return math::distance_squared_with_buffer(p, c, m);
+    return mpq3::distance_squared_with_buffer(p, c, m);
   }
   /* Check if p in edge region of ac. */
   mpq_class vb = d5 * d2 - d1 * d6;
@@ -1857,7 +1858,7 @@ static mpq_class closest_on_tri_to_point(const mpq3 &p,
     if (dbg_level > 0) {
       std::cout << "  answer = on ac at " << r << "\n";
     }
-    return math::distance_squared_with_buffer(p, r, m);
+    return mpq3::distance_squared_with_buffer(p, r, m);
   }
   /* Check if p in edge region of bc. */
   mpq_class va = d3 * d6 - d5 * d4;
@@ -1873,7 +1874,7 @@ static mpq_class closest_on_tri_to_point(const mpq3 &p,
     if (dbg_level > 0) {
       std::cout << "  answer = on bc at " << r << "\n";
     }
-    return math::distance_squared_with_buffer(p, r, m);
+    return mpq3::distance_squared_with_buffer(p, r, m);
   }
   /* p inside face region. Compute barycentric coordinates (u,v,w). */
   mpq_class denom = 1 / (va + vb + vc);
@@ -1889,7 +1890,7 @@ static mpq_class closest_on_tri_to_point(const mpq3 &p,
   if (dbg_level > 0) {
     std::cout << "  answer = inside at " << r << "\n";
   }
-  return math::distance_squared_with_buffer(p, r, m);
+  return mpq3::distance_squared_with_buffer(p, r, m);
 }
 
 static float closest_on_tri_to_point_float_dist_squared(const float3 &p,
@@ -2609,7 +2610,7 @@ static void test_tri_inside_shapes(const IMesh &tm,
   double3 test_point = calc_point_inside_tri_db(tri_test);
   /* Offset the test point a tiny bit in the tri_test normal direction. */
   tri_test.populate_plane(false);
-  double3 norm = math::normalize(tri_test.plane->norm);
+  double3 norm = tri_test.plane->norm.normalized();
   const double offset_amount = 1e-5;
   double3 offset_test_point = test_point + offset_amount * norm;
   if (dbg_level > 0) {
@@ -3001,7 +3002,7 @@ static void init_face_merge_state(FaceMergeState *fms,
       std::cout << "process tri = " << &tri << "\n";
     }
     BLI_assert(tri.plane_populated());
-    if (math::dot(norm, tri.plane->norm) <= 0.0) {
+    if (double3::dot(norm, tri.plane->norm) <= 0.0) {
       if (dbg_level > 0) {
         std::cout << "triangle has wrong orientation, skipping\n";
       }
@@ -3026,7 +3027,7 @@ static void init_face_merge_state(FaceMergeState *fms,
       }
       if (me_index == -1) {
         double3 vec = new_me.v2->co - new_me.v1->co;
-        new_me.len_squared = math::length_squared(vec);
+        new_me.len_squared = vec.length_squared();
         new_me.orig = tri.edge_orig[i];
         new_me.is_intersect = tri.is_intersect[i];
         new_me.dissolvable = (new_me.orig == NO_INDEX && !new_me.is_intersect);
@@ -3266,7 +3267,7 @@ static Vector<Face *> merge_tris_for_face(Vector<int> tris,
   bool done = false;
   double3 first_tri_normal = tm.face(tris[0])->plane->norm;
   double3 second_tri_normal = tm.face(tris[1])->plane->norm;
-  if (tris.size() == 2 && math::dot(first_tri_normal, second_tri_normal) > 0.0) {
+  if (tris.size() == 2 && double3::dot(first_tri_normal, second_tri_normal) > 0.0) {
     /* Is this a case where quad with one diagonal remained unchanged?
      * Worth special handling because this case will be very common. */
     Face &tri1 = *tm.face(tris[0]);
@@ -3331,7 +3332,7 @@ static bool approx_in_line(const double3 &a, const double3 &b, const double3 &c)
 {
   double3 vec1 = b - a;
   double3 vec2 = c - b;
-  double cos_ang = math::dot(math::normalize(vec1), math::normalize(vec2));
+  double cos_ang = double3::dot(vec1.normalized(), vec2.normalized());
   return fabs(cos_ang - 1.0) < 1e-4;
 }
 
diff --git a/source/blender/blenlib/intern/mesh_intersect.cc b/source/blender/blenlib/intern/mesh_intersect.cc
index 982759ffcff..1f150137ba3 100644
--- a/source/blender/blenlib/intern/mesh_intersect.cc
+++ b/source/blender/blenlib/intern/mesh_intersect.cc
@@ -30,13 +30,15 @@
 #  include "BLI_array.hh"
 #  include "BLI_assert.h"
 #  include "BLI_delaunay_2d.h"
+#  include "BLI_double3.hh"
+#  include "BLI_float3.hh"
 #  include "BLI_hash.hh"
 #  include "BLI_kdopbvh.h"
 #  include "BLI_map.hh"
 #  include "BLI_math_boolean.hh"
 #  include "BLI_math_mpq.hh"
-#  include "BLI_math_vec_mpq_types.hh"
-#  include "BLI_math_vec_types.hh"
+#  include "BLI_mpq2.hh"
+#  include "BLI_mpq3.hh"
 #  include "BLI_polyfill_2d.h"
 #  include "BLI_set.hh"
 #  include "BLI_span.hh"
@@ -196,14 +198,14 @@ void Face::populate_plane(bool need_exact)
       for (int i : index_range()) {
         co[i] = vert[i]->co_exact;
       }
-      normal_exact = math::cross_poly(co.as_span());
+      normal_exact = mpq3::cross_poly(co);
     }
     else {
       mpq3 tr02 = vert[0]->co_exact - vert[2]->co_exact;
       mpq3 tr12 = vert[1]->co_exact - vert[2]->co_exact;
-      normal_exact = math::cross(tr02, tr12);
+      normal_exact = mpq3::cross(tr02, tr12);
     }
-    mpq_class d_exact = -math::dot(normal_exact, vert[0]->co_exact);
+    mpq_class d_exact = -mpq3::dot(normal_exact, vert[0]->co_exact);
     plane = new Plane(normal_exact, d_exact);
   }
   else {
@@ -213,14 +215,14 @@ void Face::populate_plane(bool need_exact)
       for (int i : index_range()) {
         co[i] = vert[i]->co;
       }
-      normal = math::cross_poly(co.as_span());
+      normal = double3::cross_poly(co);
     }
     else {
       double3 tr02 = vert[0]->co - vert[2]->co;
       double3 tr12 = vert[1]->co - vert[2]->co;
-      normal = math::cross(tr02, tr12);
+      normal = double3::cross_high_precision(tr02, tr12);
     }
-    double d = -math::dot(normal, vert[0]->co);
+    double d = -double3::dot(normal, vert[0]->co);
     plane = new Plane(normal, d);
   }
 }
@@ -1096,15 +1098,15 @@ static mpq2 project_3d_to_2d(const mpq3 &p3d, int proj_axis)
  */
 static double supremum_dot_cross(const double3 &a, const double3 &b)
 {
-  double3 abs_a = math::abs(a);
-  double3 abs_b = math::abs(b);
+  double3 abs_a = double3::abs(a);
+  double3 abs_b = double3::abs(b);
   double3 c;
   /* This is dot(cross(a, b), cross(a,b)) but using absolute values for a and b
    * and always using + when operation is + or -. */
   c[0] = abs_a[1] * abs_b[2] + abs_a[2] * abs_b[1];
   c[1] = abs_a[2] * abs_b[0] + abs_a[0] * abs_b[2];
   c[2] = abs_a[0] * abs_b[1] + abs_a[1] * abs_b[0];
-  return math::dot(c, c);
+  return double3::dot(c, c);
 }
 
 /* The index of dot when inputs are plane_coords with index 1 is much higher.
@@ -1141,11 +1143,11 @@ static int filter_plane_side(const double3 &p,
                              const double3 &abs_plane_p,
                              const double3 &abs_plane_no)
 {
-  double d = math::dot(p - plane_p, plane_no);
+  double d = double3::dot(p - plane_p, plane_no);
   if (d == 0.0) {
     return 0;
   }
-  double supremum = math::dot(abs_p + abs_plane_p, abs_plane_no);
+  double supremum = double3::dot(abs_p + abs_plane_p, abs_plane_no);
   double err_bound = supremum * index_plane_side * DBL_EPSILON;
   if (fabs(d) > err_bound) {
     return d > 0 ? 1 : -1;
@@ -1176,9 +1178,9 @@ static inline mpq3 tti_interp(
   ab -= b;
   ac = a;
   ac -= c;
-  mpq_class den = math::dot_with_buffer(ab, n, dotbuf);
+  mpq_class den = mpq3::dot_with_buffer(ab, n, dotbuf);
   BLI_assert(den != 0);
-  mpq_class alpha = math::dot_with_buffer(ac, n, dotbuf) / den;
+  mpq_class alpha = mpq3::dot_with_buffer(ac, n, dotbuf) / den;
   return a - alpha * ab;
 }
 
@@ -1207,7 +1209,7 @@ static inline int tti_above(const mpq3 &a,
   n.y = ba.z * ca.x - ba.x * ca.z;
   n.z = ba.x * ca.y - ba.y * ca.x;
 
-  return sgn(math::dot_with_buffer(ad, n, dotbuf));
+  return sgn(mpq3::dot_with_buffer(ad, n, dotbuf));
 }
 
 /**
@@ -1426,11 +1428,11 @@ static ITT_value intersect_tri_tri(const IMesh &tm, int t1, int t2)
   const double3 &d_r2 = vr2->co;
   const double3 &d_n2 = tri2.plane->norm;
 
-  const double3 &abs_d_p1 = math::abs(d_p1);
-  const double3 &abs_d_q1 = math::abs(d_q1);
-  const double3 &abs_d_r1 = math::abs(d_r1);
-  const double3 &abs_d_r2 = math::abs(d_r2);
-  const double3 &abs_d_n2 = math::abs(d_n2);
+  const double3 &abs_d_p1 = double3::abs(d_p1);
+  const double3 &abs_d_q1 = double3::abs(d_q1);
+  const double3 &abs_d_r1 = double3::abs(d_r1);
+  const double3 &abs_d_r2 = double3::abs(d_r2);
+  const double3 &abs_d_n2 = double3::abs(d_n2);
 
   int sp1 = filter_plane_side(d_p1, d_r2, d_n2, abs_d_p1, abs_d_r2, abs_d_n2);
   int sq1 = filter_plane_side(d_q1, d_r2, d_n2, abs_d_q1, abs_d_r2, abs_d_n2);
@@ -1446,9 +1448,9 @@ static ITT_value intersect_tri_tri(const IMesh &tm, int t1, int t2)
   }
 
   const double3 &d_n1 = tri1.plane->norm;
-  const double3 &abs_d_p2 = math::abs(d_p2);
-  const double3 &abs_d_q2 = math::abs(d_q2);
-  const double3 &abs_d_n1 = math::abs(d_n1);
+  const double3 &abs_d_p2 = double3::abs(d_p2);
+  const double3 &abs_d_q2 = double3::abs(d_q2);
+  const double3 &abs_d_n1 = double3::abs(d_n1);
 
   int sp2 = filter_plane_side(d_p2, d_r1, d_n1, abs_d_p2, abs_d_r1, abs_d_n1);
   int sq2 = filter_plane_side(d_q2, d_r1, d_n1, abs_d_q2, abs_d_r1, abs_d_n1);
@@ -1475,17 +1477,17 @@ static ITT_value intersect_tri_tri(const IMesh &tm, int t1, int t2)
   if (sp1 == 0) {
     buf[0] = p1;
     buf[0] -= r2;
-    sp1 = sgn(math::dot_with_buffer(buf[0], n2, buf[1]));
+    sp1 = sgn(mpq3::dot_with_buffer(buf[0], n2, buf[1]));
   }
   if (sq1 == 0) {
     buf[0] = q1;
     buf[0] -= r2;
-    sq1 = sgn(math::dot_with_buffer(buf[0], n2, buf[1]));
+    sq1 = sgn(mpq3::dot_with_buffer(buf[0], n2, buf[1]));
   }
   if (sr1 == 0) {
     buf[0] = r1;
     buf[0] -= r2;
-    sr1 = sgn(math::dot_with_buffer(buf[0], n2, buf[1]));
+    sr1 = sgn(mpq3::dot_with_buffer(buf[0], n2, buf[1]));
   }
 
   if (dbg_level > 1) {
@@ -1507,17 +1509,17 @@ static ITT_value intersect_tri_tri(const IMesh &tm, int t1, int t2)
   if (sp2 == 0) {
     buf[0] = p2;
     buf[0] -= r1;
-    sp2 = sgn(math::dot_with_buffer(buf[0], n1, buf[1]));
+    sp2 = sgn(mpq3::dot_with_buffer(buf[0], n1, buf[1]));
   }
   if (sq2 == 0) {
     buf[0] = q2;
     buf[0] -= r1;
-    sq2 = sgn(math::dot_with_buffer(buf[0], n1, buf[1]));
+    sq2 = sgn(mpq3::dot_with_buffer(buf[0], n1, buf[1]));
   }
   if (sr2 == 0) {
     buf[0] = r2;
     buf[0] -= r1;
-    sr2 = sgn(math::dot_with_buffer(buf[0], n1, buf[1]));
+    sr2 = sgn(mpq3::dot_with_buffer(buf[0], n1, buf[1]));
   }
 
   if (dbg_level > 1) {
@@ -1719,7 +1721,7 @@ static CDT_data prepare_cdt_input(const IMesh &tm, int t, const Vector<ITT_value
   BLI_assert(tm.face(t)->plane_populated());
   ans.t_plane = tm.face(t)->plane;
   BLI_assert(ans.t_plane->exact_populated());
-  ans.proj_axis = math::dominant_axis(ans.t_plane->norm_exact);
+  ans.proj_axis = mpq3::dominant_axis(ans.t_plane->norm_exact);
   prepare_need_tri(ans, tm, t);
   for (const ITT_value &itt : itts) {
     switch (itt.kind) {
@@ -1755,7 +1757,7 @@ static CDT_data prepare_cdt_input_for_cluster(const IMesh &tm,
   BLI_assert(tm.face(t0)->plane_populated());
   ans.t_plane = tm.face(t0)->plane;
   BLI_assert(ans.t_plane->exact_populated());
-  ans.proj_axis = math::dominant_axis(ans.t_plane->norm_exact);
+  ans.proj_axis = mpq3::dominant_axis(ans.t_plane->norm_exact);
   for (const int t : cl) {
     prepare_need_tri(ans, tm, t);
   }
@@ -2002,9 +2004,9 @@ static bool is_quad_flip_first_third(const double3 &v1,
                                      const double3 &normal)
 {
   double3 dir_v3v1 = v3 - v1;
-  double3 tangent = math::cross(dir_v3v1, normal);
-  double dot = math::dot(v1, tangent);
-  return (math::dot(v4, tangent) >= dot) || (math::dot(v2, tangent) <= dot);
+  double3 tangent = double3::cross_high_precision(dir_v3v1, normal);
+  double dot = double3::dot(v1, tangent);
+  return (double3::dot(v4, tangent) >= dot) || (double3::dot(v2, tangent) <= dot);
 }
 
 /**
@@ -2122,7 +2124,7 @@ static Array<Face *> exact_triangulate_poly(Face *f, IMeshArena *arena)
     f->populate_plane(false);
   }
   const double3 &poly_normal = f->plane->norm;
-  int axis = math::dominant_axis(poly_normal);
+  int axis = double3::dominant_axis(poly_normal);
   /* If project down y axis as opposed to x or z, the orientation
    * of the polygon will be reversed.
    * Yet another reversal happens if the poly normal in the dominant
@@ -2201,15 +2203,15 @@ static bool face_is_degenerate(const Face *f)
   }
   double3 da = v2->co - v0->co;
   double3 db = v2->co - v1->co;
-  double3 dab = math::cross(da, db);
-  double dab_length_squared = math::length_squared(dab);
+  double3 dab = double3::cross_high_precision(da, db);
+  double dab_length_squared = dab.length_squared();
   double err_bound = supremum_dot_cross(dab, dab) * index_dot_cross * DBL_EPSILON;
   if (dab_length_squared > err_bound) {
     return false;
   }
   mpq3 a = v2->co_exact - v0->co_exact;
   mpq3 b = v2->co_exact - v1->co_exact;
-  mpq3 ab = math::cross(a, b);
+  mpq3 ab = mpq3::cross(a, b);
   if (ab.x == 0 && ab.y == 0 && ab.z == 0) {
     return true;
   }
@@ -2229,8 +2231,8 @@ static bool any_degenerate_tris_fast(const Array<Face *> triangulation)
     }
     double3 da = v2->co - v0->co;
     double3 db = v2->co - v1->co;
-    double da_length_squared = math::length_squared(da);
-    double db_length_squared = math::length_squared(db);
+    double da_length_squared = da.length_squared();
+    double db_length_squared = db.length_squared();
     if (da_length_squared == 0.0 || db_length_squared == 0.0) {
       return true;
     }
@@ -2238,8 +2240,8 @@ static bool any_degenerate_tris_fast(const Array<Face *> triangulation)
      * The triangle is almost degenerate if sin t is almost 0.
      * sin^2 t = |da x db|^2 / (|da|^2 |db|^2)
      */
-    double3 dab = math::cross(da, db);
-    double dab_length_squared = math::length_squared(dab);
+    double3 dab = double3::cross_high_precision(da, db);
+    double dab_length_squared = dab.length_squared();
     double sin_squared_t = dab_length_squared / (da_length_squared * db_length_squared);
     if (sin_squared_t < 1e-8) {
       return true;
diff --git a/source/blender/blenlib/intern/noise.cc b/source/blender/blenlib/intern/noise.cc
index 3460c1284fc..a6ad18801fd 100644
--- a/source/blender/blenlib/intern/noise.cc
+++ b/source/blender/blenlib/intern/noise.cc
@@ -50,7 +50,9 @@
 #include <cmath>
 #include <cstdint>
 
-#include "BLI_math_vec_types.hh"
+#include "BLI_float2.hh"
+#include "BLI_float3.hh"
+#include "BLI_float4.hh"
 #include "BLI_math_base_safe.h"
 #include "BLI_noise.hh"
 #include "BLI_utildefines.h"
@@ -1467,7 +1469,7 @@ void voronoi_smooth_f1(const float w,
       correctionFactor /= 1.0f + 3.0f * smoothness;
       if (r_color != nullptr) {
         const float3 cellColor = hash_float_to_float3(cellPosition + cellOffset);
-        smoothColor = math::interpolate(smoothColor, cellColor, h) - correctionFactor;
+        smoothColor = float3::interpolate(smoothColor, cellColor, h) - correctionFactor;
       }
       if (r_w != nullptr) {
         smoothPosition = mix(smoothPosition, pointPosition, h) - correctionFactor;
@@ -1590,7 +1592,7 @@ static float voronoi_distance(const float2 a,
 {
   switch (metric) {
     case NOISE_SHD_VORONOI_EUCLIDEAN:
-      return math::distance(a, b);
+      return float2::distance(a, b);
     case NOISE_SHD_VORONOI_MANHATTAN:
       return fabsf(a.x - b.x) + fabsf(a.y - b.y);
     case NOISE_SHD_VORONOI_CHEBYCHEV:
@@ -1613,7 +1615,7 @@ void voronoi_f1(const float2 coord,
                 float3 *r_color,
                 float2 *r_position)
 {
-  const float2 cellPosition = math::floor(coord);
+  const float2 cellPosition = float2::floor(coord);
   const float2 localPosition = coord - cellPosition;
 
   float minDistance = 8.0f;
@@ -1652,7 +1654,7 @@ void voronoi_smooth_f1(const float2 coord,
                        float3 *r_color,
                        float2 *r_position)
 {
-  const float2 cellPosition = math::floor(coord);
+  const float2 cellPosition = float2::floor(coord);
   const float2 localPosition = coord - cellPosition;
   const float smoothness_clamped = max_ff(smoothness, FLT_MIN);
 
@@ -1674,10 +1676,11 @@ void voronoi_smooth_f1(const float2 coord,
         correctionFactor /= 1.0f + 3.0f * smoothness;
         if (r_color != nullptr) {
           const float3 cellColor = hash_float_to_float3(cellPosition + cellOffset);
-          smoothColor = math::interpolate(smoothColor, cellColor, h) - correctionFactor;
+          smoothColor = float3::interpolate(smoothColor, cellColor, h) - correctionFactor;
         }
         if (r_position != nullptr) {
-          smoothPosition = math::interpolate(smoothPosition, pointPosition, h) - correctionFactor;
+          smoothPosition = float2::interpolate(smoothPosition, pointPosition, h) -
+                           correctionFactor;
         }
       }
     }
@@ -1701,7 +1704,7 @@ void voronoi_f2(const float2 coord,
                 float3 *r_color,
                 float2 *r_position)
 {
-  const float2 cellPosition = math::floor(coord);
+  const float2 cellPosition = float2::floor(coord);
   const float2 localPosition = coord - cellPosition;
 
   float distanceF1 = 8.0f;
@@ -1745,7 +1748,7 @@ void voronoi_f2(const float2 coord,
 
 void voronoi_distance_to_edge(const float2 coord, const float randomness, float *r_distance)
 {
-  const float2 cellPosition = math::floor(coord);
+  const float2 cellPosition = float2::floor(coord);
   const float2 localPosition = coord - cellPosition;
 
   float2 vectorToClosest = float2(0.0f, 0.0f);
@@ -1774,7 +1777,7 @@ void voronoi_distance_to_edge(const float2 coord, const float randomness, float
       const float2 perpendicularToEdge = vectorToPoint - vectorToClosest;
       if (dot_v2v2(perpendicularToEdge, perpendicularToEdge) > 0.0001f) {
         const float distanceToEdge = dot_v2v2((vectorToClosest + vectorToPoint) / 2.0f,
-                                              math::normalize(perpendicularToEdge));
+                                              perpendicularToEdge.normalized());
         minDistance = std::min(minDistance, distanceToEdge);
       }
     }
@@ -1784,7 +1787,7 @@ void voronoi_distance_to_edge(const float2 coord, const float randomness, float
 
 void voronoi_n_sphere_radius(const float2 coord, const float randomness, float *r_radius)
 {
-  const float2 cellPosition = math::floor(coord);
+  const float2 cellPosition = float2::floor(coord);
   const float2 localPosition = coord - cellPosition;
 
   float2 closestPoint = float2(0.0f, 0.0f);
@@ -1795,7 +1798,7 @@ void voronoi_n_sphere_radius(const float2 coord, const float randomness, float *
       const float2 cellOffset = float2(i, j);
       const float2 pointPosition = cellOffset +
                                    hash_float_to_float2(cellPosition + cellOffset) * randomness;
-      const float distanceToPoint = math::distance(pointPosition, localPosition);
+      const float distanceToPoint = float2::distance(pointPosition, localPosition);
       if (distanceToPoint < minDistance) {
         minDistance = distanceToPoint;
         closestPoint = pointPosition;
@@ -1814,14 +1817,14 @@ void voronoi_n_sphere_radius(const float2 coord, const float randomness, float *
       const float2 cellOffset = float2(i, j) + closestPointOffset;
       const float2 pointPosition = cellOffset +
                                    hash_float_to_float2(cellPosition + cellOffset) * randomness;
-      const float distanceToPoint = math::distance(closestPoint, pointPosition);
+      const float distanceToPoint = float2::distance(closestPoint, pointPosition);
       if (distanceToPoint < minDistance) {
         minDistance = distanceToPoint;
         closestPointToClosestPoint = pointPosition;
       }
     }
   }
-  *r_radius = math::distance(closestPointToClosestPoint, closestPoint) / 2.0f;
+  *r_radius = float2::distance(closestPointToClosestPoint, closestPoint) / 2.0f;
 }
 
 /* **** 3D Voronoi **** */
@@ -1833,7 +1836,7 @@ static float voronoi_distance(const float3 a,
 {
   switch (metric) {
     case NOISE_SHD_VORONOI_EUCLIDEAN:
-      return math::distance(a, b);
+      return float3::distance(a, b);
     case NOISE_SHD_VORONOI_MANHATTAN:
       return fabsf(a.x - b.x) + fabsf(a.y - b.y) + fabsf(a.z - b.z);
     case NOISE_SHD_VORONOI_CHEBYCHEV:
@@ -1857,7 +1860,7 @@ void voronoi_f1(const float3 coord,
                 float3 *r_color,
                 float3 *r_position)
 {
-  const float3 cellPosition = math::floor(coord);
+  const float3 cellPosition = float3::floor(coord);
   const float3 localPosition = coord - cellPosition;
 
   float minDistance = 8.0f;
@@ -1899,7 +1902,7 @@ void voronoi_smooth_f1(const float3 coord,
                        float3 *r_color,
                        float3 *r_position)
 {
-  const float3 cellPosition = math::floor(coord);
+  const float3 cellPosition = float3::floor(coord);
   const float3 localPosition = coord - cellPosition;
   const float smoothness_clamped = max_ff(smoothness, FLT_MIN);
 
@@ -1922,10 +1925,10 @@ void voronoi_smooth_f1(const float3 coord,
           correctionFactor /= 1.0f + 3.0f * smoothness;
           if (r_color != nullptr) {
             const float3 cellColor = hash_float_to_float3(cellPosition + cellOffset);
-            smoothColor = math::interpolate(smoothColor, cellColor, h) - correctionFactor;
+            smoothColor = float3::interpolate(smoothColor, cellColor, h) - correctionFactor;
           }
           if (r_position != nullptr) {
-            smoothPosition = math::interpolate(smoothPosition, pointPosition, h) -
+            smoothPosition = float3::interpolate(smoothPosition, pointPosition, h) -
                              correctionFactor;
           }
         }
@@ -1951,7 +1954,7 @@ void voronoi_f2(const float3 coord,
                 float3 *r_color,
                 float3 *r_position)
 {
-  const float3 cellPosition = math::floor(coord);
+  const float3 cellPosition = float3::floor(coord);
   const float3 localPosition = coord - cellPosition;
 
   float distanceF1 = 8.0f;
@@ -1997,7 +2000,7 @@ void voronoi_f2(const float3 coord,
 
 void voronoi_distance_to_edge(const float3 coord, const float randomness, float *r_distance)
 {
-  const float3 cellPosition = math::floor(coord);
+  const float3 cellPosition = float3::floor(coord);
   const float3 localPosition = coord - cellPosition;
 
   float3 vectorToClosest = float3(0.0f, 0.0f, 0.0f);
@@ -2029,7 +2032,7 @@ void voronoi_distance_to_edge(const float3 coord, const float randomness, float
         const float3 perpendicularToEdge = vectorToPoint - vectorToClosest;
         if (dot_v3v3(perpendicularToEdge, perpendicularToEdge) > 0.0001f) {
           const float distanceToEdge = dot_v3v3((vectorToClosest + vectorToPoint) / 2.0f,
-                                                math::normalize(perpendicularToEdge));
+                                                perpendicularToEdge.normalized());
           minDistance = std::min(minDistance, distanceToEdge);
         }
       }
@@ -2040,7 +2043,7 @@ void voronoi_distance_to_edge(const float3 coord, const float randomness, float
 
 void voronoi_n_sphere_radius(const float3 coord, const float randomness, float *r_radius)
 {
-  const float3 cellPosition = math::floor(coord);
+  const float3 cellPosition = float3::floor(coord);
   const float3 localPosition = coord - cellPosition;
 
   float3 closestPoint = float3(0.0f, 0.0f, 0.0f);
@@ -2052,7 +2055,7 @@ void voronoi_n_sphere_radius(const float3 coord, const float randomness, float *
         const float3 cellOffset = float3(i, j, k);
         const float3 pointPosition = cellOffset +
                                      hash_float_to_float3(cellPosition + cellOffset) * randomness;
-        const float distanceToPoint = math::distance(pointPosition, localPosition);
+        const float distanceToPoint = float3::distance(pointPosition, localPosition);
         if (distanceToPoint < minDistance) {
           minDistance = distanceToPoint;
           closestPoint = pointPosition;
@@ -2073,7 +2076,7 @@ void voronoi_n_sphere_radius(const float3 coord, const float randomness, float *
         const float3 cellOffset = float3(i, j, k) + closestPointOffset;
         const float3 pointPosition = cellOffset +
                                      hash_float_to_float3(cellPosition + cellOffset) * randomness;
-        const float distanceToPoint = math::distance(closestPoint, pointPosition);
+        const float distanceToPoint = float3::distance(closestPoint, pointPosition);
         if (distanceToPoint < minDistance) {
           minDistance = distanceToPoint;
           closestPointToClosestPoint = pointPosition;
@@ -2081,7 +2084,7 @@ void voronoi_n_sphere_radius(const float3 coord, const float randomness, float *
       }
     }
   }
-  *r_radius = math::distance(closestPointToClosestPoint, closestPoint) / 2.0f;
+  *r_radius = float3::distance(closestPointToClosestPoint, closestPoint) / 2.0f;
 }
 
 /* **** 4D Voronoi **** */
@@ -2093,7 +2096,7 @@ static float voronoi_distance(const float4 a,
 {
   switch (metric) {
     case NOISE_SHD_VORONOI_EUCLIDEAN:
-      return math::distance(a, b);
+      return float4::distance(a, b);
     case NOISE_SHD_VORONOI_MANHATTAN:
       return fabsf(a.x - b.x) + fabsf(a.y - b.y) + fabsf(a.z - b.z) + fabsf(a.w - b.w);
     case NOISE_SHD_VORONOI_CHEBYCHEV:
@@ -2118,7 +2121,7 @@ void voronoi_f1(const float4 coord,
                 float3 *r_color,
                 float4 *r_position)
 {
-  const float4 cellPosition = math::floor(coord);
+  const float4 cellPosition = float4::floor(coord);
   const float4 localPosition = coord - cellPosition;
 
   float minDistance = 8.0f;
@@ -2163,7 +2166,7 @@ void voronoi_smooth_f1(const float4 coord,
                        float3 *r_color,
                        float4 *r_position)
 {
-  const float4 cellPosition = math::floor(coord);
+  const float4 cellPosition = float4::floor(coord);
   const float4 localPosition = coord - cellPosition;
   const float smoothness_clamped = max_ff(smoothness, FLT_MIN);
 
@@ -2188,10 +2191,10 @@ void voronoi_smooth_f1(const float4 coord,
             correctionFactor /= 1.0f + 3.0f * smoothness;
             if (r_color != nullptr) {
               const float3 cellColor = hash_float_to_float3(cellPosition + cellOffset);
-              smoothColor = math::interpolate(smoothColor, cellColor, h) - correctionFactor;
+              smoothColor = float3::interpolate(smoothColor, cellColor, h) - correctionFactor;
             }
             if (r_position != nullptr) {
-              smoothPosition = math::interpolate(smoothPosition, pointPosition, h) -
+              smoothPosition = float4::interpolate(smoothPosition, pointPosition, h) -
                                correctionFactor;
             }
           }
@@ -2218,7 +2221,7 @@ void voronoi_f2(const float4 coord,
                 float3 *r_color,
                 float4 *r_position)
 {
-  const float4 cellPosition = math::floor(coord);
+  const float4 cellPosition = float4::floor(coord);
   const float4 localPosition = coord - cellPosition;
 
   float distanceF1 = 8.0f;
@@ -2267,7 +2270,7 @@ void voronoi_f2(const float4 coord,
 
 void voronoi_distance_to_edge(const float4 coord, const float randomness, float *r_distance)
 {
-  const float4 cellPosition = math::floor(coord);
+  const float4 cellPosition = float4::floor(coord);
   const float4 localPosition = coord - cellPosition;
 
   float4 vectorToClosest = float4(0.0f, 0.0f, 0.0f, 0.0f);
@@ -2304,7 +2307,7 @@ void voronoi_distance_to_edge(const float4 coord, const float randomness, float
           const float4 perpendicularToEdge = vectorToPoint - vectorToClosest;
           if (dot_v4v4(perpendicularToEdge, perpendicularToEdge) > 0.0001f) {
             const float distanceToEdge = dot_v4v4((vectorToClosest + vectorToPoint) / 2.0f,
-                                                  math::normalize(perpendicularToEdge));
+                                                  float4::normalize(perpendicularToEdge));
             minDistance = std::min(minDistance, distanceToEdge);
           }
         }
@@ -2316,7 +2319,7 @@ void voronoi_distance_to_edge(const float4 coord, const float randomness, float
 
 void voronoi_n_sphere_radius(const float4 coord, const float randomness, float *r_radius)
 {
-  const float4 cellPosition = math::floor(coord);
+  const float4 cellPosition = float4::floor(coord);
   const float4 localPosition = coord - cellPosition;
 
   float4 closestPoint = float4(0.0f, 0.0f, 0.0f, 0.0f);
@@ -2330,7 +2333,7 @@ void voronoi_n_sphere_radius(const float4 coord, const float randomness, float *
           const float4 pointPosition = cellOffset +
                                        hash_float_to_float4(cellPosition + cellOffset) *
                                            randomness;
-          const float distanceToPoint = math::distance(pointPosition, localPosition);
+          const float distanceToPoint = float4::distance(pointPosition, localPosition);
           if (distanceToPoint < minDistance) {
             minDistance = distanceToPoint;
             closestPoint = pointPosition;
@@ -2354,7 +2357,7 @@ void voronoi_n_sphere_radius(const float4 coord, const float randomness, float *
           const float4 pointPosition = cellOffset +
                                        hash_float_to_float4(cellPosition + cellOffset) *
                                            randomness;
-          const float distanceToPoint = math::distance(closestPoint, pointPosition);
+          const float distanceToPoint = float4::distance(closestPoint, pointPosition);
           if (distanceToPoint < minDistance) {
             minDistance = distanceToPoint;
             closestPointToClosestPoint = pointPosition;
@@ -2363,7 +2366,7 @@ void voronoi_n_sphere_radius(const float4 coord, const float randomness, float *
       }
     }
   }
-  *r_radius = math::distance(closestPointToClosestPoint, closestPoint) / 2.0f;
+  *r_radius = float4::distance(closestPointToClosestPoint, closestPoint) / 2.0f;
 }
 
 /** \} */
diff --git a/source/blender/blenlib/tests/BLI_delaunay_2d_test.cc b/source/blender/blenlib/tests/BLI_delaunay_2d_test.cc
index eac3faa6d15..70e3a99e57a 100644
--- a/source/blender/blenlib/tests/BLI_delaunay_2d_test.cc
+++ b/source/blender/blenlib/tests/BLI_delaunay_2d_test.cc
@@ -21,9 +21,10 @@ extern "C" {
 #define DO_RANDOM_TESTS 0
 
 #include "BLI_array.hh"
+#include "BLI_double2.hh"
 #include "BLI_math_boolean.hh"
 #include "BLI_math_mpq.hh"
-#include "BLI_math_vec_mpq_types.hh"
+#include "BLI_mpq2.hh"
 #include "BLI_vector.hh"
 
 #include "BLI_delaunay_2d.h"
diff --git a/source/blender/blenlib/tests/BLI_math_vec_types_test.cc b/source/blender/blenlib/tests/BLI_math_vec_types_test.cc
deleted file mode 100644
index 8aa1f90fde2..00000000000
--- a/source/blender/blenlib/tests/BLI_math_vec_types_test.cc
+++ /dev/null
@@ -1,149 +0,0 @@
-/* Apache License, Version 2.0 */
-
-#include "testing/testing.h"
-
-#include "BLI_math_vec_types.hh"
-
-namespace blender::tests {
-
-using namespace blender::math;
-
-TEST(math_vec_types, ScalarConstructorUnsigned)
-{
-  float2 u(5u);
-  EXPECT_EQ(u[0], 5.0f);
-  EXPECT_EQ(u[1], 5.0f);
-}
-
-TEST(math_vec_types, ScalarConstructorInt)
-{
-  float2 i(-5);
-  EXPECT_EQ(i[0], -5.0f);
-  EXPECT_EQ(i[1], -5.0f);
-}
-
-TEST(math_vec_types, ScalarConstructorFloat)
-{
-  float2 f(5.2f);
-  EXPECT_FLOAT_EQ(f[0], 5.2f);
-  EXPECT_FLOAT_EQ(f[1], 5.2f);
-}
-
-TEST(math_vec_types, ScalarConstructorDouble)
-{
-  float2 d(5.2);
-  EXPECT_FLOAT_EQ(d[0], 5.2f);
-  EXPECT_FLOAT_EQ(d[1], 5.2f);
-}
-
-TEST(math_vec_types, MultiScalarConstructorVec2)
-{
-  int2 i(5.5f, -1.8);
-  EXPECT_EQ(i[0], 5);
-  EXPECT_EQ(i[1], -1);
-}
-
-TEST(math_vec_types, MultiScalarConstructorVec3)
-{
-  int3 i(5.5f, -1.8, 6u);
-  EXPECT_EQ(i[0], 5);
-  EXPECT_EQ(i[1], -1);
-  EXPECT_EQ(i[2], 6);
-}
-
-TEST(math_vec_types, MultiScalarConstructorVec4)
-{
-  int4 i(5.5f, -1.8, 6u, 0.888f);
-  EXPECT_EQ(i[0], 5);
-  EXPECT_EQ(i[1], -1);
-  EXPECT_EQ(i[2], 6);
-  EXPECT_EQ(i[3], 0);
-}
-
-TEST(math_vec_types, MixedScalarVectorConstructorVec3)
-{
-  float3 fl_v2(float2(5.5f), 1.8f);
-  EXPECT_FLOAT_EQ(fl_v2[0], 5.5f);
-  EXPECT_FLOAT_EQ(fl_v2[1], 5.5f);
-  EXPECT_FLOAT_EQ(fl_v2[2], 1.8f);
-
-  float3 v2_fl(1.8f, float2(5.5f));
-  EXPECT_FLOAT_EQ(v2_fl[0], 1.8f);
-  EXPECT_FLOAT_EQ(v2_fl[1], 5.5f);
-  EXPECT_FLOAT_EQ(v2_fl[2], 5.5f);
-}
-
-TEST(math_vec_types, MixedScalarVectorConstructorVec4)
-{
-  int4 v2_fl_fl(float2(1), 2, 3);
-  EXPECT_EQ(v2_fl_fl[0], 1);
-  EXPECT_EQ(v2_fl_fl[1], 1);
-  EXPECT_EQ(v2_fl_fl[2], 2);
-  EXPECT_EQ(v2_fl_fl[3], 3);
-
-  float4 fl_v2_fl(1, int2(2), 3);
-  EXPECT_EQ(fl_v2_fl[0], 1);
-  EXPECT_EQ(fl_v2_fl[1], 2);
-  EXPECT_EQ(fl_v2_fl[2], 2);
-  EXPECT_EQ(fl_v2_fl[3], 3);
-
-  double4 fl_fl_v2(1, 2, double2(3));
-  EXPECT_EQ(fl_fl_v2[0], 1);
-  EXPECT_EQ(fl_fl_v2[1], 2);
-  EXPECT_EQ(fl_fl_v2[2], 3);
-  EXPECT_EQ(fl_fl_v2[3], 3);
-
-  int4 v2_v2(float2(1), uint2(2));
-  EXPECT_EQ(v2_v2[0], 1);
-  EXPECT_EQ(v2_v2[1], 1);
-  EXPECT_EQ(v2_v2[2], 2);
-  EXPECT_EQ(v2_v2[3], 2);
-
-  float4 v3_fl(uint3(1), 2);
-  EXPECT_EQ(v3_fl[0], 1);
-  EXPECT_EQ(v3_fl[1], 1);
-  EXPECT_EQ(v3_fl[2], 1);
-  EXPECT_EQ(v3_fl[3], 2);
-
-  uint4 fl_v3(1, float3(2));
-  EXPECT_EQ(fl_v3[0], 1);
-  EXPECT_EQ(fl_v3[1], 2);
-  EXPECT_EQ(fl_v3[2], 2);
-  EXPECT_EQ(fl_v3[3], 2);
-}
-
-TEST(math_vec_types, ComponentMasking)
-{
-  int4 i(0, 1, 2, 3);
-  float2 f2 = float2(i);
-  EXPECT_EQ(f2[0], 0.0f);
-  EXPECT_EQ(f2[1], 1.0f);
-}
-
-TEST(math_vec_types, PointerConversion)
-{
-  float array[3] = {1.0f, 2.0f, 3.0f};
-  float3 farray(array);
-  EXPECT_EQ(farray[0], 1.0f);
-  EXPECT_EQ(farray[1], 2.0f);
-  EXPECT_EQ(farray[2], 3.0f);
-}
-
-TEST(math_vec_types, PointerArrayConversion)
-{
-  float array[1][3] = {{1.0f, 2.0f, 3.0f}};
-  float(*ptr)[3] = array;
-  float3 fptr(ptr);
-  EXPECT_EQ(fptr[0], 1.0f);
-  EXPECT_EQ(fptr[1], 2.0f);
-  EXPECT_EQ(fptr[2], 3.0f);
-}
-
-TEST(math_vec_types, VectorTypeConversion)
-{
-  double2 d(int2(float2(5.75f, -1.57f)));
-  EXPECT_EQ(d[0], 5.0);
-  EXPECT_EQ(d[1], -1.0);
-}
-
-}  // namespace blender::tests
diff --git a/source/blender/blenlib/tests/BLI_memory_utils_test.cc b/source/blender/blenlib/tests/BLI_memory_utils_test.cc
index 74e54151a06..207f310d902 100644
--- a/source/blender/blenlib/tests/BLI_memory_utils_test.cc
+++ b/source/blender/blenlib/tests/BLI_memory_utils_test.cc
@@ -1,6 +1,6 @@
 /* Apache License, Version 2.0 */
 
-#include "BLI_math_vec_types.hh"
+#include "BLI_float3.hh"
 #include "BLI_memory_utils.hh"
 #include "BLI_strict_flags.h"
 #include "testing/testing.h"
diff --git a/source/blender/blenlib/tests/BLI_mesh_boolean_test.cc b/source/blender/blenlib/tests/BLI_mesh_boolean_test.cc
index 2b8fb3dbea4..d759f0c3be4 100644
--- a/source/blender/blenlib/tests/BLI_mesh_boolean_test.cc
+++ b/source/blender/blenlib/tests/BLI_mesh_boolean_test.cc
@@ -11,8 +11,8 @@
 #include "BLI_array.hh"
 #include "BLI_map.hh"
 #include "BLI_math_mpq.hh"
-#include "BLI_math_vec_mpq_types.hh"
 #include "BLI_mesh_boolean.hh"
+#include "BLI_mpq3.hh"
 #include "BLI_vector.hh"
 
 #ifdef WITH_GMP
diff --git a/source/blender/blenlib/tests/BLI_mesh_intersect_test.cc b/source/blender/blenlib/tests/BLI_mesh_intersect_test.cc
index d2d76593129..68111fb8eb1 100644
--- a/source/blender/blenlib/tests/BLI_mesh_intersect_test.cc
+++ b/source/blender/blenlib/tests/BLI_mesh_intersect_test.cc
@@ -10,8 +10,8 @@
 
 #include "BLI_array.hh"
 #include "BLI_math_mpq.hh"
-#include "BLI_math_vec_mpq_types.hh"
 #include "BLI_mesh_intersect.hh"
+#include "BLI_mpq3.hh"
 #include "BLI_task.h"
 #include "BLI_vector.hh"