BLI: Refactor vector types & functions to use templates

This patch implements the vector types (i.e:`float2`) by making heavy
usage of templating. All vector functions are now outside of the vector
classes (inside the `blender::math` namespace) and are not vector size
dependent for the most part.

In the ongoing effort to make shaders less GL centric, we are aiming
to share more code between GLSL and C++ to avoid code duplication.

####Motivations:
 - We are aiming to share UBO and SSBO structures between GLSL and C++.
 This means we will use many of the existing vector types and others
 we currently don't have (uintX, intX). All these variations were
 asking for many more code duplication.
 - Deduplicate existing code which is duplicated for each vector size.
 - We also want to share small functions. Which means that vector
 functions should be static and not in the class namespace.
 - Reduce friction to use these types in new projects due to their
 incompleteness.
 - The current state of the `BLI_(float|double|mpq)(2|3|4).hh` is a
 bit of a let down. Most clases are incomplete, out of sync with each
 others with different codestyles, and some functions that should be
 static are not (i.e: `float3::reflect()`).

####Upsides:
 - Still support `.x, .y, .z, .w` for readability.
 - Compact, readable and easilly extendable.
 - All of the vector functions are available for all the vectors types
 and can be restricted to certain types. Also template specialization
 let us define exception for special class (like mpq).
 - With optimization ON, the compiler unroll the loops and performance
 is the same.

####Downsides:
 - Might impact debugability. Though I would arge that the bugs are
 rarelly caused by the vector class itself (since the operations are
 quite trivial) but by the type conversions.
 - Might impact compile time. I did not saw a significant impact since
 the usage is not really widespread.
 - Functions needs to be rewritten to support arbitrary vector length.
 For instance, one can't call `len_squared_v3v3` in
 `math::length_squared()` and call it a day.
 - Type cast does not work with the template version of the `math::`
 vector functions. Meaning you need to manually cast `float *` and
 `(float *)[3]` to `float3` for the function calls.
 i.e: `math::distance_squared(float3(nearest.co), positions[i]);`
 - Some parts might loose in readability:
 `float3::dot(v1.normalized(), v2.normalized())`
 becoming
 `math::dot(math::normalize(v1), math::normalize(v2))`
 But I propose, when appropriate, to use
 `using namespace blender::math;` on function local or file scope to
 increase readability.
 `dot(normalize(v1), normalize(v2))`

####Consideration:
 - Include back `.length()` method. It is quite handy and is more C++
 oriented.
 - I considered the GLM library as a candidate for replacement. It felt
 like too much for what we need and would be difficult to extend / modify
 to our needs.
 - I used Macros to reduce code in operators declaration and potential
 copy paste bugs. This could reduce debugability and could be reverted.
 - This touches `delaunay_2d.cc` and the intersection code. I would like
 to know @howardt opinion on the matter.
 - The `noexcept` on the copy constructor of `mpq(2|3)` is being removed.
 But according to @JacquesLucke it is not a real problem for now.

I would like to give a huge thanks to @JacquesLucke who helped during this
and pushed me to reduce the duplication further.

Reviewed By: brecht, sergey, JacquesLucke

Differential Revision: https://developer.blender.org/D13791

45 files changed, 1457 insertions, 1877 deletions

