diff --git a/src/Magnum/Math/Test/FunctionsTest.cpp b/src/Magnum/Math/Test/FunctionsTest.cpp
index e7834c7bb..114068016 100644
--- a/src/Magnum/Math/Test/FunctionsTest.cpp
+++ b/src/Magnum/Math/Test/FunctionsTest.cpp
@@ -63,6 +63,8 @@ struct FunctionsTest: TestSuite::Tester {
     void sqrt();
     void sqrtInverted();
     void lerp();
+    void lerpLimits();
+    void lerpInfinity();
     void lerpBool();
     void lerpInverted();
     void select();
@@ -128,6 +130,8 @@ FunctionsTest::FunctionsTest() {
               &FunctionsTest::sqrt,
               &FunctionsTest::sqrtInverted,
               &FunctionsTest::lerp,
+              &FunctionsTest::lerpLimits,
+              &FunctionsTest::lerpInfinity,
               &FunctionsTest::lerpBool,
               &FunctionsTest::lerpInverted,
               &FunctionsTest::select,
@@ -403,6 +407,44 @@ void FunctionsTest::lerp() {
     CORRADE_COMPARE(Math::lerp(2.0_usec, 5.0_usec, 0.5f), 3.5_usec);
 }
 
+template<class T, class U> T lerpOptimized(const T& a, const T& b, U t) {
+    /* One multiplication and two additions, while `T((U(1) - t)*a + t*b)` is
+       two multiplications, addition and subtraction. Doesn't correctly
+       preserve boundary values. */
+    return t*(b - a) + a;
+}
+
+void FunctionsTest::lerpLimits() {
+    CORRADE_COMPARE(Math::lerp(1.0e10f, 1.0e-5f, 0.0f), 1.0e10f);
+    CORRADE_COMPARE(Math::lerp(1.0e10f, 1.0e-5f, 1.0f), 1.0e-5f);
+    CORRADE_COMPARE(Math::lerp(1.0e-5f, 1.0e10f, 0.0f), 1.0e-5f);
+    CORRADE_COMPARE(Math::lerp(1.0e-5f, 1.0e10f, 1.0f), 1.0e10f);
+
+    CORRADE_COMPARE(lerpOptimized(1.0e10f, 1.0e-5f, 0.0f), 1.0e10f);
+    {
+        CORRADE_EXPECT_FAIL("\"Optimized\" version of a lerp doesn't correctly preserve boundary values with wildly different magnitudes.");
+        CORRADE_COMPARE(lerpOptimized(1.0e10f, 1.0e-5f, 1.0f), 1.0e-5f);
+    }
+}
+
+void FunctionsTest::lerpInfinity() {
+    CORRADE_COMPARE(Math::lerp(Constants::inf(), 0.0f, 0.0f), Constants::inf());
+    CORRADE_COMPARE(Math::lerp(0.0f, Constants::inf(), 1.0f), Constants::inf());
+    {
+        CORRADE_EXPECT_FAIL("Lerp with infinity doesn't correctly preserve the other boundary value.");
+        CORRADE_COMPARE(Math::lerp(Constants::inf(), 0.0f, 1.0f), 0.0f);
+        CORRADE_COMPARE(Math::lerp(0.0f, Constants::inf(), 0.0f), 0.0f);
+    }
+
+    CORRADE_COMPARE(lerpOptimized(0.0f, Constants::inf(), 1.0f), Constants::inf());
+    {
+        CORRADE_EXPECT_FAIL("\"Optimized\" version of a lerp doesn't correctly preserve boundary values if an infinity is present.");
+        CORRADE_COMPARE(lerpOptimized(Constants::inf(), 0.0f, 0.0f), Constants::inf());
+        CORRADE_COMPARE(lerpOptimized(Constants::inf(), 0.0f, 1.0f), 0.0f);
+        CORRADE_COMPARE(lerpOptimized(0.0f, Constants::inf(), 0.0f), 0.0f);
+    }
+}
+
 void FunctionsTest::lerpBool() {
     /* Scalar interpolation phase */
     CORRADE_COMPARE(Math::lerp(Vector3i{1, 2, 3}, Vector3i{5, 6, 7}, true), (Vector3i{5, 6, 7}));
diff --git a/src/Magnum/Math/Vector.h b/src/Magnum/Math/Vector.h
index 8aee3ba35..21aaa3697 100644
--- a/src/Magnum/Math/Vector.h
+++ b/src/Magnum/Math/Vector.h
@@ -74,6 +74,9 @@ namespace Implementation {
     template<std::size_t, class, class> struct VectorConverter;
     /* Needed by DualQuaternion and Functions.h (to avoid dependency between them) */
     template<class T, class U> T lerp(const T& a, const T& b, U t) {
+        /* While `t*(b - a) + a` is one ALU op less, the following is
+           guaranteed to correctly preserves exact boundary values with t being
+           0 or 1. See FunctionsTest::lerpLimits() for details. */
         return T((U(1) - t)*a + t*b);
     }