diff --git a/dlib/geometry/plane.d b/dlib/geometry/plane.d
index 9eedfe55..27a56bf8 100644
--- a/dlib/geometry/plane.d
+++ b/dlib/geometry/plane.d
@@ -91,7 +91,7 @@ struct Plane
 
     void normalize()
     {
-        float len = sqrt(x * x + y * y + z * z);
+        float len = hypot3(x, y, z);
         x /= len;
         y /= len;
         z /= len;
@@ -101,7 +101,7 @@ struct Plane
     Plane normalized()
     {
         Plane res;
-        float len = sqrt(x * x + y * y + z * z);
+        float len = hypot3(x, y, z);
         return Plane(x / len, y / len, z / len, d / len);
     }
 
diff --git a/dlib/geometry/support.d b/dlib/geometry/support.d
index 1d44d740..32475234 100644
--- a/dlib/geometry/support.d
+++ b/dlib/geometry/support.d
@@ -64,7 +64,7 @@ Vector3f subBox(Vector3f dir, Vector3f halfSize)
 Vector3f supCylinder(Vector3f dir, float radius, float height)
 {
     Vector3f result;
-    float sigma = sqrt((dir.x * dir.x + dir.z * dir.z));
+    float sigma = hypot(dir.x, dir.z);
 
     if (sigma > 0.0f)
     {
@@ -91,7 +91,7 @@ Vector3f supCone(Vector3f dir, float radius, float height)
     len = sqrt(len);
     float half_h = height * 0.5;
 
-    float sin_a = radius / sqrt(radius * radius + 4.0f * half_h * half_h);
+    float sin_a = radius / hypot(radius, 2.0f * half_h);
 
     if (dir[2] > len * sin_a)
         return Vector3f(0.0f, 0.0f, half_h);
diff --git a/dlib/image/filters/lens.d b/dlib/image/filters/lens.d
index b7db8fe7..26d4eb05 100644
--- a/dlib/image/filters/lens.d
+++ b/dlib/image/filters/lens.d
@@ -55,7 +55,7 @@ SuperImage lensDistortion(
     float halfWidth = cast(float)img.width / 2.0f;
     float halfHeight = cast(float)img.height / 2.0f;
 
-    float correctionRadius = sqrt(cast(float)(img.width ^^ 2 + img.height ^^ 2)) / strength;
+    float correctionRadius = hypot(cast(float)img.width, img.height) / strength;
 
     foreach(y; 0..img.height)
     foreach(x; 0..img.width)
@@ -63,7 +63,7 @@ SuperImage lensDistortion(
         float newX = x - halfWidth;
         float newY = y - halfHeight;
 
-        float distance = sqrt(newX ^^ 2 + newY ^^ 2);
+        float distance = hypot(newX, newY);
         float r = distance / correctionRadius;
 
         float theta;
diff --git a/dlib/math/complex.d b/dlib/math/complex.d
index 7ed63dc4..6649ae60 100644
--- a/dlib/math/complex.d
+++ b/dlib/math/complex.d
@@ -146,7 +146,7 @@ struct Complex(T)
 
     T magnitude()
     {
-        return sqrt(re * re + im * im);
+        return hypot(re, im);
     }
 
     T norm()
@@ -184,7 +184,7 @@ unittest
 /// Complex abs
 T abs(T)(Complex!T x)
 {
-    return sqrt(x.re * x.re + x.im * x.im);
+    return hypot(x.re, x.im);
 }
 
 ///
diff --git a/dlib/math/utils.d b/dlib/math/utils.d
index db9f11b3..6fef591a 100644
--- a/dlib/math/utils.d
+++ b/dlib/math/utils.d
@@ -146,6 +146,23 @@ unittest
     assert(min3(3, 2, 1) == 1);
 }
 
+/**
+ * Calculates the length of the hypotenuse for triangle defined in 3D
+ */
+T hypot3(T) (T x, T y, T z) nothrow
+{
+    // std.math.algebraic also provides hypot() version with 3 arguments
+    // but internally for now it calls not fast mlib's hypot()
+    return hypot(hypot(x, y), z);
+}
+
+///
+unittest
+{
+    const len = hypot3!double(-2, 6, 3);
+    assert(isConsiderZero(len - 7));
+}
+
 /**
  * Limit to given range
  */
diff --git a/dlib/math/vector.d b/dlib/math/vector.d
index 352dc3e8..ca18ca4c 100644
--- a/dlib/math/vector.d
+++ b/dlib/math/vector.d
@@ -490,18 +490,24 @@ struct Vector(T, int size)
         do
         {
             static if (isFloatingPoint!T)
+                alias CastTo = T;
+            else
+                alias CastTo = float;
+
+            static if (size == 2)
             {
-                T t = 0;
-                foreach (component; arrayof)
-                    t += component * component;
-                return sqrt(t);
+                return cast(T)hypot(cast(CastTo)x, cast(CastTo)y);
+            }
+            else static if (size == 3)
+            {
+                return cast(T)hypot3(cast(CastTo)x, cast(CastTo)y, cast(CastTo)z);
             }
             else
             {
                 T t = 0;
                 foreach (component; arrayof)
                     t += component * component;
-                return cast(T)sqrt(cast(float)t);
+                return cast(T)sqrt(cast(CastTo)t);
             }
         }
 
@@ -1087,7 +1093,17 @@ do
 {
     T dx = a.x - b.x;
     T dy = a.y - b.y;
-    return sqrt((dx * dx) + (dy * dy));
+    return hypot(dx, dy);
+}
+
+///
+unittest
+{
+    const a = Vector2f(4, 0);
+    const b = Vector2f(1, 4);
+    const l = distance(a, b);
+
+    assert(isConsiderZero(l - 5));
 }
 
 /**
@@ -1110,7 +1126,7 @@ do
     T dx = a.x - b.x;
     T dy = a.y - b.y;
     T dz = a.z - b.z;
-    return sqrt((dx * dx) + (dy * dy) + (dz * dz));
+    return hypot3(dx, dy, dz);
 }
 
 /**