diff --git a/source/blender/blenlib/BLI_array_store.h b/source/blender/blenlib/BLI_array_store.h
index 68928f53e55..0be361d4ab9 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,
-                                       size_t data_len,
+                                       const 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 202ae9a9786..50fc2ce909b 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, size_t arr_stride, const uint *order, void *arr_temp);
+    void *arr, uint arr_len, const 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],
-                                   size_t elem_size,
+                                   const 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 e79d44fd934..7f577443cf7 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, size_t new_count);
+void BLI_buffer_resize(BLI_Buffer *buffer, const 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, size_t new_count);
  * - Ignored (malloc'd).
  * - Cleared (when #BLI_BUFFER_USE_CALLOC is set).
  */
-void BLI_buffer_reinit(BLI_Buffer *buffer, size_t new_count);
+void BLI_buffer_reinit(BLI_Buffer *buffer, const 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 db0df95499f..1ee0be64cee 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_mpq2.hh"
+#  include "BLI_math_vec_mpq_types.hh"
+#  include "BLI_math_vec_types.hh"
 #  include "BLI_vector.hh"
 
 namespace blender::meshintersect {
diff --git a/source/blender/blenlib/BLI_double2.hh b/source/blender/blenlib/BLI_double2.hh
deleted file mode 100644
index 0abff01ab2f..00000000000
--- a/source/blender/blenlib/BLI_double2.hh
+++ /dev/null
@@ -1,143 +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
- */
-
-#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
deleted file mode 100644
index ab258c9121b..00000000000
--- a/source/blender/blenlib/BLI_double3.hh
+++ /dev/null
@@ -1,246 +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
- */
-
-#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 28cb5f6d84b..0022823b3de 100644
--- a/source/blender/blenlib/BLI_fileops.h
+++ b/source/blender/blenlib/BLI_fileops.h
@@ -155,7 +155,8 @@ 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, size_t maxncpy) ATTR_WARN_UNUSED_RESULT ATTR_NONNULL();
+char *BLI_current_working_dir(char *dir, const 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
deleted file mode 100644
index bb4229db86e..00000000000
--- a/source/blender/blenlib/BLI_float2.hh
+++ /dev/null
@@ -1,218 +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
-
-#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
deleted file mode 100644
index 765f524fb31..00000000000
--- a/source/blender/blenlib/BLI_float3.hh
+++ /dev/null
@@ -1,320 +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
-
-#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
deleted file mode 100644
index 5b487f6d029..00000000000
--- a/source/blender/blenlib/BLI_float4.hh
+++ /dev/null
@@ -1,138 +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
-
-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 b7f839f4ddf..81c969d02d0 100644
--- a/source/blender/blenlib/BLI_float4x4.hh
+++ b/source/blender/blenlib/BLI_float4x4.hh
@@ -16,8 +16,9 @@
 
 #pragma once
 
-#include "BLI_float3.hh"
 #include "BLI_math_matrix.h"
+#include "BLI_math_vec_types.hh"
+#include "BLI_math_vector.h"
 
 namespace blender {
 
@@ -63,7 +64,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 = -float3::cross(forward, up);
+    const float3 cross = -math::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 a35c743c80b..8b32c09b56b 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(size_t elem_size);
+GSQueue *BLI_gsqueue_new(const 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 f73d1f22502..64852b95ae4 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,
-                              size_t bytes_size,
+                              const 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,
-                                        size_t string_offset,
+                                        const 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,
-                               size_t bytes_size,
+                               const 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 20fd00b2aa4..8cf93c82dec 100644
--- a/source/blender/blenlib/BLI_math_boolean.hh
+++ b/source/blender/blenlib/BLI_math_boolean.hh
@@ -21,13 +21,11 @@
  * \brief Math vector functions needed specifically for mesh intersect and boolean.
  */
 
-#include "BLI_double2.hh"
-#include "BLI_double3.hh"
+#include "BLI_math_vec_types.hh"
 
 #ifdef WITH_GMP
 #  include "BLI_math_mpq.hh"
-#  include "BLI_mpq2.hh"
-#  include "BLI_mpq3.hh"
+#  include "BLI_math_vec_mpq_types.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
new file mode 100644
index 00000000000..36eb0cac83c
--- /dev/null
+++ b/source/blender/blenlib/BLI_math_vec_mpq_types.hh
@@ -0,0 +1,91 @@
+/*
+ * 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
new file mode 100644
index 00000000000..52aacd294e4
--- /dev/null
+++ b/source/blender/blenlib/BLI_math_vec_types.hh
@@ -0,0 +1,566 @@
+/*
+ * 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
new file mode 100644
index 00000000000..e7d765df842
--- /dev/null
+++ b/source/blender/blenlib/BLI_math_vector.hh
@@ -0,0 +1,399 @@
+/*
+ * 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 bcfe2efc5e5..4ac4712bc8c 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(size_t bufsize,
+struct MemArena *BLI_memarena_new(const 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, size_t align) ATTR_NONNULL(1);
+void BLI_memarena_use_align(struct MemArena *ma, const 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 09d8f646905..79e25e26040 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, size_t size);
+bool BLI_memory_is_zero(const void *arr, const 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 37691017c12..9a5be79b61e 100644
--- a/source/blender/blenlib/BLI_memory_utils.hh
+++ b/source/blender/blenlib/BLI_memory_utils.hh
@@ -557,13 +557,4 @@ 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 71a8abb822f..0ebee6f16a8 100644
--- a/source/blender/blenlib/BLI_mesh_intersect.hh
+++ b/source/blender/blenlib/BLI_mesh_intersect.hh
@@ -28,12 +28,11 @@
 #  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_mpq3.hh"
+#  include "BLI_math_vec_mpq_types.hh"
+#  include "BLI_math_vec_types.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
deleted file mode 100644
index 18bc8821d9c..00000000000
--- a/source/blender/blenlib/BLI_mpq2.hh
+++ /dev/null
@@ -1,184 +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_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
deleted file mode 100644
index b9eda2ad7e1..00000000000
--- a/source/blender/blenlib/BLI_mpq3.hh
+++ /dev/null
@@ -1,297 +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 <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 4f68ef17ca2..297c65c250a 100644
--- a/source/blender/blenlib/BLI_noise.hh
+++ b/source/blender/blenlib/BLI_noise.hh
@@ -16,9 +16,7 @@
 
 #pragma once
 
-#include "BLI_float2.hh"
-#include "BLI_float3.hh"
-#include "BLI_float4.hh"
+#include "BLI_math_vec_types.hh"
 
 namespace blender::noise {
 
diff --git a/source/blender/blenlib/BLI_path_util.h b/source/blender/blenlib/BLI_path_util.h
index 658cc0c3825..16f479cb3b8 100644
--- a/source/blender/blenlib/BLI_path_util.h
+++ b/source/blender/blenlib/BLI_path_util.h
@@ -75,15 +75,16 @@ 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, size_t dirlen, size_t filelen);
+void BLI_split_dirfile(
+    const char *string, char *dir, char *file, const size_t dirlen, const 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, size_t dirlen);
+void BLI_split_dir_part(const char *string, char *dir, const 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, size_t filelen);
+void BLI_split_file_part(const char *string, char *file, const 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.
@@ -93,7 +94,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, size_t maxlen, const char *__restrict file)
+void BLI_path_append(char *__restrict dst, const size_t maxlen, const char *__restrict file)
     ATTR_NONNULL();
 /**
  * Simple appending of filename to dir, does not check for valid path!
@@ -103,7 +104,7 @@ void BLI_path_append(char *__restrict dst, size_t maxlen, const char *__restrict
  * that de-duplicates separators and can handle an arbitrary number of paths.
  */
 void BLI_join_dirfile(char *__restrict dst,
-                      size_t maxlen,
+                      const size_t maxlen,
                       const char *__restrict dir,
                       const char *__restrict file) ATTR_NONNULL();
 /**
@@ -113,7 +114,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, size_t dst_len, const char *path_first, ...)
+size_t BLI_path_join(char *__restrict dst, const size_t dst_len, const char *path_first, ...)
     ATTR_NONNULL(1, 3) ATTR_SENTINEL(0);
 /**
  * Like Python's `os.path.basename()`
@@ -163,12 +164,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, size_t maxlen);
+bool BLI_path_program_extensions_add_win32(char *name, const size_t maxlen);
 #endif
 /**
  * Search for a binary (executable)
  */
-bool BLI_path_program_search(char *fullname, size_t maxlen, const char *name);
+bool BLI_path_program_search(char *fullname, const size_t maxlen, const char *name);
 
 /**
  * \return true when `str` end with `ext` (case insensitive).
@@ -352,7 +353,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, size_t maxlen) ATTR_NONNULL();
+bool BLI_path_abs_from_cwd(char *path, const 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 cc9e9b374d7..667d6df8996 100644
--- a/source/blender/blenlib/BLI_rand.hh
+++ b/source/blender/blenlib/BLI_rand.hh
@@ -20,9 +20,8 @@
 
 #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 eb4e69a42d4..653f5f61c9e 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(size_t elem_size,
+BLI_Stack *BLI_stack_new_ex(const size_t elem_size,
                             const char *description,
-                            size_t chunk_size) ATTR_WARN_UNUSED_RESULT ATTR_NONNULL();
+                            const 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(size_t elem_size, const char *description) ATTR_WARN_UNUSED_RESULT
+BLI_Stack *BLI_stack_new(const 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 8177545911c..a82e97914db 100644
--- a/source/blender/blenlib/BLI_string.h
+++ b/source/blender/blenlib/BLI_string.h
@@ -42,7 +42,8 @@ extern "C" {
  * \param len: The number of bytes to duplicate
  * \retval Returns the duplicated string
  */
-char *BLI_strdupn(const char *str, size_t len) ATTR_MALLOC ATTR_WARN_UNUSED_RESULT ATTR_NONNULL();
+char *BLI_strdupn(const char *str, const size_t len) ATTR_MALLOC ATTR_WARN_UNUSED_RESULT
+    ATTR_NONNULL();
 
 /**
  * Duplicates the cstring \a str into a newly mallocN'd
@@ -73,7 +74,8 @@ char *BLI_strdupcat(const char *__restrict str1,
  * the size of dst)
  * \retval Returns dst
  */
-char *BLI_strncpy(char *__restrict dst, const char *__restrict src, size_t maxncpy) ATTR_NONNULL();
+char *BLI_strncpy(char *__restrict dst, const char *__restrict src, const size_t maxncpy)
+    ATTR_NONNULL();
 
 /**
  * Like BLI_strncpy but ensures dst is always padded by given char,
@@ -105,7 +107,7 @@ char *BLI_strncpy_ensure_pad(char *__restrict dst,
  */
 size_t BLI_strncpy_rlen(char *__restrict dst,
                         const char *__restrict src,
-                        size_t maxncpy) ATTR_WARN_UNUSED_RESULT ATTR_NONNULL();
+                        const 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();
@@ -184,7 +186,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,
-                                 size_t string_len,
+                                 const size_t string_len,
                                  const char *replace_table[][2],
                                  int replace_table_len);
 
@@ -233,7 +235,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, size_t dst_maxncpy)
+size_t BLI_str_escape(char *__restrict dst, const char *__restrict src, const size_t dst_maxncpy)
     ATTR_NONNULL();
 /**
  * This roughly matches C and Python's string escaping with double quotes - `"`.
@@ -249,9 +251,9 @@ size_t BLI_str_escape(char *__restrict dst, const char *__restrict src, size_t d
  */
 size_t BLI_str_unescape_ex(char *__restrict dst,
                            const char *__restrict src,
-                           size_t src_maxncpy,
+                           const size_t src_maxncpy,
                            /* Additional arguments. */
-                           size_t dst_maxncpy,
+                           const size_t dst_maxncpy,
                            bool *r_is_complete) ATTR_NONNULL();
 /**
  * See #BLI_str_unescape_ex doc-string.
@@ -263,7 +265,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, size_t src_maxncpy)
+size_t BLI_str_unescape(char *__restrict dst, const char *__restrict src, const size_t src_maxncpy)
     ATTR_NONNULL();
 
 /**
@@ -357,10 +359,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, size_t maxlen) ATTR_WARN_UNUSED_RESULT ATTR_NONNULL();
+size_t BLI_strnlen(const char *str, const size_t maxlen) ATTR_WARN_UNUSED_RESULT 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();
+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();
 /**
  * Strip white-space from end of the string.
  */
@@ -477,7 +479,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,
-                                size_t len,
+                                const 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 108a2f5fc7d..82622d442fb 100644
--- a/source/blender/blenlib/BLI_string_utf8.h
+++ b/source/blender/blenlib/BLI_string_utf8.h
@@ -110,12 +110,14 @@ 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, size_t outbuf_len) ATTR_NONNULL(2);
+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_as_utf32(char32_t *__restrict dst_w,
                              const char *__restrict src_c,
-                             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);
+                             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);
 /**
  * \return The UTF-32 len in UTF-8.
  */
@@ -160,20 +162,21 @@ 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, size_t maxlen, size_t *r_len_bytes)
+size_t BLI_strnlen_utf8_ex(const char *strc, const 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, size_t maxlen) ATTR_NONNULL(1) ATTR_WARN_UNUSED_RESULT;
+size_t BLI_strnlen_utf8(const char *strc, const 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,
-                                 size_t maxncpy) ATTR_NONNULL(1, 2);
+                                 const size_t maxncpy) ATTR_NONNULL(1, 2);
 size_t BLI_strncpy_wchar_from_utf8(wchar_t *__restrict dst,
                                    const char *__restrict src,
-                                   size_t maxncpy) ATTR_NONNULL(1, 2);
+                                   const 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 818bfe8182b..fd3918ff217 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, size_t str_len);
+void BLI_string_split_suffix(const char *string, char *r_body, char *r_suf, const size_t str_len);
 /**
  * `"a.b.c"` -> (`"a."`, `"b.c"`).
  */
