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

30 files changed, 1396 insertions, 1822 deletions

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_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_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_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_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_utildefines.h b/source/blender/blenlib/BLI_utildefines.h
index 35d4158de59..c4b31810669 100644
--- a/source/blender/blenlib/BLI_utildefines.h
+++ b/source/blender/blenlib/BLI_utildefines.h
@@ -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/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"