-void BLI_string_split_prefix(const char *string, char *r_pre, char *r_body, size_t str_len);
+void BLI_string_split_prefix(const char *string, char *r_pre, char *r_body, const 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,
-                                 size_t name_len);
+                                 const 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 1cd18dc86ab..f0349e289ac 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,
-                                     size_t maxncpy,
+                                     const 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, size_t maxncpy, double time_seconds)
+size_t BLI_timecode_string_from_time_simple(char *str, const size_t maxncpy, double time_seconds)
     ATTR_NONNULL();
 
 /**
@@ -72,7 +72,7 @@ size_t BLI_timecode_string_from_time_simple(char *str, size_t maxncpy, double ti
  * \note in some cases this is used to print non-seconds values.
  */
 size_t BLI_timecode_string_from_time_seconds(char *str,
-                                             size_t maxncpy,
+                                             const 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 35d4158de59..9fe092fe525 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, size_t arr_size);
+extern bool BLI_memory_is_zero(const void *arr, const size_t arr_size);
 #endif
 
 #define MEMCMP_STRUCT_AFTER_IS_ZERO(struct_var, member) \
@@ -840,6 +840,15 @@ extern bool BLI_memory_is_zero(const void *arr, 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 3958fd8e2d2..90c6760019a 100644
--- a/source/blender/blenlib/CMakeLists.txt
+++ b/source/blender/blenlib/CMakeLists.txt
@@ -192,8 +192,6 @@ 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
@@ -207,9 +205,6 @@ 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
@@ -258,6 +253,8 @@ 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
@@ -267,8 +264,6 @@ 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
@@ -444,6 +439,7 @@ 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 03b0b11297b..90569d87c41 100644
--- a/source/blender/blenlib/intern/BLI_mempool_private.h
+++ b/source/blender/blenlib/intern/BLI_mempool_private.h
@@ -54,9 +54,8 @@ typedef struct ParallelMempoolTaskData {
  *
  * See #BLI_task_parallel_mempool implementation for detailed usage example.
  */
-ParallelMempoolTaskData *mempool_iter_threadsafe_create(BLI_mempool *pool,
-                                                        size_t num_iter) ATTR_WARN_UNUSED_RESULT
-    ATTR_NONNULL();
+ParallelMempoolTaskData *mempool_iter_threadsafe_create(BLI_mempool *pool, const 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 53e881a9fc7..842e6cb6135 100644
--- a/source/blender/blenlib/intern/delaunay_2d.cc
+++ b/source/blender/blenlib/intern/delaunay_2d.cc
@@ -25,11 +25,10 @@
 #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_mpq2.hh"
+#include "BLI_math_vec_mpq_types.hh"
 #include "BLI_set.hh"
 #include "BLI_task.hh"
 #include "BLI_vector.hh"
@@ -38,6 +37,8 @@
 
 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. */
 
@@ -788,11 +789,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 (vec2<mpq_class>::dot(exact_ab, exact_ac) < 0) {
+  if (dot(exact_ab, exact_ac) < 0) {
     return false;
   }
   vec2<mpq_class> exact_bc = c.exact - b.exact;
-  return vec2<mpq_class>::dot(exact_bc, exact_ac) >= 0;
+  return dot(exact_bc, exact_ac) >= 0;
 }
 #endif
 
@@ -801,11 +802,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 (vec2<double>::dot(ab, ac) < 0) {
+  if (dot(ab, ac) < 0) {
     return false;
   }
   vec2<double> bc = c.approx - b.approx;
-  return vec2<double>::dot(bc, ac) >= 0;
+  return dot(bc, ac) >= 0;
 }
 
 template<> CDTVert<double>::CDTVert(const vec2<double> &pt)
@@ -1081,7 +1082,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(vec2<T>::interpolate(*a, *b, lambda));
+  CDTVert<T> *v = this->add_vert(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];
@@ -1704,16 +1705,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 = vec2<T>::isect_seg_seg(va->co.exact, vb->co.exact, curco.exact, v2->co.exact);
+  auto isect = isect_seg_seg<vec2<T>>(va->co.exact, vb->co.exact, curco.exact, v2->co.exact);
   T &lambda = isect.lambda;
   switch (isect.kind) {
-    case vec2<T>::isect_result::LINE_LINE_CROSS: {
+    case isect_result<vec2<T>>::LINE_LINE_CROSS: {
 #ifdef WITH_GMP
       if (!std::is_same<T, mpq_class>::value) {
 #else
       if (true) {
 #endif
-        double len_ab = vec2<double>::distance(va->co.approx, vb->co.approx);
+        double len_ab = distance(va->co.approx, vb->co.approx);
         if (lambda * len_ab <= epsilon) {
           fill_crossdata_for_through_vert(va, se_vcva, cd, cd_next);
         }
@@ -1735,7 +1736,7 @@ void fill_crossdata_for_intersect(const FatCo<T> &curco,
       }
       break;
     }
-    case vec2<T>::isect_result::LINE_LINE_EXACT: {
+    case isect_result<vec2<T>>::LINE_LINE_EXACT: {
       if (lambda == 0) {
         fill_crossdata_for_through_vert(va, se_vcva, cd, cd_next);
       }
@@ -1750,7 +1751,7 @@ void fill_crossdata_for_intersect(const FatCo<T> &curco,
       }
       break;
     }
-    case vec2<T>::isect_result::LINE_LINE_NONE: {
+    case isect_result<vec2<T>>::LINE_LINE_NONE: {
 #ifdef WITH_GMP
       if (std::is_same<T, mpq_class>::value) {
         BLI_assert(false);
@@ -1766,9 +1767,9 @@ void fill_crossdata_for_intersect(const FatCo<T> &curco,
       }
       break;
     }
-    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)) {
+    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)) {
         fill_crossdata_for_through_vert(va, se_vcva, cd, cd_next);
       }
       else {
@@ -1845,7 +1846,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 = vec2<T>::interpolate(va->co.exact, vb->co.exact, cd->lambda);
+  vec2<T> curco = 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;
@@ -2386,7 +2387,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 = vec2<double>::distance_squared(co1, co2);
+      dissolvable_edges[i].len_squared = distance_squared(co1, co2);
       i++;
     }
   }
@@ -2569,18 +2570,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 = vec2<T>::isect_seg_seg(ray_end.exact,
+          auto isect = isect_seg_seg<vec2<T>>(ray_end.exact,
                                               mid.exact,
                                               e->symedges[0].vert->co.exact,
                                               e->symedges[1].vert->co.exact);
           switch (isect.kind) {
-            case vec2<T>::isect_result::LINE_LINE_CROSS: {
+            case isect_result<vec2<T>>::LINE_LINE_CROSS: {
               hits++;
               break;
             }
-            case vec2<T>::isect_result::LINE_LINE_EXACT:
-            case vec2<T>::isect_result::LINE_LINE_NONE:
-            case vec2<T>::isect_result::LINE_LINE_COLINEAR:
+            case isect_result<vec2<T>>::LINE_LINE_EXACT:
+            case isect_result<vec2<T>>::LINE_LINE_NONE:
+            case isect_result<vec2<T>>::LINE_LINE_COLINEAR:
               break;
           }
         }
diff --git a/source/blender/blenlib/intern/math_boolean.cc b/source/blender/blenlib/intern/math_boolean.cc
index c16755868aa..0bae3c23f79 100644
--- a/source/blender/blenlib/intern/math_boolean.cc
+++ b/source/blender/blenlib/intern/math_boolean.cc
@@ -18,15 +18,10 @@
  * \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_mpq2.hh"
-#include "BLI_mpq3.hh"
+#include "BLI_math_vec_types.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 223c0e273f0..6fab6c9a383 100644
--- a/source/blender/blenlib/intern/math_vec.cc
+++ b/source/blender/blenlib/intern/math_vec.cc
@@ -18,89 +18,83 @@
  * \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_mpq.hh"
-#include "BLI_mpq2.hh"
-#include "BLI_mpq3.hh"
+#include "BLI_math_vec_mpq_types.hh"
+#include "BLI_math_vector.hh"
 #include "BLI_span.hh"
 #include "BLI_utildefines.h"
 
-namespace blender {
+namespace blender::math {
 
-float2::isect_result float2::isect_seg_seg(const float2 &v1,
-                                           const float2 &v2,
-                                           const float2 &v3,
-                                           const float2 &v4)
+template<>
+isect_result<float2> isect_seg_seg(const float2 &v1,
+                                   const float2 &v2,
+                                   const float2 &v3,
+                                   const float2 &v4)
 {
-  float2::isect_result ans;
+  isect_result<float2> 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.mu = 0.0f;
-    ans.kind = float2::isect_result::LINE_LINE_COLINEAR;
+    ans.kind = isect_result<float2>::LINE_LINE_COLINEAR;
   }
   else {
     ans.lambda = ((v1[1] - v3[1]) * (v4[0] - v3[0]) - (v1[0] - v3[0]) * (v4[1] - v3[1])) / div;
-    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;
+    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;
       }
       else {
-        ans.kind = float2::isect_result::LINE_LINE_CROSS;
+        ans.kind = isect_result<float2>::LINE_LINE_CROSS;
       }
     }
     else {
-      ans.kind = float2::isect_result::LINE_LINE_NONE;
+      ans.kind = isect_result<float2>::LINE_LINE_NONE;
     }
   }
   return ans;
 }
 
-double2::isect_result double2::isect_seg_seg(const double2 &v1,
-                                             const double2 &v2,
-                                             const double2 &v3,
-                                             const double2 &v4)
+template<>
+isect_result<double2> isect_seg_seg(const double2 &v1,
+                                    const double2 &v2,
+                                    const double2 &v3,
+                                    const double2 &v4)
 {
-  double2::isect_result ans;
+  isect_result<double2> 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 = double2::isect_result::LINE_LINE_COLINEAR;
+    ans.kind = isect_result<double2>::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 = double2::isect_result::LINE_LINE_EXACT;
+        ans.kind = isect_result<double2>::LINE_LINE_EXACT;
       }
       else {
-        ans.kind = double2::isect_result::LINE_LINE_CROSS;
+        ans.kind = isect_result<double2>::LINE_LINE_CROSS;
       }
     }
     else {
-      ans.kind = double2::isect_result::LINE_LINE_NONE;
+      ans.kind = isect_result<double2>::LINE_LINE_NONE;
     }
   }
   return ans;
 }
 
 #ifdef WITH_GMP
-mpq2::isect_result mpq2::isect_seg_seg(const mpq2 &v1,
-                                       const mpq2 &v2,
-                                       const mpq2 &v3,
-                                       const mpq2 &v4)
+template<>
+isect_result<mpq2> isect_seg_seg(const mpq2 &v1, const mpq2 &v2, const mpq2 &v3, const mpq2 &v4)
 {
-  mpq2::isect_result ans;
+  isect_result<mpq2> 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 = mpq2::isect_result::LINE_LINE_COLINEAR;
+    ans.kind = isect_result<mpq2>::LINE_LINE_COLINEAR;
   }
   else {
     ans.lambda = ((v1[1] - v3[1]) * (v4[0] - v3[0]) - (v1[0] - v3[0]) * (v4[1] - v3[1])) / div;
@@ -109,66 +103,21 @@ mpq2::isect_result mpq2::isect_seg_seg(const mpq2 &v1,
     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 = mpq2::isect_result::LINE_LINE_EXACT;
+        ans.kind = isect_result<mpq2>::LINE_LINE_EXACT;
       }
       else {
-        ans.kind = mpq2::isect_result::LINE_LINE_CROSS;
+        ans.kind = isect_result<mpq2>::LINE_LINE_CROSS;
       }
     }
     else {
-      ans.kind = mpq2::isect_result::LINE_LINE_NONE;
+      ans.kind = isect_result<mpq2>::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)
 {
@@ -176,20 +125,6 @@ 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
+}  // namespace blender::math
diff --git a/source/blender/blenlib/intern/mesh_boolean.cc b/source/blender/blenlib/intern/mesh_boolean.cc
index ce4db0c6b9d..a3eae1896d3 100644
--- a/source/blender/blenlib/intern/mesh_boolean.cc
+++ b/source/blender/blenlib/intern/mesh_boolean.cc
@@ -28,8 +28,6 @@
 #  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"
@@ -37,8 +35,9 @@
 #  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"
@@ -1633,13 +1632,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 = mpq3::cross(abscissa, ordinate);
+  mpq3 normal = math::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 = normal.length_squared();
+  mpq_class nlen2 = math::length_squared(normal);
   mpq_class max_abs_slope = -1;
   Edge esort;
   const Vector<Edge> &edges = tmtopo.vert_edges(closestp);
@@ -1648,12 +1647,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 - (mpq3::dot(evec, normal) / nlen2) * normal;
+    mpq3 proj_evec = evec - (math::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 = mpq3::dot(proj_evec, abscissa);
-    mpq_class evec_o = mpq3::dot(proj_evec, ordinate);
+    mpq_class evec_a = math::dot(proj_evec, abscissa);
+    mpq_class evec_o = math::dot(proj_evec, ordinate);
     if (dbg_level > 0) {
       std::cout << "e = " << e << "\n";
       std::cout << "v_other = " << v_other << "\n";
@@ -1791,8 +1790,8 @@ static mpq_class closest_on_tri_to_point(const mpq3 &p,
   ap = p;
   ap -= a;
 
-  mpq_class d1 = mpq3::dot_with_buffer(ab, ap, m);
-  mpq_class d2 = mpq3::dot_with_buffer(ac, ap, m);
+  mpq_class d1 = math::dot_with_buffer(ab, ap, m);
+  mpq_class d2 = math::dot_with_buffer(ac, ap, m);
   if (d1 <= 0 && d2 <= 0) {
     /* Barycentric coordinates (1,0,0). */
     *r_edge = -1;
@@ -1800,13 +1799,13 @@ static mpq_class closest_on_tri_to_point(const mpq3 &p,
     if (dbg_level > 0) {
       std::cout << "  answer = a\n";
     }
-    return mpq3::distance_squared_with_buffer(p, a, m);
+    return math::distance_squared_with_buffer(p, a, m);
   }
   /* Check if p in vertex region outside b. */
   bp = p;
   bp -= b;
-  mpq_class d3 = mpq3::dot_with_buffer(ab, bp, m);
-  mpq_class d4 = mpq3::dot_with_buffer(ac, bp, m);
+  mpq_class d3 = math::dot_with_buffer(ab, bp, m);
+  mpq_class d4 = math::dot_with_buffer(ac, bp, m);
   if (d3 >= 0 && d4 <= d3) {
     /* Barycentric coordinates (0,1,0). */
     *r_edge = -1;
@@ -1814,7 +1813,7 @@ static mpq_class closest_on_tri_to_point(const mpq3 &p,
     if (dbg_level > 0) {
       std::cout << "  answer = b\n";
     }
-    return mpq3::distance_squared_with_buffer(p, b, m);
+    return math::distance_squared_with_buffer(p, b, m);
   }
   /* Check if p in region of ab. */
   mpq_class vc = d1 * d4 - d3 * d2;
@@ -1829,13 +1828,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 mpq3::distance_squared_with_buffer(p, r, m);
+    return math::distance_squared_with_buffer(p, r, m);
   }
   /* Check if p in vertex region outside c. */
   cp = p;
   cp -= c;
-  mpq_class d5 = mpq3::dot_with_buffer(ab, cp, m);
-  mpq_class d6 = mpq3::dot_with_buffer(ac, cp, m);
+  mpq_class d5 = math::dot_with_buffer(ab, cp, m);
+  mpq_class d6 = math::dot_with_buffer(ac, cp, m);
   if (d6 >= 0 && d5 <= d6) {
     /* Barycentric coordinates (0,0,1). */
     *r_edge = -1;
@@ -1843,7 +1842,7 @@ static mpq_class closest_on_tri_to_point(const mpq3 &p,
     if (dbg_level > 0) {
       std::cout << "  answer = c\n";
     }
-    return mpq3::distance_squared_with_buffer(p, c, m);
+    return math::distance_squared_with_buffer(p, c, m);
   }
   /* Check if p in edge region of ac. */
   mpq_class vb = d5 * d2 - d1 * d6;
@@ -1858,7 +1857,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 mpq3::distance_squared_with_buffer(p, r, m);
+    return math::distance_squared_with_buffer(p, r, m);
   }
   /* Check if p in edge region of bc. */
   mpq_class va = d3 * d6 - d5 * d4;
@@ -1874,7 +1873,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 mpq3::distance_squared_with_buffer(p, r, m);
+    return math::distance_squared_with_buffer(p, r, m);
   }
   /* p inside face region. Compute barycentric coordinates (u,v,w). */
   mpq_class denom = 1 / (va + vb + vc);
@@ -1890,7 +1889,7 @@ static mpq_class closest_on_tri_to_point(const mpq3 &p,
   if (dbg_level > 0) {
     std::cout << "  answer = inside at " << r << "\n";
   }
-  return mpq3::distance_squared_with_buffer(p, r, m);
+  return math::distance_squared_with_buffer(p, r, m);
 }
 
 static float closest_on_tri_to_point_float_dist_squared(const float3 &p,
@@ -2610,7 +2609,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 = tri_test.plane->norm.normalized();
+  double3 norm = math::normalize(tri_test.plane->norm);
   const double offset_amount = 1e-5;
   double3 offset_test_point = test_point + offset_amount * norm;
   if (dbg_level > 0) {
@@ -3002,7 +3001,7 @@ static void init_face_merge_state(FaceMergeState *fms,
       std::cout << "process tri = " << &tri << "\n";
     }
     BLI_assert(tri.plane_populated());
-    if (double3::dot(norm, tri.plane->norm) <= 0.0) {
+    if (math::dot(norm, tri.plane->norm) <= 0.0) {
       if (dbg_level > 0) {
         std::cout << "triangle has wrong orientation, skipping\n";
       }
@@ -3027,7 +3026,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 = vec.length_squared();
+        new_me.len_squared = math::length_squared(vec);
         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);
@@ -3267,7 +3266,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 && double3::dot(first_tri_normal, second_tri_normal) > 0.0) {
+  if (tris.size() == 2 && math::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]);
@@ -3332,7 +3331,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 = double3::dot(vec1.normalized(), vec2.normalized());
+  double cos_ang = math::dot(math::normalize(vec1), math::normalize(vec2));
   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 1f150137ba3..982759ffcff 100644
--- a/source/blender/blenlib/intern/mesh_intersect.cc
+++ b/source/blender/blenlib/intern/mesh_intersect.cc
@@ -30,15 +30,13 @@
 #  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_mpq2.hh"
-#  include "BLI_mpq3.hh"
+#  include "BLI_math_vec_mpq_types.hh"
+#  include "BLI_math_vec_types.hh"
 #  include "BLI_polyfill_2d.h"
 #  include "BLI_set.hh"
 #  include "BLI_span.hh"
@@ -198,14 +196,14 @@ void Face::populate_plane(bool need_exact)
       for (int i : index_range()) {
         co[i] = vert[i]->co_exact;
       }
-      normal_exact = mpq3::cross_poly(co);
+      normal_exact = math::cross_poly(co.as_span());
     }
     else {
       mpq3 tr02 = vert[0]->co_exact - vert[2]->co_exact;
       mpq3 tr12 = vert[1]->co_exact - vert[2]->co_exact;
-      normal_exact = mpq3::cross(tr02, tr12);
+      normal_exact = math::cross(tr02, tr12);
     }
-    mpq_class d_exact = -mpq3::dot(normal_exact, vert[0]->co_exact);
+    mpq_class d_exact = -math::dot(normal_exact, vert[0]->co_exact);
     plane = new Plane(normal_exact, d_exact);
   }
   else {
@@ -215,14 +213,14 @@ void Face::populate_plane(bool need_exact)
       for (int i : index_range()) {
         co[i] = vert[i]->co;
       }
-      normal = double3::cross_poly(co);
+      normal = math::cross_poly(co.as_span());
     }
     else {
       double3 tr02 = vert[0]->co - vert[2]->co;
       double3 tr12 = vert[1]->co - vert[2]->co;
-      normal = double3::cross_high_precision(tr02, tr12);
+      normal = math::cross(tr02, tr12);
     }
-    double d = -double3::dot(normal, vert[0]->co);
+    double d = -math::dot(normal, vert[0]->co);
     plane = new Plane(normal, d);
   }
 }
@@ -1098,15 +1096,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 = double3::abs(a);
-  double3 abs_b = double3::abs(b);
+  double3 abs_a = math::abs(a);
+  double3 abs_b = math::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 double3::dot(c, c);
+  return math::dot(c, c);
 }
 
 /* The index of dot when inputs are plane_coords with index 1 is much higher.
@@ -1143,11 +1141,11 @@ static int filter_plane_side(const double3 &p,
                              const double3 &abs_plane_p,
                              const double3 &abs_plane_no)
 {
-  double d = double3::dot(p - plane_p, plane_no);
+  double d = math::dot(p - plane_p, plane_no);
   if (d == 0.0) {
     return 0;
   }
-  double supremum = double3::dot(abs_p + abs_plane_p, abs_plane_no);
+  double supremum = math::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;
@@ -1178,9 +1176,9 @@ static inline mpq3 tti_interp(
   ab -= b;
   ac = a;
   ac -= c;
-  mpq_class den = mpq3::dot_with_buffer(ab, n, dotbuf);
+  mpq_class den = math::dot_with_buffer(ab, n, dotbuf);
   BLI_assert(den != 0);
-  mpq_class alpha = mpq3::dot_with_buffer(ac, n, dotbuf) / den;
+  mpq_class alpha = math::dot_with_buffer(ac, n, dotbuf) / den;
   return a - alpha * ab;
 }
 
@@ -1209,7 +1207,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(mpq3::dot_with_buffer(ad, n, dotbuf));
+  return sgn(math::dot_with_buffer(ad, n, dotbuf));
 }
 
 /**
@@ -1428,11 +1426,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 = 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);
+  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);
 
   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);
@@ -1448,9 +1446,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 = double3::abs(d_p2);
-  const double3 &abs_d_q2 = double3::abs(d_q2);
-  const double3 &abs_d_n1 = double3::abs(d_n1);
+  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);
 
   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);
@@ -1477,17 +1475,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(mpq3::dot_with_buffer(buf[0], n2, buf[1]));
+    sp1 = sgn(math::dot_with_buffer(buf[0], n2, buf[1]));
   }
   if (sq1 == 0) {
     buf[0] = q1;
     buf[0] -= r2;
-    sq1 = sgn(mpq3::dot_with_buffer(buf[0], n2, buf[1]));
+    sq1 = sgn(math::dot_with_buffer(buf[0], n2, buf[1]));
   }
   if (sr1 == 0) {
     buf[0] = r1;
     buf[0] -= r2;
-    sr1 = sgn(mpq3::dot_with_buffer(buf[0], n2, buf[1]));
+    sr1 = sgn(math::dot_with_buffer(buf[0], n2, buf[1]));
   }
 
   if (dbg_level > 1) {
@@ -1509,17 +1507,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(mpq3::dot_with_buffer(buf[0], n1, buf[1]));
+    sp2 = sgn(math::dot_with_buffer(buf[0], n1, buf[1]));
   }
   if (sq2 == 0) {
     buf[0] = q2;
     buf[0] -= r1;
-    sq2 = sgn(mpq3::dot_with_buffer(buf[0], n1, buf[1]));
+    sq2 = sgn(math::dot_with_buffer(buf[0], n1, buf[1]));
   }
   if (sr2 == 0) {
     buf[0] = r2;
     buf[0] -= r1;
-    sr2 = sgn(mpq3::dot_with_buffer(buf[0], n1, buf[1]));
+    sr2 = sgn(math::dot_with_buffer(buf[0], n1, buf[1]));
   }
 
   if (dbg_level > 1) {
@@ -1721,7 +1719,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 = mpq3::dominant_axis(ans.t_plane->norm_exact);
+  ans.proj_axis = math::dominant_axis(ans.t_plane->norm_exact);
   prepare_need_tri(ans, tm, t);
   for (const ITT_value &itt : itts) {
     switch (itt.kind) {
@@ -1757,7 +1755,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 = mpq3::dominant_axis(ans.t_plane->norm_exact);
+  ans.proj_axis = math::dominant_axis(ans.t_plane->norm_exact);
   for (const int t : cl) {
     prepare_need_tri(ans, tm, t);
   }
@@ -2004,9 +2002,9 @@ static bool is_quad_flip_first_third(const double3 &v1,
                                      const double3 &normal)
 {
   double3 dir_v3v1 = v3 - v1;
-  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);
+  double3 tangent = math::cross(dir_v3v1, normal);
+  double dot = math::dot(v1, tangent);
+  return (math::dot(v4, tangent) >= dot) || (math::dot(v2, tangent) <= dot);
 }
 
 /**
@@ -2124,7 +2122,7 @@ static Array<Face *> exact_triangulate_poly(Face *f, IMeshArena *arena)
     f->populate_plane(false);
   }
   const double3 &poly_normal = f->plane->norm;
-  int axis = double3::dominant_axis(poly_normal);
+  int axis = math::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
@@ -2203,15 +2201,15 @@ static bool face_is_degenerate(const Face *f)
   }
   double3 da = v2->co - v0->co;
   double3 db = v2->co - v1->co;
-  double3 dab = double3::cross_high_precision(da, db);
-  double dab_length_squared = dab.length_squared();
+  double3 dab = math::cross(da, db);
+  double dab_length_squared = math::length_squared(dab);
   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 = mpq3::cross(a, b);
+  mpq3 ab = math::cross(a, b);
   if (ab.x == 0 && ab.y == 0 && ab.z == 0) {
     return true;
   }
@@ -2231,8 +2229,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 = da.length_squared();
-    double db_length_squared = db.length_squared();
+    double da_length_squared = math::length_squared(da);
+    double db_length_squared = math::length_squared(db);
     if (da_length_squared == 0.0 || db_length_squared == 0.0) {
       return true;
     }
@@ -2240,8 +2238,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 = double3::cross_high_precision(da, db);
-    double dab_length_squared = dab.length_squared();
+    double3 dab = math::cross(da, db);
+    double dab_length_squared = math::length_squared(dab);
     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 a6ad18801fd..3460c1284fc 100644
--- a/source/blender/blenlib/intern/noise.cc
+++ b/source/blender/blenlib/intern/noise.cc
@@ -50,9 +50,7 @@
 #include <cmath>
 #include <cstdint>
 
-#include "BLI_float2.hh"
-#include "BLI_float3.hh"
-#include "BLI_float4.hh"
+#include "BLI_math_vec_types.hh"
 #include "BLI_math_base_safe.h"
 #include "BLI_noise.hh"
 #include "BLI_utildefines.h"
@@ -1469,7 +1467,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 = float3::interpolate(smoothColor, cellColor, h) - correctionFactor;
+        smoothColor = math::interpolate(smoothColor, cellColor, h) - correctionFactor;
       }
       if (r_w != nullptr) {
         smoothPosition = mix(smoothPosition, pointPosition, h) - correctionFactor;
@@ -1592,7 +1590,7 @@ static float voronoi_distance(const float2 a,
 {
   switch (metric) {
     case NOISE_SHD_VORONOI_EUCLIDEAN:
-      return float2::distance(a, b);
+      return math::distance(a, b);
     case NOISE_SHD_VORONOI_MANHATTAN:
       return fabsf(a.x - b.x) + fabsf(a.y - b.y);
     case NOISE_SHD_VORONOI_CHEBYCHEV:
@@ -1615,7 +1613,7 @@ void voronoi_f1(const float2 coord,
                 float3 *r_color,
                 float2 *r_position)
 {
-  const float2 cellPosition = float2::floor(coord);
+  const float2 cellPosition = math::floor(coord);
   const float2 localPosition = coord - cellPosition;
 
   float minDistance = 8.0f;
@@ -1654,7 +1652,7 @@ void voronoi_smooth_f1(const float2 coord,
                        float3 *r_color,
                        float2 *r_position)
 {
-  const float2 cellPosition = float2::floor(coord);
+  const float2 cellPosition = math::floor(coord);
   const float2 localPosition = coord - cellPosition;
   const float smoothness_clamped = max_ff(smoothness, FLT_MIN);
 
@@ -1676,11 +1674,10 @@ 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 = float3::interpolate(smoothColor, cellColor, h) - correctionFactor;
+          smoothColor = math::interpolate(smoothColor, cellColor, h) - correctionFactor;
         }
         if (r_position != nullptr) {
-          smoothPosition = float2::interpolate(smoothPosition, pointPosition, h) -
-                           correctionFactor;
+          smoothPosition = math::interpolate(smoothPosition, pointPosition, h) - correctionFactor;
         }
       }
     }
@@ -1704,7 +1701,7 @@ void voronoi_f2(const float2 coord,
                 float3 *r_color,
                 float2 *r_position)
 {
-  const float2 cellPosition = float2::floor(coord);
+  const float2 cellPosition = math::floor(coord);
   const float2 localPosition = coord - cellPosition;
 
   float distanceF1 = 8.0f;
@@ -1748,7 +1745,7 @@ void voronoi_f2(const float2 coord,
 
 void voronoi_distance_to_edge(const float2 coord, const float randomness, float *r_distance)
 {
-  const float2 cellPosition = float2::floor(coord);
+  const float2 cellPosition = math::floor(coord);
   const float2 localPosition = coord - cellPosition;
 
   float2 vectorToClosest = float2(0.0f, 0.0f);
@@ -1777,7 +1774,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,
-                                              perpendicularToEdge.normalized());
+                                              math::normalize(perpendicularToEdge));
         minDistance = std::min(minDistance, distanceToEdge);
       }
     }
@@ -1787,7 +1784,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 = float2::floor(coord);
+  const float2 cellPosition = math::floor(coord);
   const float2 localPosition = coord - cellPosition;
 
   float2 closestPoint = float2(0.0f, 0.0f);
@@ -1798,7 +1795,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 = float2::distance(pointPosition, localPosition);
+      const float distanceToPoint = math::distance(pointPosition, localPosition);
       if (distanceToPoint < minDistance) {
         minDistance = distanceToPoint;
         closestPoint = pointPosition;
@@ -1817,14 +1814,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 = float2::distance(closestPoint, pointPosition);
+      const float distanceToPoint = math::distance(closestPoint, pointPosition);
       if (distanceToPoint < minDistance) {
         minDistance = distanceToPoint;
         closestPointToClosestPoint = pointPosition;
       }
     }
   }
-  *r_radius = float2::distance(closestPointToClosestPoint, closestPoint) / 2.0f;
+  *r_radius = math::distance(closestPointToClosestPoint, closestPoint) / 2.0f;
 }
 
 /* **** 3D Voronoi **** */
@@ -1836,7 +1833,7 @@ static float voronoi_distance(const float3 a,
 {
   switch (metric) {
     case NOISE_SHD_VORONOI_EUCLIDEAN:
-      return float3::distance(a, b);
+      return math::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:
@@ -1860,7 +1857,7 @@ void voronoi_f1(const float3 coord,
                 float3 *r_color,
                 float3 *r_position)
 {
-  const float3 cellPosition = float3::floor(coord);
+  const float3 cellPosition = math::floor(coord);
   const float3 localPosition = coord - cellPosition;
 
   float minDistance = 8.0f;
@@ -1902,7 +1899,7 @@ void voronoi_smooth_f1(const float3 coord,
                        float3 *r_color,
                        float3 *r_position)
 {
-  const float3 cellPosition = float3::floor(coord);
+  const float3 cellPosition = math::floor(coord);
   const float3 localPosition = coord - cellPosition;
   const float smoothness_clamped = max_ff(smoothness, FLT_MIN);
 
@@ -1925,10 +1922,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 = float3::interpolate(smoothColor, cellColor, h) - correctionFactor;
+            smoothColor = math::interpolate(smoothColor, cellColor, h) - correctionFactor;
           }
           if (r_position != nullptr) {
-            smoothPosition = float3::interpolate(smoothPosition, pointPosition, h) -
+            smoothPosition = math::interpolate(smoothPosition, pointPosition, h) -
                              correctionFactor;
           }
         }
@@ -1954,7 +1951,7 @@ void voronoi_f2(const float3 coord,
                 float3 *r_color,
                 float3 *r_position)
 {
-  const float3 cellPosition = float3::floor(coord);
+  const float3 cellPosition = math::floor(coord);
   const float3 localPosition = coord - cellPosition;
 
   float distanceF1 = 8.0f;
@@ -2000,7 +1997,7 @@ void voronoi_f2(const float3 coord,
 
 void voronoi_distance_to_edge(const float3 coord, const float randomness, float *r_distance)
 {
-  const float3 cellPosition = float3::floor(coord);
+  const float3 cellPosition = math::floor(coord);
   const float3 localPosition = coord - cellPosition;
 
   float3 vectorToClosest = float3(0.0f, 0.0f, 0.0f);
@@ -2032,7 +2029,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,
-                                                perpendicularToEdge.normalized());
+                                                math::normalize(perpendicularToEdge));
           minDistance = std::min(minDistance, distanceToEdge);
         }
       }
@@ -2043,7 +2040,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 = float3::floor(coord);
+  const float3 cellPosition = math::floor(coord);
   const float3 localPosition = coord - cellPosition;
 
   float3 closestPoint = float3(0.0f, 0.0f, 0.0f);
@@ -2055,7 +2052,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 = float3::distance(pointPosition, localPosition);
+        const float distanceToPoint = math::distance(pointPosition, localPosition);
         if (distanceToPoint < minDistance) {
           minDistance = distanceToPoint;
           closestPoint = pointPosition;
@@ -2076,7 +2073,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 = float3::distance(closestPoint, pointPosition);
+        const float distanceToPoint = math::distance(closestPoint, pointPosition);
         if (distanceToPoint < minDistance) {
           minDistance = distanceToPoint;
           closestPointToClosestPoint = pointPosition;
@@ -2084,7 +2081,7 @@ void voronoi_n_sphere_radius(const float3 coord, const float randomness, float *
       }
     }
   }
-  *r_radius = float3::distance(closestPointToClosestPoint, closestPoint) / 2.0f;
+  *r_radius = math::distance(closestPointToClosestPoint, closestPoint) / 2.0f;
 }
 
 /* **** 4D Voronoi **** */
@@ -2096,7 +2093,7 @@ static float voronoi_distance(const float4 a,
 {
   switch (metric) {
     case NOISE_SHD_VORONOI_EUCLIDEAN:
-      return float4::distance(a, b);
+      return math::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:
@@ -2121,7 +2118,7 @@ void voronoi_f1(const float4 coord,
                 float3 *r_color,
                 float4 *r_position)
 {
-  const float4 cellPosition = float4::floor(coord);
+  const float4 cellPosition = math::floor(coord);
   const float4 localPosition = coord - cellPosition;
 
   float minDistance = 8.0f;
@@ -2166,7 +2163,7 @@ void voronoi_smooth_f1(const float4 coord,
                        float3 *r_color,
                        float4 *r_position)
 {
-  const float4 cellPosition = float4::floor(coord);
+  const float4 cellPosition = math::floor(coord);
   const float4 localPosition = coord - cellPosition;
   const float smoothness_clamped = max_ff(smoothness, FLT_MIN);
 
@@ -2191,10 +2188,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 = float3::interpolate(smoothColor, cellColor, h) - correctionFactor;
+              smoothColor = math::interpolate(smoothColor, cellColor, h) - correctionFactor;
             }
             if (r_position != nullptr) {
-              smoothPosition = float4::interpolate(smoothPosition, pointPosition, h) -
+              smoothPosition = math::interpolate(smoothPosition, pointPosition, h) -
                                correctionFactor;
             }
           }
@@ -2221,7 +2218,7 @@ void voronoi_f2(const float4 coord,
                 float3 *r_color,
                 float4 *r_position)
 {
-  const float4 cellPosition = float4::floor(coord);
+  const float4 cellPosition = math::floor(coord);
   const float4 localPosition = coord - cellPosition;
 
   float distanceF1 = 8.0f;
@@ -2270,7 +2267,7 @@ void voronoi_f2(const float4 coord,
 
 void voronoi_distance_to_edge(const float4 coord, const float randomness, float *r_distance)
 {
-  const float4 cellPosition = float4::floor(coord);
+  const float4 cellPosition = math::floor(coord);
   const float4 localPosition = coord - cellPosition;
 
   float4 vectorToClosest = float4(0.0f, 0.0f, 0.0f, 0.0f);
@@ -2307,7 +2304,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,
-                                                  float4::normalize(perpendicularToEdge));
+                                                  math::normalize(perpendicularToEdge));
             minDistance = std::min(minDistance, distanceToEdge);
           }
         }
@@ -2319,7 +2316,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 = float4::floor(coord);
+  const float4 cellPosition = math::floor(coord);
   const float4 localPosition = coord - cellPosition;
 
   float4 closestPoint = float4(0.0f, 0.0f, 0.0f, 0.0f);
@@ -2333,7 +2330,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 = float4::distance(pointPosition, localPosition);
+          const float distanceToPoint = math::distance(pointPosition, localPosition);
           if (distanceToPoint < minDistance) {
             minDistance = distanceToPoint;
             closestPoint = pointPosition;
@@ -2357,7 +2354,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 = float4::distance(closestPoint, pointPosition);
+          const float distanceToPoint = math::distance(closestPoint, pointPosition);
           if (distanceToPoint < minDistance) {
             minDistance = distanceToPoint;
             closestPointToClosestPoint = pointPosition;
@@ -2366,7 +2363,7 @@ void voronoi_n_sphere_radius(const float4 coord, const float randomness, float *
       }
     }
   }
-  *r_radius = float4::distance(closestPointToClosestPoint, closestPoint) / 2.0f;
+  *r_radius = math::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 70e3a99e57a..eac3faa6d15 100644
--- a/source/blender/blenlib/tests/BLI_delaunay_2d_test.cc
+++ b/source/blender/blenlib/tests/BLI_delaunay_2d_test.cc
@@ -21,10 +21,9 @@ 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_mpq2.hh"
+#include "BLI_math_vec_mpq_types.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
new file mode 100644
index 00000000000..8aa1f90fde2
--- /dev/null
+++ b/source/blender/blenlib/tests/BLI_math_vec_types_test.cc
@@ -0,0 +1,149 @@
+/* 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 207f310d902..74e54151a06 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_float3.hh"
+#include "BLI_math_vec_types.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 d759f0c3be4..2b8fb3dbea4 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 68111fb8eb1..d2d76593129 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"